毎回やる儀式

バッターのデータ，batter.csvをb.dataというオブジェクトに入れているとします。 できてないひとは，次の通りやってみよう。

b.data <- read.csv("batter.csv")

今日使うパッケージ

階層線形モデルって何？

線形モデルは今までいろいろやってきて，しかも前回それを「一般化」して様々な分布に対応させました。今回はさらに，モデルの中に階層構造を入れることができます。

階層構造というのは，ネストされたデータとも言います。イメージはマトリョーシカ。 「ある都道府県の中の，ある地区の，ある学校のなかの一人一人のデータ」というのは，都道府県，地区，学校，という枠の中でネストされている，と表現します。

たとえば学習時間が成績を予測する，という簡単な回帰分析を考えた時に，地区や学校の影響なしで学習時間と成績の関係をみたいでしょう。しかし地区や学校によって平均が違ったり，回帰係数が違ったりする可能性もあるので，データをそのまま使うのはちょっと・・・というときにこのモデルを使います。

ちなみに別名として 混合モデル(Mixed Model),マルチレベルモデル と言われることもあります。

階層線形モデルへのルート

階層線形モデルは，回帰分析の「階層」への拡張，という理解の仕方もできます。

また，「反復測定分散分析の拡張」という理解の仕方もできます。反復測定ANOVA(群内要因ともいう)のデータは，個人を単位にネストされていると考えることができるからです。この時に，個人差を抜きにして実験の効果を見たいというのは，まさに上の階層データの例の別側面だからです。このほかに，実験刺激群の違いを考慮してモデリングしたい，という時にも使えます。

このように，階層線形モデルは，社会的なデータから実験室実験のデータにまで，広く応用することができるものです。

階層線形モデルのその後

ちなみに，様々な入り口から階層線形モデルに入ってくることができますが，その次には「そういった違いを生んだ要因はどういう特徴をしているか」という展開を迎えることになります（階層ベイズ）。こうしたアプローチについては，次回以降触れることになります。

実際にやってみましょう

実際にやってみましょう。まずは回帰。普通の重回帰。 年棒が1000万円あがったらホームランが何本増えるか，見てみましょう。

HR.lm <- lm(HR~pay,data=b.data)
summary(HR.lm)
## 
## Call:
## lm(formula = HR ~ pay, data = b.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.262  -6.813  -2.671   3.508  29.305 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.6419     1.8671   3.557 0.000641 ***
## pay           0.2567     0.1033   2.484 0.015142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.466 on 78 degrees of freedom
## Multiple R-squared:  0.0733, Adjusted R-squared:  0.06142 
## F-statistic: 6.169 on 1 and 78 DF,  p-value: 0.01514

図でも確認しておきます。

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.2
g <- ggplot(b.data,aes(x=pay,y=HR)) + geom_point()
g <- g + geom_smooth(method="lm")
g

そしてこれが，チームごとに同じ傾向があるかを見てみたいと思います。 先に図から見ましょう。

g <- ggplot(b.data,aes(x=pay,y=HR,color=team)) + geom_point()
g

さて，これを階層線形モデルで確認します。階層線形モデルを実行するパッケージは今の所lmer4が使い勝手も良いようです。このパッケージはlmerTestパッケージから使うともっと便利になるので，そうします。

library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.2.2
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.2.2
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 3.2.2
## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step
HR.lmer <- lmer(HR~pay+(1|team),data=b.data)
summary(HR.lmer)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: HR ~ pay + (1 | team)
##    Data: b.data
## 
## REML criterion at convergence: 585.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3546 -0.6617 -0.2879  0.3941  3.0330 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team     (Intercept)  3.667   1.915   
##  Residual             86.265   9.288   
## Number of obs: 80, groups:  team, 12
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   6.5642     1.9400 39.7000   3.384  0.00162 **
## pay           0.2661     0.1034 77.9800   2.573  0.01198 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## pay -0.792

ここでRandom effectsとあるところがチームごとのばらつきによる部分です。切片は標準偏差1.915分ぐらいのばらつきがあるようです。

チームの効果だけを確認したいときは次のように書きます。

HR.lmer.null <- lmer(HR~(1|team),data=b.data)
summary(HR.lmer.null)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: HR ~ (1 | team)
##    Data: b.data
## 
## REML criterion at convergence: 588.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1235 -0.6881 -0.2784  0.4173  2.8067 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team     (Intercept)  0.9877  0.9938  
##  Residual             94.5552  9.7239  
## Number of obs: 80, groups:  team, 12
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   10.480      1.127 10.409   9.301 2.29e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ここで，チームによる分散が0.987,残差(residual)の分散が94.555とあります。チームによる分散を考慮することの重要性は，この比率を次のようにして計算します。

0.987/(0.987+94.555)
## [1] 0.01033054

今回，この値が0.001と小さいものであったため，実はチームごとの分散はそれほど考慮する必要がない，ということがわかります。この指標は級内相関 と呼ばれ，0.1を超えたら階層モデルをする意義がある，とされています(今回はそれほど適切な例ではないことになりますが，続けます。)

ちなみにggplot2をつかって，チームごとの回帰モデルを図示しておきましょう。

#プロット（上の続き）
g <- g + geom_smooth(method="lm",se=FALSE)
g

プロットみると，チームごとに切片も違っていますけど，傾きが違うものがあるようです。これも検討してみましょう。

library(lmerTest)
HR.lmer2 <- lmer(HR~pay+(1+pay|team),data=b.data)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0289856 (tol =
## 0.002, component 1)
summary(HR.lmer)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: HR ~ pay + (1 | team)
##    Data: b.data
## 
## REML criterion at convergence: 585.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3546 -0.6617 -0.2879  0.3941  3.0330 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team     (Intercept)  3.667   1.915   
##  Residual             86.265   9.288   
## Number of obs: 80, groups:  team, 12
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   6.5642     1.9400 39.7000   3.384  0.00162 **
## pay           0.2661     0.1034 77.9800   2.573  0.01198 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## pay -0.792

そうそう，一般化線形モデリングで説明したように，これはポアソン分布を仮定したほうがいいのでした。

HR.glmer <- glmer(HR~pay+(1+pay|team),data=b.data,family="poisson")
summary(HR.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: HR ~ pay + (1 + pay | team)
##    Data: b.data
## 
##      AIC      BIC   logLik deviance df.resid 
##    853.5    865.4   -421.7    843.5       75 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8453 -2.1163 -0.5706  1.4204  7.3184 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. Corr
##  team   (Intercept) 0.0618985 0.2488       
##         pay         0.0001665 0.0129   0.57
## Number of obs: 80, groups:  team, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.906048   0.103071  18.493  < 2e-16 ***
## pay         0.026697   0.005699   4.685  2.8e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## pay -0.195

ま，結果は今ひとつでしたが，こういうアプローチができるということです。

# モデル同士の比較をするも，最初ん回帰が一番いいみたいw
anova(HR.lmer,HR.lmer2,HR.glmer)
## refitting model(s) with ML (instead of REML)
## Data: b.data
## Models:
## object: HR ~ pay + (1 | team)
## ..2: HR ~ pay + (1 + pay | team)
## ..1: HR ~ pay + (1 + pay | team)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## object  4 592.53 602.05 -292.26   584.53                             
## ..2     5 853.48 865.39 -421.74   843.48   0.00      1          1    
## ..1     6 595.36 609.65 -291.68   583.36 260.12      1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

本日の課題