毎回やる儀式

# 野球ファイルの読み込み
baseball <- read.csv("ball2017.csv",na.strings="*",fileEncoding = "UTF-8")
# 年収単位の調整(100万円に)
baseball$pay <- baseball$pay / 100
# M1のスコアファイルの読み込み
m1 <- read.csv("m1score.csv",na.strings=".",fileEncoding = "UTF-8")

はじめに

前回は一般化線形モデル(GLM)というお話をしました。今回はHLMです。 GのつぎはHかよ，というツッコミが聞こえて来そうですね。これは正式にはHierarichal Linear Modeling,ヒエラルキーのある線形モデル，という意味です。ヒエラルキーがあるというのは，とあるカテゴリーの中に複数のものが含まれているという構造のこと。例えば・・・

というのが「ヒエラルキーのあるデータ」構造。

そんなの気にしなくていいじゃない，と思うかもしれません。例えば最初の例だと，\(30 \times 3=90\)で分析したらいいじゃないか，と思うかもしれません。

では次のようなデータを見てください。

set.seed(12345)
N <- 20
# 独立変数を一様分布からつくる。ちょっと範囲が違う
X1 <- runif(N,1,10)
X2 <- runif(N,10,20)
X3 <- runif(N,20,30)
# 線形モデルを仮定。ただし，切片が各群で異なる
Yhat1 <- 0 - 2 * X1
Yhat2 <- 50 - 2 * X2
Yhat3 <- 100 - 2 * X3
# それに平均0,SD=10の誤差がついてデータが得られたとします
sig <- 10
Y1 <- Yhat1 + rnorm(N,0,sig)
Y2 <- Yhat2 + rnorm(N,0,sig)
Y3 <- Yhat3 + rnorm(N,0,sig)
# 組み合わせます
dat <- data.frame(list(class=factor(rep(1:3,each=N),labels=c("A","B","C")),X=c(X1,X2,X3),Y=c(Y1,Y2,Y3)))
# 全体の相関
cor(dat$X,dat$Y)
## [1] 0.8150185
# 散布図
library(ggplot2)
g <- ggplot(dat,aes(x=X,y=Y))+geom_point()+geom_smooth(method="lm")
g

# 全体の回帰分析
summary(lm(Y~X,data=dat))
## 
## Call:
## lm(formula = Y ~ X, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.734 -10.297  -0.281  12.418  36.989 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -19.8554     4.6115  -4.306 6.50e-05 ***
## X             2.8638     0.2673  10.712 2.28e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.15 on 58 degrees of freedom
## Multiple R-squared:  0.6643, Adjusted R-squared:  0.6585 
## F-statistic: 114.8 on 1 and 58 DF,  p-value: 2.28e-15

この散布図を見て，「綺麗な正の相関！」と思ったのではないでしょうか。回帰係数は2.6だし。でもちょっと待って，郡ごとに色分けしてみると・・・

g <- g + geom_point(aes(color=class))+geom_smooth(method="lm")
g

傾向は明らかに逆じゃありませんか？ そういえば，データを作った時は，\(Y=\beta_{0j} - 2 X\)だったので，回帰係数はマイナスのはずなのです。それがプラスになって出ちゃってる。群ごとに違うデータのはずなのに，まとめてやってしまうと，間違った結論に到達する可能性がある，というのがこれでお分りいただけると思います。

もちろん群ごとに分けても変わらないデータ構造というのもありますが，階層データはグループごとの違いの可能性を考えておいた方がいいでしょう。 こういう場合，適切な分析は次のように群ごとにするべきですね。

g <- ggplot(dat,aes(x=X,y=Y,color=class)) + geom_point()+geom_smooth(method="lm")
g

階層モデル

階層モデルは，群ごとに傾きや切片が異なることをモデル化するものです。 野球のデータを例にして考えてみましょう。 たとえば，どのチームに所属していても，ヒットをたくさん打つ選手ほど報酬が高い選手である，という傾向があると仮定します。 問題は，チームごとに「裕福かどうか」が違うので，ベースとなる基本給=切片が違うし，一本あたりの報酬増加率＝傾きも違う，と考えてみましょう。モデルで表現すると，

\[ Pay_{ij} = \beta_{0j} + \beta_{1j} Hit_i \]

です。ここで\(i\)は個人の番号，\(j\)は所属球団をさします。手持ちのデータのプロットを見てみましょう。

g <- ggplot(baseball,aes(x=AB,y=pay,color=team))+geom_point()+geom_smooth(method='lm',se=F)
g

まずデータセットを準備します。

# 必要な変数だけ抜き出す
base2 <- subset(baseball,select=c("team","AB","pay"))
# チーム名をチームを表す番号に書き換える
base2$team <- as.numeric(base2$team)
# データの一部を確認
head(base2)
##   team  AB pay
## 1    2 504 100
## 2    6 542 100
## 3    8 344 110
## 4   11 469 150
## 5    5 164  82
## 6    3 478 400

次に分析のコードを準備します。次のコードをコピペしてHLM1.stanという名前で保存してください。

data{
  int<lower=0> N;                 // プレイヤーの数
  int<lower=0> K;                 // チームの数
  real X[N];                      // 独立変数
  real Y[N];                      // 従属変数
  int<lower=0,upper=K> teamID[N]; // プレイヤーがどこのチームに属しているかを表す識別子
}

parameters{
  real beta0[K];          // 切片
  real beta1[K];          // 傾き
  real<lower=0> sig;      // 誤差分散
}

model{
  for(n in 1:N){
    Y[n] ~ normal(beta0[teamID[n]] + beta1[teamID[n]] * X[n] ,  sig);
  }
}

これを次のコードで呼び出して分析します。

# データの作成
datastan <- list(N=nrow(base2),K=max(base2$team),X=base2$AB,Y=base2$pay,teamID=base2$team)
# rstanの準備
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# stanモデルのコンパイル
model <- stan_model("HLM1.stan")
# モデルフィット
fit <- sampling(model,datastan)
# 結果出力
fit
## Inference for Stan model: HLM1.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##              mean se_mean     sd    2.5%     25%     50%     75%   97.5%
## beta0[1]     1.14    2.55 161.45 -313.90 -108.59    0.89  114.06  318.50
## beta0[2]   -32.29    2.32 146.43 -318.52 -130.04  -32.95   65.04  248.68
## beta0[3]    81.74    1.84 116.09 -143.98    5.68   80.10  155.47  319.50
## beta0[4]   -93.25    2.97 178.76 -432.79 -212.01  -91.27   23.49  263.00
## beta0[5]    62.32    2.03 128.16 -189.09  -24.28   62.35  145.14  314.52
## beta0[6]    79.40    2.22 140.31 -191.52  -14.14   82.40  169.91  357.44
## beta0[7]  -141.95    2.74 173.10 -482.32 -262.06 -142.84  -22.26  194.86
## beta0[8]    92.65    1.58 100.22 -113.12   26.13   95.79  159.39  290.52
## beta0[9]  -244.29    3.50 221.08 -678.77 -389.09 -245.70  -99.24  187.59
## beta0[10]  182.36    2.19 138.44  -81.77   85.10  182.63  274.74  455.64
## beta0[11]   55.48    2.60 164.46 -273.62  -52.08   58.01  164.81  374.03
## beta0[12]   63.51    1.49  93.97 -118.14    1.64   63.22  124.94  244.01
## beta1[1]     0.32    0.01   0.34   -0.36    0.09    0.32    0.55    0.99
## beta1[2]     0.39    0.01   0.38   -0.37    0.13    0.39    0.64    1.13
## beta1[3]     0.38    0.00   0.28   -0.18    0.19    0.38    0.56    0.92
## beta1[4]     0.60    0.01   0.42   -0.21    0.32    0.61    0.88    1.43
## beta1[5]     0.10    0.01   0.39   -0.66   -0.15    0.10    0.37    0.85
## beta1[6]     0.02    0.01   0.35   -0.64   -0.21    0.02    0.25    0.71
## beta1[7]     0.80    0.01   0.41   -0.03    0.52    0.79    1.08    1.59
## beta1[8]     0.03    0.00   0.24   -0.45   -0.13    0.03    0.19    0.52
## beta1[9]     1.22    0.01   0.54    0.14    0.85    1.22    1.56    2.31
## beta1[10]    0.12    0.01   0.33   -0.53   -0.11    0.12    0.34    0.76
## beta1[11]    0.22    0.01   0.45   -0.69   -0.08    0.21    0.51    1.11
## beta1[12]    0.32    0.00   0.26   -0.19    0.14    0.32    0.50    0.82
## sig        105.61    0.17  10.52   87.36   98.21  104.95  112.10  128.23
## lp__      -397.00    0.13   4.24 -406.44 -399.58 -396.69 -393.97 -389.73
##           n_eff Rhat
## beta0[1]   4000    1
## beta0[2]   4000    1
## beta0[3]   4000    1
## beta0[4]   3620    1
## beta0[5]   4000    1
## beta0[6]   4000    1
## beta0[7]   4000    1
## beta0[8]   4000    1
## beta0[9]   4000    1
## beta0[10]  4000    1
## beta0[11]  4000    1
## beta0[12]  4000    1
## beta1[1]   4000    1
## beta1[2]   4000    1
## beta1[3]   4000    1
## beta1[4]   4000    1
## beta1[5]   4000    1
## beta1[6]   4000    1
## beta1[7]   4000    1
## beta1[8]   4000    1
## beta1[9]   4000    1
## beta1[10]  4000    1
## beta1[11]  4000    1
## beta1[12]  4000    1
## sig        4000    1
## lp__       1113    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:11 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

各球団ごとの，切片と傾きの推定値が得られました。ところが，あれれ，2.5%と97.5%のパーセンタイル点をみると，プラスやマイナスに大きく散らばっています。これは各数値についての信頼できる幅が広い（広すぎる！）ことなので，何も言ったことになりません。もう少し，傾きや切片についてモデルを組み入れていきます。

平均とそこからのズレで考える

切片と傾きが，それぞれのチーム分計算されましたが，これを「係数の全体平均+各チームの差」のような形で表現することを考えます。すなわち，切片＝切片の平均+チームの差(雪片)，傾き＝傾きの平均+チームの差(傾き)，とします。ここでチームの差,と表現したのは，平均0で正負どちらの方向にも分布しうるものとして，正規分布を仮定します。コードで表現すると次のようになります。

次のコードをコピペしてHLM2.stanという名前で保存してください。

data{
  int<lower=0> N;                 // プレイヤーの数
  int<lower=0> K;                 // チームの数
  real X[N];                      // 独立変数
  real Y[N];                      // 従属変数
  int<lower=0,upper=K> teamID[N]; // プレイヤーがどこのチームに属しているかを表す識別子
}

parameters{
  real<lower=0> mean_beta0;     // 切片の平均。平均がマイナスってことはないだろうと信じています。
  real mean_beta1;             // 傾きの平均
  real beta0[K];          // 切片
  real beta1[K];          // 傾き
  real<lower=0> sig;      // 誤差成分
  real<lower=0> sig0;     // 切片の散らばり具合
  real<lower=0> sig1;     // 傾きの散らばり具合
}

transformed parameters{
  // 平均と差の関係をわかりやすくするためにここで先に計算しておく
  real b0[K]; //計算された切片
  real b1[K]; //計算された傾き
  for(k in 1:K){
    b0[k] = mean_beta0 + beta0[k];
    b1[k] = mean_beta1 + beta1[k];
  }
}

model{
  // データの成り立ち
  for(n in 1:N){
    Y[n] ~ normal(b0[teamID[n]] + b1[teamID[n]] * X[n] ,  sig);
  }
  // モデルの仮定
  for(k in 1:K){
    beta0[k] ~ normal(0,sig0); // 平均周りを散らばる切片
    beta1[k] ~ normal(0,sig1); // 平均周りを散らばる傾き
  }
}

これを次のコードで呼び出して分析します。

# stanモデルのコンパイル
model2 <- stan_model("HLM2.stan")
## recompiling to avoid crashing R session
# モデルフィット
fit2 <- sampling(model2,datastan)
## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
# 結果出力
fit2
## Inference for Stan model: HLM2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##               mean se_mean    sd    2.5%     25%     50%     75%   97.5%
## mean_beta0   50.27    0.67 30.06    4.12   26.62   47.12   69.63  116.19
## mean_beta1    0.28    0.00  0.08    0.11    0.23    0.29    0.34    0.43
## beta0[1]     -6.29    0.51 26.04  -68.46  -18.19   -3.30    6.54   44.80
## beta0[2]    -10.81    0.65 25.90  -72.98  -23.45   -6.20    3.62   34.13
## beta0[3]     15.76    0.93 28.10  -29.34   -1.85    9.49   30.17   84.52
## beta0[4]     -3.37    0.44 24.80  -59.71  -15.63   -1.92    8.97   46.52
## beta0[5]     -9.69    0.58 25.31  -70.47  -22.09   -5.21    4.30   34.82
## beta0[6]    -14.38    0.89 28.36  -85.71  -27.72   -8.03    2.50   31.22
## beta0[7]      2.51    0.47 24.45  -47.25   -9.70    1.33   14.65   55.52
## beta0[8]    -10.93    0.75 26.63  -74.66  -23.86   -6.24    3.90   36.45
## beta0[9]     12.00    0.83 28.73  -35.88   -3.72    6.58   24.42   83.57
## beta0[10]    16.04    0.90 28.64  -28.22   -1.24    9.34   29.41   88.35
## beta0[11]    -3.51    0.38 23.80  -58.03  -15.05   -1.68    8.81   44.96
## beta0[12]     6.05    0.40 23.46  -38.15   -6.48    3.22   17.60   60.70
## beta1[1]     -0.02    0.00  0.07   -0.18   -0.06   -0.02    0.02    0.12
## beta1[2]     -0.03    0.00  0.08   -0.21   -0.08   -0.02    0.01    0.12
## beta1[3]      0.07    0.00  0.08   -0.06    0.01    0.05    0.12    0.26
## beta1[4]      0.00    0.00  0.07   -0.14   -0.04    0.00    0.04    0.16
## beta1[5]     -0.04    0.00  0.08   -0.22   -0.09   -0.03    0.01    0.11
## beta1[6]     -0.07    0.00  0.08   -0.26   -0.12   -0.05   -0.01    0.06
## beta1[7]      0.03    0.00  0.07   -0.09   -0.01    0.02    0.07    0.20
## beta1[8]     -0.07    0.00  0.08   -0.25   -0.11   -0.05   -0.01    0.06
## beta1[9]      0.07    0.00  0.09   -0.07    0.01    0.05    0.13    0.30
## beta1[10]     0.05    0.00  0.08   -0.09    0.00    0.04    0.09    0.21
## beta1[11]    -0.01    0.00  0.08   -0.18   -0.06   -0.01    0.03    0.14
## beta1[12]     0.02    0.00  0.07   -0.12   -0.02    0.02    0.06    0.20
## sig         100.85    0.17  9.07   84.88   94.46  100.46  106.43  120.41
## sig0         26.65    1.01 18.18    3.11   12.76   22.81   37.32   68.20
## sig1          0.09    0.00  0.05    0.01    0.05    0.09    0.12    0.22
## b0[1]        43.98    0.79 38.77  -29.80   18.06   42.22   68.09  123.73
## b0[2]        39.47    0.89 38.38  -34.35   14.49   37.41   62.58  117.68
## b0[3]        66.03    1.17 40.91   -0.06   35.96   62.55   91.93  154.43
## b0[4]        46.90    0.72 37.71  -22.45   21.33   44.12   70.83  124.38
## b0[5]        40.58    0.82 36.37  -30.49   16.79   38.72   63.25  116.69
## b0[6]        35.89    0.99 39.22  -45.96   11.72   35.41   60.54  113.08
## b0[7]        52.78    0.80 38.39  -14.32   25.79   50.01   76.72  134.70
## b0[8]        39.34    0.99 38.00  -35.32   14.23   37.59   63.58  116.67
## b0[9]        62.28    1.12 40.98   -4.77   32.96   57.96   88.00  154.23
## b0[10]       66.31    1.15 41.37   -0.29   35.96   61.84   91.81  159.85
## b0[11]       46.77    0.68 36.67  -20.30   21.35   44.62   70.90  122.63
## b0[12]       56.32    0.75 36.63   -5.37   29.93   53.66   78.86  134.21
## b1[1]         0.26    0.00  0.10    0.05    0.20    0.27    0.33    0.44
## b1[2]         0.25    0.00  0.11    0.02    0.18    0.26    0.32    0.43
## b1[3]         0.35    0.00  0.11    0.14    0.28    0.35    0.42    0.55
## b1[4]         0.28    0.00  0.10    0.08    0.22    0.29    0.35    0.47
## b1[5]         0.24    0.00  0.11    0.00    0.18    0.25    0.32    0.44
## b1[6]         0.21    0.00  0.11   -0.03    0.15    0.22    0.29    0.40
## b1[7]         0.32    0.00  0.10    0.12    0.25    0.32    0.38    0.50
## b1[8]         0.22    0.00  0.10   -0.01    0.15    0.22    0.29    0.40
## b1[9]         0.36    0.00  0.12    0.13    0.28    0.35    0.43    0.60
## b1[10]        0.33    0.00  0.10    0.12    0.26    0.33    0.40    0.53
## b1[11]        0.27    0.00  0.11    0.05    0.20    0.27    0.34    0.46
## b1[12]        0.31    0.00  0.10    0.10    0.24    0.31    0.37    0.50
## lp__       -405.23    0.65  9.88 -423.41 -411.91 -405.75 -399.17 -384.62
##            n_eff Rhat
## mean_beta0  2027 1.00
## mean_beta1  2068 1.00
## beta0[1]    2657 1.00
## beta0[2]    1571 1.00
## beta0[3]     906 1.00
## beta0[4]    3149 1.00
## beta0[5]    1896 1.00
## beta0[6]    1021 1.00
## beta0[7]    2714 1.00
## beta0[8]    1247 1.00
## beta0[9]    1207 1.00
## beta0[10]   1012 1.00
## beta0[11]   4000 1.00
## beta0[12]   3514 1.00
## beta1[1]    2397 1.00
## beta1[2]    2809 1.00
## beta1[3]     885 1.01
## beta1[4]    3300 1.00
## beta1[5]    1790 1.00
## beta1[6]     881 1.00
## beta1[7]    1673 1.00
## beta1[8]     993 1.00
## beta1[9]     910 1.01
## beta1[10]   1287 1.00
## beta1[11]   2835 1.00
## beta1[12]   2671 1.00
## sig         2842 1.00
## sig0         327 1.00
## sig1         398 1.01
## b0[1]       2392 1.00
## b0[2]       1848 1.00
## b0[3]       1232 1.00
## b0[4]       2713 1.00
## b0[5]       1960 1.00
## b0[6]       1572 1.00
## b0[7]       2280 1.00
## b0[8]       1474 1.00
## b0[9]       1331 1.00
## b0[10]      1289 1.00
## b0[11]      2901 1.00
## b0[12]      2383 1.00
## b1[1]       2265 1.00
## b1[2]       1489 1.00
## b1[3]       1541 1.00
## b1[4]       2911 1.00
## b1[5]       1893 1.00
## b1[6]       1103 1.00
## b1[7]       2267 1.00
## b1[8]       1048 1.00
## b1[9]       1356 1.00
## b1[10]      1653 1.00
## b1[11]      2751 1.00
## b1[12]      2470 1.00
## lp__         233 1.00
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:48 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

この結果からわかることはなんでしょう？まず切片や傾きの情報に注目してみます。

# 切片と傾きの平均
print(fit2,pars=c("mean_beta0","mean_beta1"))
## Inference for Stan model: HLM2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##             mean se_mean    sd 2.5%   25%   50%   75%  97.5% n_eff Rhat
## mean_beta0 50.27    0.67 30.06 4.12 26.62 47.12 69.63 116.19  2027    1
## mean_beta1  0.28    0.00  0.08 0.11  0.23  0.29  0.34   0.43  2068    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:48 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

切片の平均が51(00万円),傾きの平均が0.28，ということがわかります。ヒット一本打つことで，年収が28万円上がるわけです。あくまでも平均で，個人差はありますでしょうが。

つぎにチームごとの特徴を見てみましょう。平均を置いたことで，相対比較になっています。

# 切片情報
print(fit2,pars="beta0")
## Inference for Stan model: HLM2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##             mean se_mean    sd   2.5%    25%   50%   75% 97.5% n_eff Rhat
## beta0[1]   -6.29    0.51 26.04 -68.46 -18.19 -3.30  6.54 44.80  2657    1
## beta0[2]  -10.81    0.65 25.90 -72.98 -23.45 -6.20  3.62 34.13  1571    1
## beta0[3]   15.76    0.93 28.10 -29.34  -1.85  9.49 30.17 84.52   906    1
## beta0[4]   -3.37    0.44 24.80 -59.71 -15.63 -1.92  8.97 46.52  3149    1
## beta0[5]   -9.69    0.58 25.31 -70.47 -22.09 -5.21  4.30 34.82  1896    1
## beta0[6]  -14.38    0.89 28.36 -85.71 -27.72 -8.03  2.50 31.22  1021    1
## beta0[7]    2.51    0.47 24.45 -47.25  -9.70  1.33 14.65 55.52  2714    1
## beta0[8]  -10.93    0.75 26.63 -74.66 -23.86 -6.24  3.90 36.45  1247    1
## beta0[9]   12.00    0.83 28.73 -35.88  -3.72  6.58 24.42 83.57  1207    1
## beta0[10]  16.04    0.90 28.64 -28.22  -1.24  9.34 29.41 88.35  1012    1
## beta0[11]  -3.51    0.38 23.80 -58.03 -15.05 -1.68  8.81 44.96  4000    1
## beta0[12]   6.05    0.40 23.46 -38.15  -6.48  3.22 17.60 60.70  3514    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:48 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
# 傾き情報
print(fit2,pars="beta1")
## Inference for Stan model: HLM2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## beta1[1]  -0.02       0 0.07 -0.18 -0.06 -0.02  0.02  0.12  2397 1.00
## beta1[2]  -0.03       0 0.08 -0.21 -0.08 -0.02  0.01  0.12  2809 1.00
## beta1[3]   0.07       0 0.08 -0.06  0.01  0.05  0.12  0.26   885 1.01
## beta1[4]   0.00       0 0.07 -0.14 -0.04  0.00  0.04  0.16  3300 1.00
## beta1[5]  -0.04       0 0.08 -0.22 -0.09 -0.03  0.01  0.11  1790 1.00
## beta1[6]  -0.07       0 0.08 -0.26 -0.12 -0.05 -0.01  0.06   881 1.00
## beta1[7]   0.03       0 0.07 -0.09 -0.01  0.02  0.07  0.20  1673 1.00
## beta1[8]  -0.07       0 0.08 -0.25 -0.11 -0.05 -0.01  0.06   993 1.00
## beta1[9]   0.07       0 0.09 -0.07  0.01  0.05  0.13  0.30   910 1.01
## beta1[10]  0.05       0 0.08 -0.09  0.00  0.04  0.09  0.21  1287 1.00
## beta1[11] -0.01       0 0.08 -0.18 -0.06 -0.01  0.03  0.14  2835 1.00
## beta1[12]  0.02       0 0.07 -0.12 -0.02  0.02  0.06  0.20  2671 1.00
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:48 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

切片は相対的に10番目のチームがよくて，6番目のチームが悪そうですね。10番目，6番目は・・・っと。どっちもパリーグか。

levels(baseball$team)
##  [1] "DeNA"         "オリックス"   "ソフトバンク" "ヤクルト"    
##  [5] "ロッテ"       "楽天"         "巨人"         "広島"        
##  [9] "阪神"         "西武"         "中日"         "日本ハム"

ヒットごとによる上昇率の良さは阪神で，悪いのは楽天と広島，ですね。 (チームの経営，運営方針というより，ここのデータが全選手ではなく成績上位の選手だから，という特徴があるかもしれません。念のため。)

こうしたチームごとの散らばりを表現しているのが次のパラメータ。

print(fit2,pars=c("sig","sig0","sig1"))
## Inference for Stan model: HLM2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean    sd  2.5%   25%    50%    75%  97.5% n_eff Rhat
## sig  100.85    0.17  9.07 84.88 94.46 100.46 106.43 120.41  2842 1.00
## sig0  26.65    1.01 18.18  3.11 12.76  22.81  37.32  68.20   327 1.00
## sig1   0.09    0.00  0.05  0.01  0.05   0.09   0.12   0.22   398 1.01
## 
## Samples were drawn using NUTS(diag_e) at Tue Dec 19 14:17:48 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

一つ目のsigがこのモデルで説明できない部分の散らばりです。sig0，sig1は係数の散らばりです。sigがsig0,sig1よりも大きな値になっています。すなわち，チームごとの差異よりも個別の違いのほうが大きいということで，やはり「実力主義」の世界なんだなあ，という感じです。

一般化線形モデルと階層線形モデル

前回の「一般化線形モデル」でお話ししたのは，正規分布してなさそうなデータ生成メカニズムに基づく分布モデルでも，パラメータをうまく置き換えてモデリングするという話でした。これは適用する確率分布に対する制限から解き放たれた「一般化」があったわけですが，階層線形モデルもこの考え方を応用できます。すなわち，階層化された一般化線形モデルというのも可能です。

例えば，チームごとの三振の数をHRの数で予測するモデルを考えると次のようになります。 まず三振の数はカウントデータなのでポアソン分布を使います。ポアソン分布のパラメータは正の数しかとりませんから，expで変換する必要があります。そこに注意して，後は平均と分散で係数の散らばりを表現します。

data{
  int<lower=0> N;                 // プレイヤーの数
  int<lower=0> K;                 // チームの数
  real X[N];                      // 独立変数
  int<lower=0> Y[N];              // 従属変数。整数型であることに注意
  int<lower=0,upper=K> teamID[N]; // プレイヤーがどこのチームに属しているかを表す識別子
}

parameters{
  real mean_beta0;        // 切片の平均
  real mean_beta1;        // 傾きの平均
  real beta0[K];          // 切片
  real beta1[K];          // 傾き
  real<lower=0> sig0;     // 切片の散らばり具合
  real<lower=0> sig1;     // 傾きの散らばり具合
}

transformed parameters{
  // 平均と差の関係をわかりやすくするためにここで先に計算しておく
  real b0[K]; //計算された切片
  real b1[K]; //計算された傾き
  for(k in 1:K){
    b0[k] = mean_beta0 + beta0[k];
    b1[k] = mean_beta1 + beta1[k];
  }
}

model{
  // データの成り立ち
  for(n in 1:N){
    // ポアソン分布を使っていることに注意。パラメータは一つで，それを直接モデリングしている。
    Y[n] ~ poisson(exp(b0[teamID[n]] + b1[teamID[n]] * X[n]));
  }
  // モデルの仮定
  for(k in 1:K){
    beta0[k] ~ normal(0,sig0); // 平均周りを散らばる切片
    beta1[k] ~ normal(0,sig1); // 平均周りを散らばる傾き
  }
}

これを呼び出すコードは次のようになります。実行内容・結果はともかく，何をやろうとしているかのイメージができるようになってください。

# データの作成
datastan <- list(N=nrow(base2),K=max(base2$team),X=baseball$HR,Y=baseball$K,teamID=base2$team)
# stanモデルのコンパイル
model3 <- stan_model("HLM3.stan")
# モデルフィット
fit3 <- sampling(model3,datastan)
# 結果出力
fit3

本日の課題

野球のデータを使って，「年収(pay)」を「ヒットの数(AB)」と「ホームランの数(HR)」の二変数で予測するモデルを作ります。

  1. 階層性を考えないモデル，すなわちチームごとの違いを考慮しないモデルを作って結果を示してください。
  2. 階層性を考え，各チームごとに違う傾き，切片を計算するモデルを作って結果を示してください。
  3. 階層性を考え，各チームごとの傾き・切片の平均値とそこからの相対的なズレで考えるモデルを作って結果を示してください。