{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "UIp998OHnZSN"
},
"source": [
"##### Copyright 2018 The TensorFlow Probability Authors.\n",
"\n",
"Licensed under the Apache License, Version 2.0 (the \"License\");"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"id": "li5wNGR6naj0"
},
"outputs": [],
"source": [
"#@title Licensed under the Apache License, Version 2.0 (the \"License\"); { display-mode: \"form\" }\n",
"# you may not use this file except in compliance with the License.\n",
"# You may obtain a copy of the License at\n",
"#\n",
"# https://www.apache.org/licenses/LICENSE-2.0\n",
"#\n",
"# Unless required by applicable law or agreed to in writing, software\n",
"# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
"# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
"# See the License for the specific language governing permissions and\n",
"# limitations under the License."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ykgJW69K4TNL"
},
"source": [
"# {TF Probability、R、Stan} における線形混合効果回帰\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "epdzLgKQRfAI"
},
"source": [
"## 1 はじめに"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "bMstt69hR44D"
},
"source": [
"このコラボでは、線形混合効果回帰モデルを人気のあるトイデータセットに適合させます。R の `lme4`、Stan の混合効果パッケージ、および TensorFlow Probability (TFP) プリミティブを使用して、これを 3 回適合させます。そして、これらからほぼ同じ適合パラメータと事後分布を得られることを示します。\n",
"\n",
"主な結論として、TFP には HLM のようなモデルを適合させるために必要な一般的な要素があり、`lme4` や `rstanarm` などの他のソフトウェアパッケージと一致する結果を生成します。このコラボは、比較したパッケージの計算効率を正確に反映したものではありません。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "0axKjgZvRtL9"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import os\n",
"from six.moves import urllib\n",
"import numpy as np\n",
"import pandas as pd\n",
"import warnings\n",
"\n",
"from matplotlib import pyplot as plt\n",
"import seaborn as sns\n",
"\n",
"from IPython.core.pylabtools import figsize\n",
"figsize(11, 9)\n",
"\n",
"import tensorflow.compat.v1 as tf\n",
"import tensorflow_datasets as tfds\n",
"import tensorflow_probability as tfp"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "IFC9r-h0XlQ3"
},
"source": [
"## 2 階層線形モデル\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "wqS-HhvhDlno"
},
"source": [
"R、Stan、TFP を比較するために、[階層線形モデル](https://en.wikipedia.org/wiki/Multilevel_model) (HLM) を[ラドンデータセット](http://www.stat.columbia.edu/~gelman/arm/examples/radon/)に適合させます。このデータセットはゲルマンら (559 ページ、第 2 版、250 ページ、第 3 版) による[*ベイジアンデータ分析*](http://www.stat.columbia.edu/~gelman/book/)で有名になったものです。\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "fAD8am2a4TaY"
},
"source": [
"次の生成モデルを前提としています。\n",
"\n",
"$$\\begin{align*} \\text{for } & c=1\\ldots \\text{NumCounties}:\\ & \\beta_c \\sim \\text{Normal}\\left(\\text{loc}=0, \\text{scale}=\\sigma_C \\right) \\ \\text{for } & i=1\\ldots \\text{NumSamples}:\\ &\\eta_i = \\underbrace{\\omega_0 + \\omega_1 \\text{Floor}i}\\text{fixed effects} + \\underbrace{\\beta_{ \\text{County}i} \\log( \\text{UraniumPPM}{\\text{County}i}))}\\text{random effects} \\ &\\log(\\text{Radon}_i) \\sim \\text{Normal}(\\text{loc}=\\eta_i , \\text{scale}=\\sigma_N) \\end{align*}$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "5styKLl_MyWu"
},
"source": [
"R の `lme4` 「チルダ表記」では、このモデルは次と同等です。\n",
"\n",
"> `log_radon ~ 1 + floor + (0 + log_uranium_ppm | county)`\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "SpurhP2gP_T8"
},
"source": [
"${\\beta_c}_{c=1}^\\text{NumCounties}$ の事後分布 (証拠を条件とする) を使用して、$\\omega, \\sigma_C, \\sigma_N$ の MLE を見つけます。"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Nj6adNwgPTUP"
},
"source": [
"本質的に同じモデルですが、ランダム切片があります。*[付録 A](#scrollTo=tsXhZ4rtNUXL)* を参照してください。\n",
"\n",
"HLM のより一般的な仕様については、*[付録 B](#scrollTo=H0w7ofFvNsxi)* を参照してください。"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "LR0ZC0dE4MWb"
},
"source": [
"## 3 データマンジング"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "OdboSs9G3JlE"
},
"source": [
"このセクションでは、[radon データセット](http://www.stat.columbia.edu/~gelman/arm/examples/radon/)を取得し、想定されるモデルに準拠するように最小限の前処理を行います。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "4LjOBqLDV0IQ"
},
"outputs": [],
"source": [
"def load_and_preprocess_radon_dataset(state='MN'):\n",
" \"\"\"Preprocess Radon dataset as done in \"Bayesian Data Analysis\" book.\n",
" \n",
" We filter to Minnesota data (919 examples) and preprocess to obtain the\n",
" following features:\n",
" - `log_uranium_ppm`: Log of soil uranium measurements.\n",
" - `county`: Name of county in which the measurement was taken.\n",
" - `floor`: Floor of house (0 for basement, 1 for first floor) on which the\n",
" measurement was taken.\n",
"\n",
" The target variable is `log_radon`, the log of the Radon measurement in the\n",
" house.\n",
" \"\"\"\n",
" ds = tfds.load('radon', split='train')\n",
" radon_data = tfds.as_dataframe(ds)\n",
" radon_data.rename(lambda s: s[9:] if s.startswith('feat') else s, axis=1, inplace=True)\n",
" df = radon_data[radon_data.state==state.encode()].copy()\n",
"\n",
" # For any missing or invalid activity readings, we'll use a value of `0.1`.\n",
" df['radon'] = df.activity.apply(lambda x: x if x > 0. else 0.1)\n",
" # Make county names look nice. \n",
" df['county'] = df.county.apply(lambda s: s.decode()).str.strip().str.title()\n",
" # Remap categories to start from 0 and end at max(category).\n",
" county_name = sorted(df.county.unique())\n",
" df['county'] = df.county.astype(\n",
" pd.api.types.CategoricalDtype(categories=county_name)).cat.codes\n",
" county_name = list(map(str.strip, county_name))\n",
"\n",
" df['log_radon'] = df['radon'].apply(np.log)\n",
" df['log_uranium_ppm'] = df['Uppm'].apply(np.log)\n",
" df = df[['idnum', 'log_radon', 'floor', 'county', 'log_uranium_ppm']]\n",
"\n",
" return df, county_name"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "hJE3-eC0I-Lm"
},
"outputs": [],
"source": [
"radon, county_name = load_and_preprocess_radon_dataset()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "nV-IAEW2FIqX"
},
"outputs": [],
"source": [
"# We'll use the following directory to store our preprocessed dataset.\n",
"CACHE_DIR = os.path.join(os.sep, 'tmp', 'radon')\n",
"\n",
"# Save processed data. (So we can later read it in R.)\n",
"if not tf.gfile.Exists(CACHE_DIR):\n",
" tf.gfile.MakeDirs(CACHE_DIR)\n",
"with tf.gfile.Open(os.path.join(CACHE_DIR, 'radon.csv'), 'w') as f:\n",
" radon.to_csv(f, index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ubvf1vHenyCx"
},
"source": [
"### 3.1 データを調査する"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "r39MiQV2zUz0"
},
"source": [
"このセクションでは、提案されたモデルが合理的である理由を理解するために、`radon` データセットを調べます。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "GRCyjhSknu9z"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" log_radon | \n",
" floor | \n",
" county | \n",
" log_uranium_ppm | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0.788457 | \n",
" 1 | \n",
" 0 | \n",
" -0.689048 | \n",
"
\n",
" \n",
" | 1 | \n",
" 0.788457 | \n",
" 0 | \n",
" 0 | \n",
" -0.689048 | \n",
"
\n",
" \n",
" | 2 | \n",
" 1.064711 | \n",
" 0 | \n",
" 0 | \n",
" -0.689048 | \n",
"
\n",
" \n",
" | 3 | \n",
" 0.000000 | \n",
" 0 | \n",
" 0 | \n",
" -0.689048 | \n",
"
\n",
" \n",
" | 4 | \n",
" 1.131402 | \n",
" 0 | \n",
" 1 | \n",
" -0.847313 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" grpvar | \n",
" term | \n",
" grp | \n",
" condval | \n",
" condsd | \n",
"
\n",
" \n",
" \n",
" \n",
" | 1 | \n",
" county | \n",
" log_uranium_ppm | \n",
" 0 | \n",
" 0.667653 | \n",
" 0.465584 | \n",
"
\n",
" \n",
" | 2 | \n",
" county | \n",
" log_uranium_ppm | \n",
" 1 | \n",
" 0.697805 | \n",
" 0.123133 | \n",
"
\n",
" \n",
" | 3 | \n",
" county | \n",
" log_uranium_ppm | \n",
" 2 | \n",
" -0.010856 | \n",
" 0.847489 | \n",
"
\n",
" \n",
" | 4 | \n",
" county | \n",
" log_uranium_ppm | \n",
" 3 | \n",
" -0.068872 | \n",
" 0.422883 | \n",
"
\n",
" \n",
" | 5 | \n",
" county | \n",
" log_uranium_ppm | \n",
" 4 | \n",
" 0.036075 | \n",
" 0.825677 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" X.Intercept. | \n",
" floor | \n",
" b.log_uranium_ppm.county.0. | \n",
" b.log_uranium_ppm.county.1. | \n",
" b.log_uranium_ppm.county.2. | \n",
" b.log_uranium_ppm.county.3. | \n",
" b.log_uranium_ppm.county.4. | \n",
" b.log_uranium_ppm.county.5. | \n",
" b.log_uranium_ppm.county.6. | \n",
" b.log_uranium_ppm.county.7. | \n",
" ... | \n",
" b.log_uranium_ppm.county.76. | \n",
" b.log_uranium_ppm.county.77. | \n",
" b.log_uranium_ppm.county.78. | \n",
" b.log_uranium_ppm.county.79. | \n",
" b.log_uranium_ppm.county.80. | \n",
" b.log_uranium_ppm.county.81. | \n",
" b.log_uranium_ppm.county.82. | \n",
" b.log_uranium_ppm.county.83. | \n",
" b.log_uranium_ppm.county.84. | \n",
" sigma | \n",
"
\n",
" \n",
" \n",
" \n",
" | 1 | \n",
" 1.450331 | \n",
" -0.757375 | \n",
" 1.093099 | \n",
" 0.603985 | \n",
" -0.661680 | \n",
" -0.667052 | \n",
" 0.549178 | \n",
" -0.150145 | \n",
" 1.486173 | \n",
" 0.053578 | \n",
" ... | \n",
" 0.323107 | \n",
" 0.282121 | \n",
" -0.188251 | \n",
" 0.889385 | \n",
" 1.008348 | \n",
" 0.429430 | \n",
" 0.404690 | \n",
" 0.112424 | \n",
" -0.071612 | \n",
" 0.778706 | \n",
"
\n",
" \n",
" | 2 | \n",
" 1.453383 | \n",
" -0.659085 | \n",
" 0.560266 | \n",
" 0.580569 | \n",
" 0.419441 | \n",
" -0.154925 | \n",
" -1.263534 | \n",
" 0.166291 | \n",
" 1.534917 | \n",
" 1.141392 | \n",
" ... | \n",
" 0.306941 | \n",
" 0.331696 | \n",
" 0.614852 | \n",
" 0.153913 | \n",
" 1.253688 | \n",
" 0.246746 | \n",
" 0.849947 | \n",
" -1.399027 | \n",
" -1.697360 | \n",
" 0.754960 | \n",
"
\n",
" \n",
" | 3 | \n",
" 1.460351 | \n",
" -0.804513 | \n",
" 1.145694 | \n",
" 0.705337 | \n",
" 2.216364 | \n",
" -1.044031 | \n",
" 1.510397 | \n",
" -0.472675 | \n",
" 0.427448 | \n",
" 0.940869 | \n",
" ... | \n",
" 1.941637 | \n",
" 0.188311 | \n",
" -0.690044 | \n",
" 0.946370 | \n",
" 1.626093 | \n",
" 0.930756 | \n",
" 0.141125 | \n",
" 0.041392 | \n",
" -0.135538 | \n",
" 0.763605 | \n",
"
\n",
" \n",
" | 4 | \n",
" 1.438725 | \n",
" -0.565721 | \n",
" 0.204148 | \n",
" 0.810580 | \n",
" -1.982466 | \n",
" 1.044825 | \n",
" -1.339406 | \n",
" 0.643384 | \n",
" 2.379191 | \n",
" 0.295981 | \n",
" ... | \n",
" -0.832814 | \n",
" 0.246694 | \n",
" -0.801078 | \n",
" 0.072430 | \n",
" 0.108882 | \n",
" -0.203805 | \n",
" 0.443365 | \n",
" -0.901578 | \n",
" -0.063022 | \n",
" 0.750404 | \n",
"
\n",
" \n",
" | 5 | \n",
" 1.457569 | \n",
" -0.658910 | \n",
" 0.140554 | \n",
" 0.593056 | \n",
" -2.618940 | \n",
" 0.484812 | \n",
" -0.417834 | \n",
" 0.735558 | \n",
" 1.693997 | \n",
" 0.230496 | \n",
" ... | \n",
" -0.608215 | \n",
" -0.586245 | \n",
" 0.204547 | \n",
" 0.350195 | \n",
" 0.164781 | \n",
" -0.009945 | \n",
" 0.023533 | \n",
" -0.177017 | \n",
" -0.425871 | \n",
" 0.755898 | \n",
"
\n",
" \n",
"
\n",
"
5 rows × 88 columns
\n",
"
![](\"https://www.tensorflow.org/images/tf_logo_32px.png\"){kind=logo}
TensorFlow.org で表示![](\"https://www.tensorflow.org/images/colab_logo_32px.png\"){kind=logo}
Google Colab で実行![](\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\")
GitHubでソースを表示![](\"https://www.tensorflow.org/images/download_logo_32px.png\"){kind=logo}
ノートブックをダウンロード