理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

混合モデル(マルチレベルモデル)の基礎

投稿日:2021.5.2

忘れないように基本的な部分のみ貼っておきますが、名称だけでも統一していただきたいです・・・

  • 線形混合モデル ( linear\ mixed\ model)
  • 線形混合効果モデル ( linear\ mixed\ effect\ model)
  • 階層線形モデル ( hierarchical\  linear\ model )
  • マルチレベルモデル ( multilevel \ model )・・・etc

サンプルはパッケージlme4の「cbpp」を使用します
概要:牛肺疫(CBPP:contagious bovine pleuropneumonia)は、マイコプラズマを原因とするウシおよびスイギュウの感染症。本データセットは、エチオピアのBoji地区にある商業用の牛群(15グループ)における血清学的なCBPP発生状況を示したもの。これら15グループすべての牛から四半期ごとに血液サンプルを採取した追跡調査。 ​

library(lme4)
data(package="lme4") #パッケージのサンプルを見る

head(cbpp)
#    herd incidence size period
#  1    1         2   14      1
#  2    1         3   12      2
#  3    1         4    9      3
#  4    1         0    5      4
#  5    2         3   22      1
#  6    2         1   18      2

herd:牛のグループ名(15グループ)
incidence:各グループの各期間におけるCBPP発症数。
size:各グループに飼育されている牛飼養頭数
period:調査時期

サンプルの目的とは異なりますが「牛飼養頭数とCBPP発症頭数の関連性について」という設定で進めます。

全体の散布図と回帰直線

library(lattice)
xyplot(incidence  ~ size,  data = cbpp, col=1, type = c("p","r"))

f:id:yoshida931:20210502225944p:plain:w400

fit1 <- lm(incidence  ~ size,  data = cbpp)
summary(fit1)

#             Estimate Std. Error t value Pr(>|t|)   
# (Intercept) -0.31106    0.73554  -0.423  0.67405   
# size         0.13827    0.04389   3.150  0.00266 **

牛の数が1頭増えるごとに、発病する牛は0.13827頭増加する傾向にある(p=0.0027)。つまり牛の飼養頭数が10頭増すとCBPP感染牛が1.4頭増える傾向にあることが示唆されている。

調査時期を考慮した場合はどうなるか?

xyplot(incidence  ~ size, data = cbpp, group = period, pch = 19, type = c("p","r"), auto.key = list(pch = 19, corner=c(0,1), border = T, padding.text = 1.5, columns = 1), par.settings = simpleTheme(pch = 16))

f:id:yoshida931:20210502230145p:plain:w400

調査時期ごとの単回帰分析の結果

md <- lmList(incidence ~ size|period, data = cbpp)
summary(md)  
   
# Coefficients:
#    (Intercept) 
#     Estimate Std. Error    t value  Pr(>|t|)
# 1 -0.5048610   1.731796 -0.2915246 0.7719063
# 2 -0.1375479   1.323172 -0.1039531 0.9176397
# 3  1.3450519   1.390118  0.9675809 0.3381039
# 4 -0.1213868   1.066066 -0.1138642 0.9098204
#    size 
#      Estimate Std. Error    t value    Pr(>|t|)
# 1  0.24666516 0.08851373  2.7867445 0.007603578
# 2  0.08927203 0.07871405  1.1341308 0.262372119
# 3 -0.02452146 0.08995962 -0.2725829 0.786343078
# 4  0.05534211 0.07412426  0.7466127 0.458938830

両変数の関連に有意な傾向を示したのは、時期1のみである。時期3は有意ではないが、やや下降傾向にある。各調査の回帰直線を考慮した上で、モデルを構築する必要があります。

単回帰モデル

fit <- lm(incidence  ~ size,  data = cbpp)
summary(fit)
#            Estimate Std. Error t value Pr(>|t|)   
#(Intercept) -0.31106    0.73554  -0.423  0.67405   
#size         0.13827    0.04389   3.150  0.00266 **

期間をランダム効果として取り入れたモデル
ただし傾きと切片に共分散を仮定(共分散≠0)

fitm1 <-  lmer(incidence  ~ size + (size|period),  data = cbpp)
summary(fitm1)
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  period   (Intercept) 0.24363  0.4936        
#           size        0.01047  0.1023   -1.00
#  Residual             4.40888  2.0997        
# Number of obs: 56, groups:  period, 4
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept) -0.01094    0.68611  -0.016
# size         0.10386    0.06422   1.617

sizeは全期間を通して0.1023程度のバラツキがある. 各時期(period)の傾向を考慮したことで、傾きは0.13827から0.10386に減少しているが各グループの牛の数に応じて発症数は増加傾向にあることが理解できる.  時期(period)を考慮した場合、牛の数が1頭増えるごとに、発病する牛は0.10386頭増加する傾向にある.

それぞれの時期における推定値

ranef(fitm1)$period

#   (Intercept)        size
# 1  -0.6635655  0.13756593
# 2   0.1108923 -0.02298944
# 3   0.2900289 -0.06012683
# 4   0.2626443 -0.05444966

# 例)時期1の場合の発症数の推定値
# incidence  = -0.01094 + (-0.01094)*size + (-0.6635655) + 0.13756593*size
#              =(-0.01094 - 0.6635655)+(-0.01094 + 0.13756593)*size

他の混合モデルの例
切片を固定、傾きをランダムにしたモデル

(  )の中:0=切片固定、1=切片ランダム
fitm2 <-  lmer(incidence  ~ size + (0 + size|period),  data = cbpp)
summary(fitm2)

# Random effects:
#  Groups   Name Variance Std.Dev.
#  period   size 0.006005 0.0775  
#  Residual      4.430191 2.1048  
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.08343    0.64411   0.130
# size         0.10009    0.05508   1.817

傾きを固定、切片をランダムにしたモデル

fitm3 <-  lmer(incidence  ~ size + (1|period),  data = cbpp)
summary(fitm3)

# Random effects:
#  Groups   Name        Variance Std.Dev.
#  period   (Intercept) 1.506    1.227   
#  Residual             4.800    2.191   
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.24068    0.91873   0.262
# size         0.09878    0.04136   2.389

傾き、切片をランダムにしたモデル
最初に作ったfitm1と同じモデル

fitm4 <-  lmer(incidence  ~ size + (1 + size|period),  data = cbpp)
summary(fitm4)

#  Groups   Name        Variance Std.Dev. Corr 
#  period   (Intercept) 0.24363  0.4936        
#           size        0.01047  0.1023   -1.00
#  Residual             4.40888  2.0997        
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept) -0.01094    0.68611  -0.016
# size         0.10386    0.06422   1.617

null model(傾き=0、切片固定のモデル)

fitm5 <-  lmer(incidence  ~ 1 + (1|period),  data = cbpp)
summary(fitm5)

# Random effects:
#  Groups   Name        Variance Std.Dev.
#  period   (Intercept) 2.219    1.490   
#  Residual             5.124    2.264   
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)    1.714      0.804   2.131

null modelと比較

lm(incidence  ~ 1 , data = cbpp)

傾きと切片をそれぞれ独立したランダム変数としたモデル

lmer(incidence ~ size + (1|period) + (0 + size|period), data = cbpp)

グループを考慮した場合はどうなるか?

各グループ最大4回調査を実施していますので、各時期の頭数と発症数の関連性を見てみます。

xyplot(incidence  ~ size, data = cbpp, group = herd, pch = 19, type = c("p","r"), auto.key = list(pch = 19, corner=c(0,1), border = T, padding.text = 1.5, columns = 1), par.settings = simpleTheme(pch = 16))

f:id:yoshida931:20210505084307p:plain:w500
15グループあるので非常に複雑で、しかも切片・傾きともにランダム(バラバラ)であることには間違いないようです・・・といった具合にマルチレベルの分析は果てしなく続きます。