ランダム効果は傾きにも付与できる

　過分散なデータを対処する術として、切片にランダム効果を付与するGLMMを紹介した。切片以外にも、説明変数に対してランダム効果を仮定することも可能である。つまり、説明変数にかかる係数にランダムな変動がかかるということなので、線形回帰で考えればランダム効果に相当するグループごとに傾きが変わるようなデータになるということになる。この辺りは言葉で説明するよりもデータを見たほうが早いだろうから、さっそく解析をやってみよう。

具体的な解析

　では、以下のようなデータが得られているとき、どうやって解析していくかを確認してゆこう。group変数には、測定した各地点A,B,C……が格納されていると考えてほしい。

------------------------------------------------------

library(plotn)

library(glmmTMB)

library(MASS)

g <- 6

n <- 10

group <- rep(LETTERS[1:g], each = n)

x <- c(9.8, -1.3, 3.8, -9.4, -14, -13.8, 1.2, -8.5, 14.8, 

       0.1, 4.4, -11.8, 1.7, -12.7, -3.9, -6.5, 5, 3.7, 

       -8.9, -9, -9.7, -14.7, -9.6, 8.1, -10.5, 10.8, 11.3,

       -13.3, 7.8, -0.5, -14.6, 14.9, 7.5, 2.7, -1.6, -8.3,

       11.1, -4.7, -3.9, 0.3, 2.1, 0.7, -6.4, -8.5, -2.1, 

       7.2, -13.5, -9.7, 0.4, -7.1, -14.1, 6, 9.7, -6.1, 13,

       8.6, 7.1, 4.6, 3.6, -13)

y <- c(65.6, -0.2, 24.7, -51.6, -82.8, -82.2, 8.3, -50.3,

       96, 3.6, -2.1, 8.8, 2.1, 5.5, 2.5, 7.9, -0.1, -1.4,

       1, 1.8, 9.5, 17.7, 15.2, -7.4, 10.9, -6.7, -7.3, 

       23.2, -6, 1.3, 8.6, -22.3, -9.9, -8.7, -1.4, 2.5, 

       -19.3, 0.1, 5.4, -0.6, -5, -5, -1.9, -3.5, -2.6, 

       -4.7, 0.1, -4.4, -4.5, -1.8, -1.6, -0.8, -0.3, 0.6,

       0.9, -0.5, 1, 4.2, 1.6, -2)

d <- data.frame(x = x, group = group, y = y)

plotn(y ~ x, data = d)#図1の描画

------------------------------------------------------

なにやら、一部のデータで外れ値が出ている感じがする。そこでgroupごとに色分けして、図示してみる。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")#図2の描画

------------------------------------------------------

明らかにグループごとに傾きが異なっている。こういう時には、GLMMによっての解析が活きるときである。glmmTMBパッケージ内の、glmmTMB関数を使って解析してみよう。下記のようにすると、切片と説明変数xに対してランダム効果を仮定することができる。

------------------------------------------------------

fit1 <- glmmTMB(y ~ x + (x|group), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル

------------------------------------------------------

もし、切片はランダム効果がなくxにだけランダム効果を仮定するのであれば、(0 + x|group)とすることで目的に沿った解析が可能である。また、(x|group)とすれば切片とxの間のランダム効果に相関がある可能性を考慮した解析になるが（多変量正規分布からランダム効果が生成したと仮定する）、切片とxの間のランダム効果が独立に決まると仮定すると(1|group) + (0 + x|group)という表記になる。この辺りはちょっと込み入ったことを考えていることになるので、今は忘れてもらってもよい。

　では、解析結果を眺めてみる。

------------------------------------------------------

fit2 <- summary(fit1)

fit2

##  Family: gaussian  ( identity )

## Formula:          y ~ x + (x | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    345.6    358.1   -166.8    333.6       54 

## 

## Random effects:

## 

## Conditional model:

##  Groups   Name        Variance Std.Dev. Corr 

##  group    (Intercept) 8.952    2.992         

##           x           6.443    2.538    0.56 

##  Residual             6.228    2.496         

## Number of obs: 60, groups:  group, 6

## 

## Dispersion estimate for gaussian family (sigma^2): 6.23 

## 

## Conditional model:

##             Estimate Std. Error z value Pr(>|z|)

## (Intercept) -0.04344    1.26888  -0.034    0.973

## x            0.58966    1.03706   0.569    0.570

------------------------------------------------------

以前紹介したようにRandom effectsの項がランダム効果であり、ランダム効果をもたらす正規分布の標準偏差が推定されている。今回だと、groupランダム変数の標準偏差は、切片に対して2.992、xに対して2.538、そして切片とxの相関は0.56程度の多変量正規分布から生成されている、と推定したようだ。

　Conditional modelが固定効果の推定値である。切片は-0.043、xの係数は0.590程度のようだ。このデータはxとはあまり関係のないデータのようだ。

　ちなみに普通の回帰分析をしてみると下記の通り。

------------------------------------------------------

summary(lm(y ~ x))

## 

## Call:

## lm(formula = y ~ x)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -72.908  -6.571   1.037   8.257  85.732 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)  

## (Intercept) -0.09181    3.16000  -0.029   0.9769  

## x            0.70002    0.36035   1.943   0.0569 .

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 24.12 on 58 degrees of freedom

## Multiple R-squared:  0.06109,    Adjusted R-squared:  0.0449 

## F-statistic: 3.774 on 1 and 58 DF,  p-value: 0.05692

------------------------------------------------------

普通の回帰分析ではxに関して有意ではないが、5%近くになっており、GLMMとは傾向が異なっていることがわかる。

　では、回帰曲線を描いてみよう。固定効果だけを使った回帰曲線は下記の通り。

------------------------------------------------------

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(y ~ x, data = d)#図3の描画

xx <- seq(min(d$x), max(d$x), length = 200)

yy <- b0e + b1e * xx

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

下記のようにglmmTMBパッケージなどでは、ranef関数が実装されており、ランダム効果の推定値を取り出すことができる（固定効果を取り出したいときはfixef関数もある）。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)          x

## A   3.8259934  5.5878326

## B   0.9388431 -0.9957598

## C   2.5544760 -1.6537775

## D  -3.8665347 -1.6802266

## E  -3.5885883 -0.7562610

## F   0.1358130 -0.5018164

------------------------------------------------------

これを使うことで、各グループごとの回帰曲線も描くことができる。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")#図4の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- (b0e + ranef(fit1)$cond$group[i,1]) + (b1e + ranef(fit1)$cond$group[i,2]) * xx

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

ランダム効果によって切片やxの傾きがばらつくことで、様々な傾向が検出されていることがわかるだろう。

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 3#固定効果：正規分布の標準偏差

is <- 3#切片ランダム効果：正規分布の標準偏差

ss <- 2#傾きランダム効果：正規分布の標準偏差

g <- 6

n <- 10

x <- round(runif(n*g, -15, 15), digits = 1)

group <- rep(LETTERS[1:g], each = n)

i_sds <- rep(rnorm(g, mean = 0, sd = is), each = n)

s_sds <- rep(rnorm(g, mean = 0, sd = ss), each = n)

y <- rnorm(n = n*g, mean = (b0+i_sds) + (b1+s_sds)*x, sd = s)

------------------------------------------------------

つまり、傾きは0であることが正解である。また、ランダム効果の正規分布の標準偏差は切片は3、傾きは2が正解だったわけであるが、今回はそれぞれ2.992、2.538と妥当な値が推定されている。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。また、通常の線形回帰も行って、係数などの比較も行おう。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 3

is <- 3

ss <- 2

g <- 6

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

p1_l <- NULL

p2_l <- NULL

b0e_l <- NULL

b1e_l <- NULL

is_e <- NULL

ss_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(n*g, -15, 15), digits = 1)

  i_sds <- rep(rnorm(g, mean = 0, sd = is), each = n)

  s_sds <- rep(rnorm(g, mean = 0, sd = ss), each = n)

  y <- rnorm(n = n*g, mean = (b0+i_sds) + (b1+s_sds)*x, sd = s)

  d <- data.frame(x = x, group = group, y = y)

  fit1 <- glmmTMB(y ~ x + (x|group), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル

  fit2 <- summary(fit1)

  fit3 <- lm(y ~ x, data = d)

  fit4 <- summary(fit3)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  p1_l <- c(p1_l, coef(fit4)[1,4])

  p2_l <- c(p2_l, coef(fit4)[2,4])

  b0e_l <- c(b0e_l, coef(fit4)[1,1])

  b1e_l <- c(b1e_l, coef(fit4)[2,1])

  is_e <- c(is_e, as.numeric(fit2$varcor$cond$group[1,1]))

  ss_e <- c(ss_e, as.numeric(fit2$varcor$cond$group[2,2]))

}

------------------------------------------------------

シミュレーションできれいなデータを生成してもエラーがかなり出るのだが、やはり混合モデルの収束の難しさが出ていると考えてよいだろう。エラーが出た分のデータを除いて解析した結果は下記の通り。

------------------------------------------------------

histn(p1[!is.na(p1)], xlab = "P value", ylab = "頻度")#図5の描画

histn(p2[!is.na(p2)], xlab = "P value", ylab = "頻度")#図6の描画

sum(p1 < 0.05, na.rm = T)

## [1] 452

sum(p2 < 0.05, na.rm = T)

## [1] 132

sum(!is.na(p1))

## [1] 998

sum(!is.na(p2))

## [1] 998

------------------------------------------------------

切片の検出力は45.3%、xの係数の危険率は13.2%程度となった。実はGLMMでもちょっと危険率は高めである。一方で、通常の線形回帰を確認してみよう。

------------------------------------------------------

histn(p1_l, xlab = "P value", ylab = "頻度")#図7の描画

histn(p2_l, xlab = "P value", ylab = "頻度")#図8の描画

sum(p1_l < 0.05, na.rm = T)

## [1] 217

sum(p2_l < 0.05, na.rm = T)

## [1] 587

------------------------------------------------------

切片の検出力は21.7%、xの係数の危険率は58.7%程度となった。非常に高い危険率となり、GLMMのほうが信用度が高いだろう。併せて切片の検出力も落ちている。

　さらに推定値のヒストグラムや平均、分散を確認してみよう。

------------------------------------------------------

histn(b0e[!is.na(p1)], xlab = "b0", ylab = "頻度")#図9の描画

histn(b1e[!is.na(p2)], xlab = "b1", ylab = "頻度")#図10の描画

mean(b0e[!is.na(p1)])

## [1] 1.972333

mean(b1e[!is.na(p2)])

## [1] 0.02653856

var(b0e[!is.na(p1)])

## [1] 1.717112

var(b1e[!is.na(p2)])

## [1] 0.647029

------------------------------------------------------

　一方、線形回帰の推定値は下記の通り。

------------------------------------------------------

histn(b0e_l, xlab = "b0", ylab = "頻度")#図11の描画

histn(b1e_l, xlab = "b1", ylab = "頻度")#図12の描画

mean(b0e_l)

## [1] 1.860626

mean(b1e_l)

## [1] 0.03875636

var(b0e_l)

## [1] 6.302271

var(b1e_l)

## [1] 0.7061915

------------------------------------------------------

平均的な推定値はGLMM、線形回帰の両方で似たり寄ったりである。注目するべきは、切片の推定値のばらつきだろうか。GLMMのほうが精度良く推定できることがわかる。

　最後、切片とxの係数のランダム効果の分散を確認しておこう。

------------------------------------------------------

histn(is_e[!is.na(p1)], ylab = "頻度")#図13の描画

histn(ss_e[!is.na(p2)], ylab = "頻度")#図14の描画

mean(is_e[!is.na(p1)])

## [1] 7.494322

mean(ss_e[!is.na(p2)])

## [1] 3.346453

------------------------------------------------------

それぞれの分散の平均は切片では7.49、xの係数では3.35であり、元の設定の3^2 = 9、2^2 = 4の5/6となっていることがわかる。

　次は以下のような例を考えてみよう。

------------------------------------------------------

g <- 6

n <- 10

group <- rep(LETTERS[1:g], each = n)

x <- c(-8.2, 11.4, 7.2, -7.1, -5.5, -9.1, 1.3, 7.7, -5.1,

       0.4, 7.3, 11.7, 0, -0.7, 5.6, -11.5, 10.4, -10.9,

       5.4, -13.4, -10.6, 14.5, 14.7, 2.5, 13.8, 1, 7.9, 

       14.8, 4, 8.7, -5.2, -7.9, 12.6, 5.2, -12.3, -2.9, 

       8.9, -10.9, 1.3, -7.8, -11.6, -12.9, -2.3, -14.8, 

       -1.2, 6.4, -4.8, 11, 8.5, 2.8, 8.7, -8.9, -12.9, 

       -12.9, -11.5, 13, -6.1, 14.8, 13.5, 9.4)

y <- c(-1.4, 4.1, -0.5, 3.8, 1.6, 1.7, 3.3, 0.5, -1, -0.7,

       -23.4, -26.3, 0.2, -3.7, -16.3, 22.5, -28.7, 25.1, 

       -12.4, 33.7, 53.7, -57.5, -58.9, -3.5, -53.4, 3.5,

       -31.5, -59.2, -13, -32.3, 16.2, 15.5, -18.6, -8.5, 

       19, 9.8, -11.8, 18.3, 4.9, 12.5, 37.2, 48.5, 5.6, 

       53.9, 2.2, -27.4, 10.5, -46.1, -41.2, -15.2, -50, 

       44.3, 60, 56, 48, -65.3, 18.3, -74.9, -74.9, -51.7)

d <- data.frame(x = x, group = group, y = y)

plotn(y ~ x, data = d)#図15の描画

------------------------------------------------------

groupごとに色分けして、図示してみる。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")#図16の描画

------------------------------------------------------

このデータも明らかにグループごとに傾きが異なっている。glmmTMB関数を使って解析してみよう。

------------------------------------------------------

fit1 <- glmmTMB(y ~ x + (x|group), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: gaussian  ( identity )

## Formula:          y ~ x + (x | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    360.9    373.4   -174.4    348.9       54 

## 

## Random effects:

## 

## Conditional model:

##  Groups   Name        Variance Std.Dev. Corr 

##  group    (Intercept) 15.971   3.996         

##           x            2.945   1.716    0.31 

##  Residual              8.421   2.902         

## Number of obs: 60, groups:  group, 6

## 

## Dispersion estimate for gaussian family (sigma^2): 8.42 

## 

## Conditional model:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  -0.3596     1.6810  -0.214    0.831    

## x            -2.8494     0.7020  -4.059 4.93e-05 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------

今回だと、groupランダム変数の標準偏差は、切片に対して3.996、xに対して1.716、そして切片とxの相関は0.31程度の多変量正規分布から生成されている、と推定したようだ。

　Conditional modelが固定効果の推定値である。切片は-0.036、xの係数は-2.85程度のようだ。このデータはxとはあまり関係のないデータのようだ。

　ちなみに普通の回帰分析をしてみると下記の通り。

------------------------------------------------------

summary(lm(y ~ x))

## 

## Call:

## lm(formula = y ~ x)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -28.596  -9.084   0.792   9.274  43.600 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  -2.5620     1.9409   -1.32    0.192    

## x            -3.2401     0.2084  -15.54   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 15 on 58 degrees of freedom

## Multiple R-squared:  0.8064, Adjusted R-squared:  0.8031 

## F-statistic: 241.6 on 1 and 58 DF,  p-value: < 2.2e-16

------------------------------------------------------

　では、回帰曲線を描いてみよう。固定効果だけを使った回帰曲線は下記の通り。

------------------------------------------------------

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(y ~ x, data = d)#図17の描画

xx <- seq(min(d$x), max(d$x), length = 200)

yy <- b0e + b1e * xx

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

各グループごとの回帰曲線も描くと下記の通り。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")#図18の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- (b0e + ranef(fit1)$cond$group[i,1]) + (b1e + ranef(fit1)$cond$group[i,2]) * xx

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 0

b1 <- -3

s <- 3#固定効果：正規分布の標準偏差

is <- 3#切片ランダム効果：正規分布の標準偏差

ss <- 2#傾きランダム効果：正規分布の標準偏差

g <- 6

n <- 10

x <- round(runif(n*g, -15, 15), digits = 1)

group <- rep(LETTERS[1:g], each = n)

i_sds <- rep(rnorm(g, mean = 0, sd = is), each = n)

s_sds <- rep(rnorm(g, mean = 0, sd = ss), each = n)

y <- rnorm(n = n*g, mean = (b0+i_sds) + (b1+s_sds)*x, sd = s)

------------------------------------------------------

つまり、切片は0、傾きは-3であることが正解である。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。また、通常の線形回帰も行って、係数などの比較も行おう。

------------------------------------------------------

b0 <- 0

b1 <- -3

s <- 3

is <- 3

ss <- 2

g <- 6

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

p1_l <- NULL

p2_l <- NULL

b0e_l <- NULL

b1e_l <- NULL

is_e <- NULL

ss_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(n*g, -15, 15), digits = 1)

  i_sds <- rep(rnorm(g, mean = 0, sd = is), each = n)

  s_sds <- rep(rnorm(g, mean = 0, sd = ss), each = n)

  y <- rnorm(n = n*g, mean = (b0+i_sds) + (b1+s_sds)*x, sd = s)

  d <- data.frame(x = x, group = group, y = y)

  fit1 <- glmmTMB(y ~ x + (x|group), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル

  fit2 <- summary(fit1)

  fit3 <- lm(y ~ x, data = d)

  fit4 <- summary(fit3)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  p1_l <- c(p1_l, coef(fit4)[1,4])

  p2_l <- c(p2_l, coef(fit4)[2,4])

  b0e_l <- c(b0e_l, coef(fit4)[1,1])

  b1e_l <- c(b1e_l, coef(fit4)[2,1])

  is_e <- c(is_e, as.numeric(fit2$varcor$cond$group[1,1]))

  ss_e <- c(ss_e, as.numeric(fit2$varcor$cond$group[2,2]))

}

------------------------------------------------------

エラーが出た分のデータを除いて解析した結果は下記の通り。

------------------------------------------------------

histn(p1[!is.na(p1)], xlab = "P value", ylab = "頻度")#図19の描画

histn(p2[!is.na(p2)], xlab = "P value", ylab = "頻度")#図20の描画

sum(p1 < 0.05, na.rm = T)

## [1] 129

sum(p2 < 0.05, na.rm = T)

## [1] 946

sum(!is.na(p1))

## [1] 997

sum(!is.na(p2))

## [1] 997

------------------------------------------------------

切片の危険率は12.9%、xの係数の検出力は94.9%程度となった。やはりGLMMでもちょっと危険率は高めである。一方で、通常の線形回帰を確認してみよう。

------------------------------------------------------

histn(p1_l, xlab = "P value", ylab = "頻度")#図21の描画

histn(p2_l, xlab = "P value", ylab = "頻度")#図22の描画

sum(p1_l < 0.05, na.rm = T)

## [1] 109

sum(p2_l < 0.05, na.rm = T)

## [1] 998

------------------------------------------------------

切片の危険率は10.9%、xの係数の検出力は99.8%程度となった。今回の例ではGLMMとは大差がないことがわかる。

　さらに推定値のヒストグラムや平均、分散を確認してみよう。

------------------------------------------------------

histn(b0e[!is.na(p1)], xlab = "b0", ylab = "頻度")#図23の描画

histn(b1e[!is.na(p2)], xlab = "b1", ylab = "頻度")#図24の描画

mean(b0e[!is.na(p1)])

## [1] 0.001394154

mean(b1e[!is.na(p2)])

## [1] -2.988343

var(b0e[!is.na(p1)])

## [1] 1.632033

var(b1e[!is.na(p2)])

## [1] 0.6853282

------------------------------------------------------

　一方、線形回帰の推定値は下記の通り。

------------------------------------------------------

histn(b0e_l, xlab = "b0", ylab = "頻度")#図25の描画

histn(b1e_l, xlab = "b1", ylab = "頻度")#図26の描画

mean(b0e_l)

## [1] -0.04458816

mean(b1e_l)

## [1] -2.996494

var(b0e_l)

## [1] 5.860064

var(b1e_l)

## [1] 0.724934

------------------------------------------------------

この場合でも、平均的な推定値はGLMM、線形回帰の両方で似たり寄ったりである。そして、切片に関してはGLMMのほうがばらつきが小さく推定できることがわかる。

　最後、切片とxの係数のランダム効果の分散を確認しておこう。

------------------------------------------------------

histn(is_e[!is.na(p1)], ylab = "頻度")#図27の描画

histn(ss_e[!is.na(p2)], ylab = "頻度")#図28の描画

mean(is_e[!is.na(p1)])

## [1] 7.131953

mean(ss_e[!is.na(p2)])

## [1] 3.381329

------------------------------------------------------

それぞれの分散の平均は切片では7.13、xの係数では3.38であり、元の設定の3^2 = 9、2^2 = 4の5/6となっていることがわかる。

　最後、シミュレーションだけだが、切片とxの係数のランダム効果間で相関がある場合を示しておく。結論を先に述べると、この状況でも先述の推定精度と大差はない。今回は、MASSパッケージのmvrnorm関数を使うことで、相関のあるランダム効果を作ることにする。今回はランダム効果の分散だけでなく共分散も記録しておこう。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 3

is <- 3

ss <- 2

r <- -0.9#相関係数

g <- 6

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

p1_l <- NULL

p2_l <- NULL

b0e_l <- NULL

b1e_l <- NULL

is_e <- NULL

ss_e <- NULL

cv_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(n*g, -15, 15), digits = 1)

  sds <- mvrnorm(g, mu = rep(0, 2), 

                 Sigma = matrix(c(is^2, r*is*ss, r*is*ss, ss^2), ncol = 2))

  i_sds <- rep(sds[,1], each = n)

  s_sds <- rep(sds[,2], each = n)

  y <- rnorm(n = n*g, mean = (b0+i_sds) + (b1+s_sds)*x, sd = s)

  d <- data.frame(x = x, group = group, y = y)

  fit1 <- glmmTMB(y ~ x + (x|group), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル

  fit2 <- summary(fit1)

  fit3 <- lm(y ~ x, data = d)

  fit4 <- summary(fit3)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  p1_l <- c(p1_l, coef(fit4)[1,4])

  p2_l <- c(p2_l, coef(fit4)[2,4])

  b0e_l <- c(b0e_l, coef(fit4)[1,1])

  b1e_l <- c(b1e_l, coef(fit4)[2,1])

  is_e <- c(is_e, as.numeric(fit2$varcor$cond$group[1,1]))

  ss_e <- c(ss_e, as.numeric(fit2$varcor$cond$group[2,2]))

  cv_e <- c(cv_e, as.numeric(fit2$varcor$cond$group[1,2]))#共分散

}

plotn(y ~ x, data = d, mode = "m", group = "group")#図29の図示

------------------------------------------------------

エラーが出た分のデータを除いて解析した結果は下記の通り。

------------------------------------------------------

histn(p1[!is.na(p1)], xlab = "P value", ylab = "頻度")#図30の描画

histn(p2[!is.na(p2)], xlab = "P value", ylab = "頻度")#図31の描画

sum(p1 < 0.05, na.rm = T)

## [1] 441

sum(p2 < 0.05, na.rm = T)

## [1] 135

sum(!is.na(p1))

## [1] 981

sum(!is.na(p2))

## [1] 981

------------------------------------------------------

切片の検出力は44.1%、xの係数の危険率は13.8%程度となった。一方で、通常の線形回帰を確認してみよう。

------------------------------------------------------

histn(p1_l, xlab = "P value", ylab = "頻度")#図32の描画

histn(p2_l, xlab = "P value", ylab = "頻度")#図33の描画

sum(p1_l < 0.05, na.rm = T)

## [1] 225

sum(p2_l < 0.05, na.rm = T)

## [1] 584

------------------------------------------------------

切片の検出力は22.5%、xの係数の危険率は58.4%程度となった。やはり、線形回帰は傾き推定に弱いことがわかる。

　さらに推定値のヒストグラムや平均、分散を確認してみよう。

------------------------------------------------------

histn(b0e[!is.na(p1)], xlab = "b0", ylab = "頻度")#図34の描画

histn(b1e[!is.na(p2)], xlab = "b1", ylab = "頻度")#図35の描画

mean(b0e[!is.na(p1)])

## [1] 2.036088

mean(b1e[!is.na(p2)])

## [1] -0.03693563

var(b0e[!is.na(p1)])

## [1] 1.822456

var(b1e[!is.na(p2)])

## [1] 0.6886941

------------------------------------------------------

　一方、線形回帰の推定値は下記の通り。

------------------------------------------------------

histn(b0e_l, xlab = "b0", ylab = "頻度")#図36の描画

histn(b1e_l, xlab = "b1", ylab = "頻度")#図37の描画

mean(b0e_l)

## [1] 2.088242

mean(b1e_l)

## [1] -0.03216513

var(b0e_l)

## [1] 5.885132

var(b1e_l)

## [1] 0.7260289

------------------------------------------------------

この場合でも、平均的な推定値はGLMM、線形回帰の両方で似たり寄ったりである。そして、切片に関してはGLMMのほうがばらつきが小さく推定できることがわかる。

　最後、切片とxの係数のランダム効果の分散を確認しておこう。

------------------------------------------------------

histn(is_e[!is.na(p1)], ylab = "頻度")#図38の描画

histn(ss_e[!is.na(p2)], ylab = "頻度")#図39の描画

histn(cv_e[!is.na(p1)], ylab = "頻度")#図40の描画

mean(is_e[!is.na(p1)])

## [1] 7.799051

mean(ss_e[!is.na(p2)])

## [1] 3.330479

mean(cv_e[!is.na(p1)])

## [1] -4.582304

------------------------------------------------------

それぞれの分散の平均は切片では7.80、xの係数では3.33であり、元の設定の3^2 = 9、2^2 = 4の5/6となっていることがわかる。また、共分散は-4.58であり、3×2×(-0.9) = 5.4の5/6となっている。

おわりに

　さて、以上のように、私もシミュレーションしてみてわかったが、傾きに関しては平均的な推定精度はGLMMも線形回帰も似たり寄ったりであることがわかってきた。切片に関してはGLMMのほうが格段に推定精度が良い。この点から、モデルの収束などを考えると、切片項にだけランダム効果を仮定するというのはかなりリーズナブルな選択になりえることがわかる。必要に応じて、通常の線形回帰も考慮する必要もあるだろう。ただし、これはあくまで全体傾向と郡内傾向が同じ場合である。次回紹介するように、全体傾向と郡内傾向が異なる場合はGLMMを実行しないと解析に失敗することになる。