今日の内容は、McElreath, R. (2016) Statistical Rethinking. CRC Press に基づく。

準備

これまでは、グリッド近似によってベイズ推定を行ってきたが、パラメタの数が増えるとグリッド近似の実行が難しくなる。そこで、これからは Stan を使ってベイズ推定を行う（詳しくはリンク先を参照）。Stan をR上で利用するために、RStan をインストールする。

install.packages('rstan', dependencies=TRUE)

次に、Stanモデルを直接書かなくても推定を実行できるようするため、Richard McElreath の rethinking パッケージをインストールする。このパッケージはCRANからではなく、github からインストールするので、devtools::install_github() を使う。（devtools パッケージがインストール済みでない場合は、install.packages('devtools', dep = TRUE) でまずインストールする。）

devtools::install_github('rmcelreath/rethinking')

最後に、いつも通り今回使うパッケージを読み込む。

library('MASS')
library('rgl')
library('dplyr')
library('ggplot2')
library('rethinking')

線形回帰モデル

統計分析のために、統計モデルを作らなければならない 回帰モデルには、以下の要素が必要である。

  1. 結果変数（応答変数）
  2. 尤度関数：結果変数がどのように生成されるか
  3. 予測変数（説明変数）
  4. 尤度関数と予測変数の関係
  5. 推定する母数（パラメタ）の事前確率

線形回帰の場合、これらの要素を元に以下のようなモデルを作る。

\[y_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Cauchy}^{+}(0, 1)\]

ただし、\(y\)は結果変数、\(x\)は予測変数、\(\alpha\), \(\beta\), \(\sigma\) が推定するパラメタである。

地球儀問題のモデル

第3回の講義で検討した地球儀問題を統計モデルで表現すると、次のようになる。

\[w \sim \mbox{Binomial}(n, \theta)\] \[\theta \sim \mbox{Uniform}(0, 1)\]

ただし、\(w\) はトスの結果「水」が得られた回数、\(n\) はトスの回数、\(\theta\) は地球儀上で水面がしめる割合（推定するパラメタ）である。

モデルのうち、「\(\sim\)」で繋がれた部分は、「確率的 (stochastic)」関係を表している。つまり、\(n\) と\(\theta\) の値がわかっても、\(w\) がどの値をとるかは決まらない。わかるのは、\(w\) が特定の二項分布に従うということだけである（よって、期待値、分散などはわかる）。

このモデルは、「\(w\) は、サイズ\(n\)で成功確率 \(\theta\) の二項分布にしたがっており」、「\(\theta\) は最小値0で最大値1の一様分布にしたがっている」ことを表現している。 このうち、第1行が尤度を、第2行が事前分布を表している。このモデルに予測変数はない。

尤度と事前分布がモデル化されているので、事後分布は以下の式で求められる。 \[p(\theta | w, n) = \frac{\mbox{Binomial}(w|n, \theta) \mbox{Uniform}(\theta|0, 1)}{\int \mbox{Binomial}(w|n, \theta) \mbox{Uniform}(\theta|0, 1)d\theta}\]

\(w=6\), \(n=9\) だとすると、この事後分布は以下のように推定できる。

w <- 6
n <- 9
theta_grid <- seq(from = 0, to = 1, length.out = 1000)
posterior <- dbinom(w, n, prob = theta_grid) * dunif(theta_grid, min = 0, max = 1)
posterior <- posterior / sum(posterior * 1/1000)  # normalizing
plot(theta_grid, posterior, type = 'l', xlab = expression(theta), col = 'blue')

正規分布モデル：身長の推定

線形回帰では、結果変数の尤度を正規分布によってモデル化する。 したがって、線形回帰を理解するためには、正規分布モデルの理解が不可欠である。 まずは、正規分布モデルを理解しよう。

データ

例として、rethinking パッケージに含まれる Howell1 というデータセットを使う。

data(Howell1)
myd <- Howell1
str(myd)
## 'data.frame':    544 obs. of  4 variables:
##  $ height: num  152 140 137 157 145 ...
##  $ weight: num  47.8 36.5 31.9 53 41.3 ...
##  $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
##  $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
summary(myd)
##      height           weight            age             male       
##  Min.   : 53.98   Min.   : 4.252   Min.   : 0.00   Min.   :0.0000  
##  1st Qu.:125.09   1st Qu.:22.008   1st Qu.:12.00   1st Qu.:0.0000  
##  Median :148.59   Median :40.058   Median :27.00   Median :0.0000  
##  Mean   :138.26   Mean   :35.611   Mean   :29.34   Mean   :0.4724  
##  3rd Qu.:157.48   3rd Qu.:47.209   3rd Qu.:43.00   3rd Qu.:1.0000  
##  Max.   :179.07   Max.   :62.993   Max.   :88.00   Max.   :1.0000

このデータセットには、身長、体重、年齢、性別（男性ダミー）の4変数が含まれている。 ここでは、身長 (height) を結果変数として利用する。

まず、結果変数の分布を確認してみよう。

hist_ht <- ggplot(myd, aes(x = height, y = ..density..)) +
    geom_histogram(binwidth = 2.5, color = 'black')
print(hist_ht)

身長が160cm付近に山があるが、非常に低い身長もたくさんあることがわかる。 低身長の1つの理由は、未成年のデータが含まれていることである。 そこで、データを18歳以上に制限しよう。

myd2 <- myd %>% filter(age >= 18)
hist_ht %+% myd2

先ほどのヒストグラムと比較すると、左右のバランスがよくなった。

統計モデル

次に、身長を結果変数とする統計モデルを作る。 上のヒストグラムから、身長の分布はなんとなく正規分布に近く見えるので、身長 (\(h\)) は正規分布から生成されていると仮定する。 正規分布のパラメタは、平均 \(\mu\) と標準偏差 \(\sigma\) で表すことにする。 すると、尤度は \[h_i \sim \mbox{Normal}(\mu, \sigma)\] と表すことができる。

ここで推定するパラメタは\(\mu\) と\(\sigma\) の2つである。 パラメタが2つあるので、事前確率は\(p(\mu, \sigma)\) という同時確率になる。 しかし、それぞれのパラメタの事前分布が独立であると仮定すれば、\(p(\mu, \sigma) = p(\mu)p(\sigma)\)　と書けるので、これを使ってモデルを作ることにする。

よって、身長を結果変数とする（予測変数のない）統計モデルは、次のようになる。 \[h_i \sim \mbox{Normal}(\mu, \sigma)\] \[\mu \sim \mbox{Normal}(168, 20)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

ここで、\(\mu\) の事前分布は168cm を中心とするばらつきの大きな正規分布である。 この事前分布によって、身長の平均値は95％の確率で128cm から208cm の間にあるだろうという事前の知識をモデルに組み込んでいる。また、正規分布を使うことによって、128cmや208cm よりは168cm のほうが平均値になりやすそうだということも考慮されている（中心を168cmにしたのは恣意的な選択）。

この事前分布を図示する。

curve(dnorm(x, mean = 168, sd = 20), from = 100, to = 240,
      xlab = expression(mu), ylab = 'prior density', col = 'red')

このように、この事前分布は値を狭い範囲に限定していない。

同様に、\(\sigma\) の事前分布を図示する。

curve(dunif(x, min = 0, max = 50), from = -10, to = 60,
      xlab = expression(sigma), ylab = 'prior density', col = 'red')

\(\sigma\) については事前の知識がないので、特定の値の密度を高くせず、範囲も広くとってある。\(\sigma > 0\) なので下限の設定は簡単だが、上限はよくわからないので少し大きめに設定してある。

これらの事前分布を使って、身長の事前分布（つまり、データを観察する前に身長がどのように分布していると考えるか）をシミュレートしてみよう。

n <- 10000
sample_mu <- rnorm(n, mean = 168, sd = 20)
sample_sigma <- runif(n, min = 0, max = 50)
prior_ht <- rnorm(n, mean = sample_mu, sd = sample_sigma)
dens(prior_ht, xlab = 'height (cm)', ylab = 'p(height)', col = 'red')

このように、なんとなく正規分布のような事前分布が得られる。

正規分布モデルのグリッド近似

グリッド近似を使い、このモデルのパラメタを推定してみよう。

mu_grid <- seq(from = 140, to = 200, length.out = 500)
sigma_grid <- seq(from = 4, to = 10, length.out = 500)
post <- expand.grid(mu = mu_grid, sigma = sigma_grid)
post$ll <- sapply(1:nrow(post), function(i) {
    sum(dnorm(myd2$height, mean = post$mu[i], sd = post$sigma[i], log = TRUE))})
post <- post %>%
    mutate(prod = ll + 
               dnorm(mu, 168, 20, log = TRUE) + 
               dunif(sigma, 0, 50, log = TRUE),
           prob = exp(prod - max(prod)))

同時事後分布を等高線図にする。

cont_ht <- ggplot(post, aes(mu, sigma, z = prob)) +
    geom_contour()
print(cont_ht)

3次元図（webでは非表示）。

plot3d(post[c('mu', 'sigma', 'prob')], col = 'skyblue',
       xlim = c(153, 156), ylim = c(6, 9))

事後分布からのサンプリング

事後分布をより詳しく調べるため、事後分布から無作為抽出を行う。 パラメタが2つあるので、\(\mu\) と\(\sigma\) をペアで選ぶ必要がある。 そこでまず、事後確率に基づいてグリッドから10,000行をランダムに復元抽出し、抽出された行にある\(\mu\)と\(\sigma\) をサンプルとして利用する。

sample_rows <- sample(1:nrow(post), size = 10000,
                      replace = TRUE, prob = post$prob)
sample_mu <- post$mu[sample_rows]
sample_sigma <- post$sigma[sample_rows]
df <- data.frame(mu = sample_mu, sigma = sample_sigma)

サンプルを散布図にしてみよう。

plt <- ggplot(df, aes(x = mu, y = sigma)) +
    geom_point(alpha = 1/5, color = 'royalblue')
print(plt)

各パラメタの周辺分布を図示する。

mgnl_mu <- ggplot(df, aes(mu)) + geom_density() +
    ylab('marginal posterior density')
mgnl_sigma <- ggplot(df, aes(sigma)) + geom_density() +
    ylab('marginal posterior density')
print(mgnl_mu)

print(mgnl_sigma)

どちらの分布も正規分布に近い。ただし、\(\sigma\)の周辺事後分布のほうが、右に歪んでいる（右裾が長い）ことがわかる。

2つのパラメタの77%HPDI（最高事後密度区間）を求めてみよう。

HPDI(sample_mu, prob = .77)
##    |0.77    0.77| 
## 153.9479 154.9098
HPDI(sample_sigma, prob = .77)
##    |0.77    0.77| 
## 7.426854 8.124248

この統計モデルとデータによれば、\(\mu\)が154cmから155cmの間にある確率と\(\sigma\)が7.4から8.1の間にある確率はどちらも約77％であることがわかる。

map を使ってパラメタを推定する

同じモデルをパッケージを使って推定してみよう。Stanモデルの書き方がわかるなら、rstan::stan() で推定すればよいが、ここでは rethinking::map() を使う。map はMAP (maximum a posteriori) を推定するための関数である。

この関数を使うために、まず alist() を使って統計モデルを記述する。 ここで使う統計モデルは、先ほどと同じ \[h_i \sim \mbox{Normal}(\mu, \sigma)\] \[\mu \sim \mbox{Normal}(168, 20)\] \[\sigma \sim \mbox{Uniform}(0, 50)\] である。これをRで記述する。

model1 <- alist(
    height ~ dnorm(mu, sigma),
    mu ~ dnorm(168, 20),
    sigma ~ dunif(0, 50)
)

このモデルを使い、map() でパラメタを推定する。

fit1 <- map(model1, data = myd2)

推定結果は rethinking::precis() で確認する。

precis(fit1)
##         Mean StdDev   5.5%  94.5%
## mu    154.60   0.41 153.94 155.26
## sigma   7.73   0.29   7.27   8.20

map() は正規分布によって近似された推定値を返す。したがって、この結果は、\(\mu\) の周辺事後分布は Normal(155, 0.4) であり、\(\sigma\) の周辺事後分布は Normal(7.7, 0.3) と推定されたことを示している。デフォルトでは89パーセンタイルの下限と上限が示されているが、77パーセンタイルに変更したいときは次のようにする。

precis(fit1, prob = .77)
##         Mean StdDev  11.5%  88.5%
## mu    154.60   0.41 154.11 155.10
## sigma   7.73   0.29   7.38   8.08

上で求めた77%HPDI と比較してみよう！

これまで使ってきたモデルでは\(\mu\) の事前分布として \(\mbox{Normal}(168, 20)\) を使ってきた。この事前分布では\(\mu\) が取り得る値の範囲が広い。したがって、この事前分布は「弱い」事前分布であると考えられる。 これに対し、\(\mu \sim \mbox{Normal}(168, 0.1)\) とすると、\(\mu\)が取り得る値が著しく制限されるので、これは「強い」事前分布である。この事前分布を使って推定を行ってみよう。

model2 <- alist(
    height ~ dnorm(mu, sigma),
    mu ~ dnorm(168, 0.1),
    sigma ~ dunif(0, 50)
)
fit2 <- map(model2, data = myd2)
precis(fit2, prob = .77)
##        Mean StdDev  11.5%  88.5%
## mu    167.8   0.10 167.68 167.92
## sigma  15.3   0.58  14.60  16.00

推定結果は、事前分布の影響を強く受けている。 また、\(\sigma\) の推定値が大きくなっている。 これはなぜだろうか？

map で推定した後のサンプリング

map で推定した場合も、事後分布について考えるにはサンプルがあったほうが便利である。 そこで、map推定後にサンプルを抽出する方法を考える。

現在考えているモデルには\(\mu\) と\(\sigma\) という2つのパラメタがあり、事後分布はこれらのパラメタの同時分布である。したがって、\(\mu\) と\(\sigma\) は同時に（ペアで）抽出する必要がある。どちらも正規分布で近似されているので、多次元正規分布 (multivariate normal distribution) からサンプリングを行えばばよい。Rでは、MASS:mvrnorm() を使えばよい。 mvrnorm() では、平均値ベクトルと分散共分散行列 (variance-covariance matrix) を指定する必要がある。 平均値には推定された平均値を使い、分散共分散行列は vcov()で求められるので、次のようにサンプルが取れる。

post <- mvrnorm(10000, mu = coef(fit1), Sigma = vcov(fit1))
head(post)
##            mu    sigma
## [1,] 154.9642 7.795044
## [2,] 154.9860 7.574889
## [3,] 154.2664 7.383163
## [4,] 153.9034 7.460735
## [5,] 154.8965 7.635642
## [6,] 154.8235 8.286587

このサンプルを使って、周辺事後分布を図示してみよう。

post <- as.data.frame(post)
mgnl_mu <- ggplot(post, aes(mu)) + 
    geom_histogram(aes(y=..density..), color = 'black')
mgnl_sigma <- ggplot(post, aes(sigma)) + 
    geom_histogram(aes(y=..density..), color = 'black')
print(mgnl_mu)

print(mgnl_sigma)

77%HPDI。

HPDI(post$mu, prob = .77)
##    |0.77    0.77| 
## 154.1041 155.0910
HPDI(post$sigma, prob = .77)
##    |0.77    0.77| 
## 7.385838 8.079943

線形回帰モデル

次に、線形回帰モデルを理解しよう。 上の正規分布モデルとの違いは、予測変数（説明変数）が統計モデルに登場するというところにある。

ここでは、例として体重 (weight [kg]) を予測変数とする。このとき、私たちは「身長を体重に回帰する (regress height on weight)」という。

統計モデル

このとき、統計モデルは次のように書ける。 \[h_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 10)\] \[\beta \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

ただし、\(x\) は体重である。 このモデルで、「=」で繋がれた部分は、「決定的」関係を示す。 \(\mu_i = \alpha + \beta x_i\) は尤度と予測変数の関係を表している。この関係が線形なので、このモデルは「線形回帰 (linear regression)」と呼ばれる。また、尤度が正規分布なので、正規線形モデル (normal-linear model) と呼ばれることもある。

推定

このようにモデルかできれば、推定は先ほどと同様の方法で行える。

model3 <- alist(
    height ~ dnorm(mu, sigma),
    mu <- alpha + beta * weight,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit3 <- map(model3, data = myd2)

推定結果の解釈

推定結果を見てみよう。

precis(fit3, prob = .77)
##         Mean StdDev  11.5%  88.5%
## alpha 109.84   1.91 107.55 112.13
## beta    0.99   0.04   0.94   1.04
## sigma   5.11   0.19   4.87   5.34

まず、alpha の推定値である 113.9 は何を意味しているだろうか。これは、体重が0の人の\(\mu_i\) は113.9cmになるということを示している（そんな人はいないので、このままではまずい）。

次に、beta の推定値である0.9 は何を意味するだろうか。これは、体重が1kg 増えるごとに、\(\mu_i\) が0.9cm 高くなることを示している。私たちは線形の関係を想定しているので、体重が10kg 増えれば身長は9cm 高くなることが予測される。ただし、データの範囲外（例えば、体重が500kgの場合の身長）について予測を行うことは危険である。

最後に、sigma の推定値である5は、\(\mu\)の分布の幅を示している。正規分布するデータの95%はプラマイ2標準偏差内に収まるはずなので、身長の95％は推定された平均値からプラスマイナス10cm程度に収まると考えられる。

このモデルでは alpha が意味のない値になってしまった。 この例に限らず、切片の値は意味がないものになることが多い。 このような事態を避けるために、予測変数は中心化 (centering) したほうがよい。

中心化する基準になるものは色々あるが、最も単純なのは算術平均で中心化する方法である。体重を標本平均で中心化して推定し直してみよう。

myd2 <- myd2 %>% mutate(weight_c = weight - mean(weight)) 
model4 <- alist(
    height ~ dnorm(mu, sigma),
    mu <- alpha + beta * weight_c,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit4 <- map(model4, data = myd2)
precis(fit4, prob = .77)
##         Mean StdDev  11.5%  88.5%
## alpha 154.48   0.27 154.16 154.81
## beta    0.91   0.04   0.85   0.96
## sigma   5.07   0.19   4.84   5.30

beta とsigma の推定値は先ほどと（ほぼ）同じだが、alpha の推定値は154.6 になった。これは、体重が平均値のときの身長の平均値であり、データの平均値に一致する。

mean(myd2$height)
## [1] 154.5971

推定結果を図示する

回帰直線の表示

結果を図示する方法は色々あるが、とりあえずデータ上に回帰直線を乗せてみよう。

scat <- ggplot(myd2, aes(weight, height)) +
    geom_point(color = 'royalblue') +
    geom_abline(intercept = coef(fit3)['alpha'],
                slope = coef(fit3)['beta'])
print(scat)

この直線のおかげで、体重と身長の関係がわかりやすくなった。

推定の不確実性

しかし、上の図からは推定の不確実性 (uncertainty) を読み取ることができない。そこで、事後分布からのサンプルを利用して図を作り直す。

post <- mvrnorm(10000, mu = coef(fit3), Sigma = vcov(fit3)) %>%
    as.data.frame()

この中からalpha とbetaのペアをランダムに100組選び、それを利用して推定の不確実性を図示する。

post100 <- post %>% sample_n(100)
scat <- ggplot(myd2, aes(weight, height)) +
    geom_abline(intercept = post100$alpha,
                slope = post100$beta,
                color = 'gray') +
    geom_point(color = 'royalblue') +
    geom_abline(intercept = coef(fit3)['alpha'],
                slope = coef(fit3)['beta'])
print(scat)

サンプルサイズの変化が推定の不確実性に与える変化を確認するために、データの一部だけを使って推定を行いその結果を図示してみよう。

まず、データセットからランダムに10ユニットだけ取り出し、線形回帰モデルでパラメタを推定する。

set.seed(2016-07-06)
myd_sml <- myd %>% sample_n(10)
fit_sml <- map(model3, data = myd_sml)
post_sml <- mvrnorm(100, mu = coef(fit_sml), Sigma = vcov(fit_sml)) %>%
    as.data.frame()
scat <- ggplot(myd_sml, aes(weight, height)) +
    geom_abline(intercept = post_sml$alpha,
                slope = post_sml$beta,
                color = 'gray') +
    geom_point(color = 'royalblue') +
    geom_abline(intercept = coef(fit_sml)['alpha'],
                slope = coef(fit_sml)['beta'])
print(scat)

同じことを、50ユニットで行う。

myd_sml <- myd %>% sample_n(50)
fit_sml <- map(model3, data = myd_sml)
post_sml <- mvrnorm(100, mu = coef(fit_sml), Sigma = vcov(fit_sml)) %>%
    as.data.frame()
scat <- ggplot(myd_sml, aes(weight, height)) +
    geom_abline(intercept = post_sml$alpha,
                slope = post_sml$beta,
                color = 'gray') +
    geom_point(color = 'royalblue') +
    geom_abline(intercept = coef(fit_sml)['alpha'],
                slope = coef(fit_sml)['beta'])
print(scat)

同様に、200ユニットで行う。

myd_sml <- myd %>% sample_n(200)
fit_sml <- map(model3, data = myd_sml)
post_sml <- mvrnorm(100, mu = coef(fit_sml), Sigma = vcov(fit_sml)) %>%
    as.data.frame()
scat <- ggplot(myd_sml, aes(weight, height)) +
    geom_abline(intercept = post_sml$alpha,
                slope = post_sml$beta,
                color = 'gray') +
    geom_point(color = 'royalblue') +
    geom_abline(intercept = coef(fit_sml)['alpha'],
                slope = coef(fit_sml)['beta'])
print(scat)

このように、サンプルサイズが大きくなると、推定の不確実性（回帰直線のぐらつき）が小さくなることがわかる。

推定区間の表示

回帰直線は、条件付き期待値の集合であると考えることができる。予測変数を\(X\)、結果変数を\(Y\) とすると、\(X=x\) で回帰直線は \(\mathbb{E}[Y|X=x]\) という値をとる。しかし、推定された回帰直線はばらつく（推定の不確実性がある）ので、各\(x\)において、期待値の推定値にはばらつきがある。

例えば、weight が 45kg のときの身長の平均値は、シミュレートしたサンプルを使って次のように求められる。

mu_w45 <- post$alpha + post$beta * 45

これを図示する。

dens(mu_w45, col = 'royalblue',
     xlab = expression(paste(mu, '|weight=45')))

この \(\mu\) は、統計モデルでは、\(\mu_i = \alpha + \beta x_i\) と表現されている。 ここで、\(\alpha\) と\(\beta\) がパラメタであり、これらは分布する。したがって、\(\alpha\) と\(\beta\) によって値が決まる\(\mu\) も分布する。よって、\(\mu\) についてもHPDI を考えることができる。 weight が45kg のときの \(\mu\) の77％HPDI は

HPDI(mu_w45, prob = .77)
##    |0.77    0.77| 
## 154.1852 154.8364

である。つまり、体重が45kgの人の平均身長の分布の中央部分77％は、154.2cmから154.8cmの間にある。

これをweightのすべての値について計算すれば、区間推定を図示することができる。 データセットにあるweightの範囲は

with(myd2, range(weight))
## [1] 31.07105 62.99259

なので、この範囲で\(\mu\)の区間を計算する。 まずは、シミュレートしたalpha とbeta のペアを使って、\(\mu\)を計算する。

weight_vec <- seq(min(myd2$weight), max(myd2$weight), length.out = 100)
mu_w <- sapply(weight_vec, function(x) post$alpha + post$beta * x)

これで、weight_vec の各要素（つまり、weight の各値）に対し、それぞれ10000個のmu が計算できた。 weight_vec の各要素について、muの平均値と93％HPDI を求める。

mu_mean <- apply(mu_w, 2, mean)
mu_HPDI <- apply(mu_w, 2, HPDI, prob = .93)  # 93% HPDI

これを、図示する。

df <- data.frame(weight = weight_vec, mu_mean,
                 lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
intrvl <- ggplot(myd2, aes(x = weight)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_point(aes(y = height), color = 'royalblue') +
    geom_line(data = df, aes(y = mu_mean))
print(intrvl)

この図で、黒い直線で示されているのが、muの推定値の平均値（回帰直線）である。回帰直線の周りにグレーで色付けされているのが、weightの各値に対応する 93%HPDI である。93％にした理由は特にない（95％と同様に恣意的な選択である。）

予測区間の表示

上で示したのは、各体重に対応する平均身長\(\mu\)の区間（推定の不確実）である。 今度は、身長そのものの予測値とその予測の不確実性を図に示そう。

私たちが利用している統計モデルでは、身長 \(h\) は次のように決まる。 \[h_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] したがって、\(h\) は\(\mu\)のばらつきだけでなく、\(\sigma\)の影響も受けることを考慮に入れる必要がある。つまり、\(h\)は、推定された\(\sigma\) を使って正規分布から無作為に抽出する必要がある。

# 各体重について、muの推定値とsigmaの推定値を使って height を無作為抽出
sim_height <- sapply(weight_vec, function(x) {
    rnorm(n = nrow(post), mean = post$alpha + post$beta*x, sd = post$sigma)
    })
# 各体重に対応する93% percentile interval を求める
height_pi <- apply(sim_height, 2, PI, prob = .93)

これを図示する。

df <- df %>%
    mutate(lb_pi = height_pi[1,],
           ub_pi = height_pi[2,])
pred_plot <- ggplot(myd2, aes(weight)) +
    geom_ribbon(data = df, aes(ymin = lb_pi, ymax = ub_pi), fill = 'gray') +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray55') +
    geom_point(aes(y = height), color = 'royalblue') +
    geom_line(data = df, aes(y = mu_mean))
print(pred_plot)

この図は、MAP推定値（直線）、\(\mu_i\) のHPDI（濃いグレー）、予測値\(\hat{h_i}\)のPI（淡いグレー）を同時に示している。予測値\(\hat{h_i}\)のPIは93％区間で示されているので、データポイントのうち約7パーセントはグレーで色付けされた領域の外側にあることが期待される。

最後に、線形回帰で通常もっとも注目される\(\beta\)の周辺事後分布を示す。

mgnl_beta <- ggplot(post, aes(beta, ..density..)) + 
    geom_histogram(color = 'black') +
    labs(x = expression(beta), y = expression(paste('p(', beta, '|h, w)')))
print(mgnl_beta)

回帰分析の結果として注目するのは\(beta\)の値（特に、統計的因果推論においては処置変数 (treatment variable) の係数のみ）であり、多くの場合はその効果が「0でない」ことをまず示したいので、次のような図で結果を報告するのが（一般的には）良いだろう。

result_beta <- mgnl_beta +
    xlim(0, max(post$beta)) +
    geom_vline(xintercept = 0, color = 'red') +
    labs(x = expression(beta), y = expression(paste('p(', beta, '|h, w)')))
print(result_beta)

ただし、この例の場合は1つ前のヒストグラムから\(\beta\)が0より大きくなる確率がほぼ1になることは明らかなので、\(\beta=0\)の線を無理に加える必要はないだろう。効果が0でないことを示したら、次にすべきことは効果量 (effect size) を明らかにすることである。今回の例では、\(\beta\) の事後分布は0.9付近に集中しているので、「体重が1kg 増加するごとに身長が0.9cm 伸びる」というのが推定される効果量である。有意かどうかよりも、効果量の大きさの方が重要なので、必ず実質的な意味づけとともに効果量を説明することが求められる。

いずれにせよ、表よりも図（ヒストグラムや散布図）のほうが結果がわかりやすいので、図を積極的に活用すべきである。



授業の内容に戻る