群内傾向を考慮しないとどうなるか？

　これまでGLMMを使うことで、主には興味がない効果＝ランダム効果を考慮しつつも、主に興味のある効果＝固定効果を適切に推定することが可能となってきた。今回は、ある条件の元、群内傾向を考慮し忘れたことによって起こる、解釈の間違いが起こる例を紹介する。ここで起こることは、例え、メインの興味ではなくても、ランダム効果として変数に組み込むことで、致命的な解釈の間違いを避けられることを意味している。もちろん、固定効果として組み込んでも同じ効能は得られるが、自由度などの観点からランダム効果としてみなすことで得られるメリットが存在している。

具体的な解析

　では、以下のようなデータが得られているとき、どうやって解析していくかを考える。このデータをまず図示してみる。

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

library(plotn)

library(glmmTMB)

x <- c(3.6, 4.7, 3.5, 4.3, 2.6, 2.7, 3.4, 3.7, 3.2, 4, 0, -2.6, -1.6, 

       -0.7, 0.1, -2, -2.7, -2.7, -2.1, -0.5, -4.2, -2.8, -2.6, -3.4, 

       -0.9, -2.5, -4.7, -0.5, -1.9, -2.5, 4.3, 1.5, 5.7, 4.2, 5.4, 4.6, 

       3.9, 5.2, 3.3, 3.4, 1.2, 0, 1.2, -0.3, 0.7, 0.7, 0.1, 1, 0.3, -0.3,

       2.7, 3.4, 1.2, 3.7, 1.6, 3.5, 3.5, 3.5, 1.3, 1.8, 2.6, 2.9, 2.7, 

       3, 2.8, 2.3, 3.4, 5.1, 2.7, 3.5, 1.4, 0.7, 0.8, 0.1, 1.4, 0.4, 1.7, 

       1.1, -1, 0.9, 3.9, 1.8, 1.5, 1, 1.8, 0.5, 0.4, 0.5, 0.3, 2.5, 3.7,

       0.8, 3.7, 2.9, 4.1, 2.7, 3.2, 2.4, 2.8, 1.3)

y <- c(-7.7, -5.7, -4.2, -5.7, -7.6, -7.5, -6.5, -6.1, -5.9, -4.6, -1.6, 

       -0.6, -1.3, -1.4, -1.7, -1.6, -0.8, -1.3, -1.5, -0.5, -2.3, -0.9, 

       -0.4, -0.6, 2.3, -0.5, -4.3, 2, -0.2, -0.3, -5.8, -12.2, -4.6, -6.8, 

       -4.7, -4.2, -7.8, -5.6, -8.9, -8.5, -1.2, -3.7, -1.6, -4.3, -3.2,

       -2.6, -2, -1.1, -3.1, -4, -5.2, -4.9, -6.8, -3.8, -6.5, -5.1, -4.4,

       -4.3, -6.7, -6.6, -7.2, -6.1, -6.9, -6.8, -6.6, -6.3, -5.6, -4.1, 

       -6.4, -5.4, -2.2, -2.9, -4.4, -3.9, -3.4, -4.1, -3.8, -3.2, -5.4, 

       -4, -2.2, -4.4, -5, -5.1, -4.1, -6.5, -6.6, -6, -5.7, -3.5, -5.7, 

       -6.6, -6, -6.9, -6.4, -4.8, -7.8, -6.1, -6.3, -6.7)

plotn(y ~ x)#図1の描画

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

データを見慣れているとこの時点で、あからさまに嫌な予感がしてくる。もちろん、現実世界のデータはこんなわかりやすくなんてないので、見抜くのは至難の業ではあるが。このデータを何も考慮せず、線形回帰してみよう。

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

summary(lm(y ~ x))

## 

## Call:

## lm(formula = y ~ x)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -7.7804 -1.3572  0.2732  1.2893  5.1240 

## 

## Coefficients:

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

## (Intercept) -3.42238    0.23662 -14.464  < 2e-16 ***

## x           -0.66483    0.08654  -7.682  1.2e-11 ***

## ---

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

## 

## Residual standard error: 1.992 on 98 degrees of freedom

## Multiple R-squared:  0.3759, Adjusted R-squared:  0.3695 

## F-statistic: 59.01 on 1 and 98 DF,  p-value: 1.202e-11

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

なるほど、xの係数の推定値は-0.66程度で、かつ有意である。では、このデータはxは増えると、yが減る、というような負の相関があるのだ……と考えるのは間違いになりえる。さて、後付けだが、このデータ、実は複数地点から取得したものだったのだ。グループの振り分けは以下のとおりである。

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

g <- 10

n <- 10

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

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

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

グループを考慮して図示すると下記の通り。

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

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

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

同じ色＝同じグループの点だけに注目するとわかるが、グループ内ではどう見ても正の相関がある。では、グループをランダム効果としてGLMMを実行してみよう。

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

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 

##    292.7    308.3   -140.3    280.7       94 

## 

## Random effects:

## 

## Conditional model:

##  Groups   Name        Variance Std.Dev. Corr  

##  group    (Intercept) 21.3513  4.6207         

##           x            0.3174  0.5633   -0.34 

##  Residual              0.4347  0.6593         

## Number of obs: 100, groups:  group, 10

## 

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

## 

## Conditional model:

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

## (Intercept)  -6.1201     1.4749  -4.149 3.33e-05 ***

## x             0.9717     0.1945   4.995 5.87e-07 ***

## ---

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

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

すると確かに、xの係数の推定値は0.97と「正になり」、有意である。全体傾向と異なり、群内の傾向はむしろxが増えるとyも増えるという傾向を示している。

　線形回帰とGLMMの結果を図示すると以下のようになる。

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

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

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

b0e_l <- coef(summary(lm(y ~ x)))[1,1]

b1e_l <- coef(summary(lm(y ~ x)))[2,1]

plotn(y ~ x, data = d)

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

yy1 <- b0e + b1e * xx

yy2 <- b0e_l + b1e_l * xx

overdraw(points(xx, yy1, type = "l", col = "red"),

         points(xx, yy2, type = "l", col = "blue"))#図3の描画

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

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

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

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]))

}

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

そもそも、なんでこんなことが起こってしまったのだろうか。データは下記のようなスクリプトから生成した。

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

b0 <- -3

b1 <- 1

s1 <- 1

s2 <- 0.7

is <- 3

ss <- 0.5

sign <- -1

g <- 10

n <- 10

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

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

x <- round(rnorm(n*g, mean = sign*i_sds, sd = s1), digits = 1)

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

y <- rnorm(n = n*g,

           mean = (b1 + s_sds) * 

             (x - rep(tapply(x, group, mean), each = n)) + (b0 + i_sds),

           sd = s2)

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

これを見てもすぐには理解ができないだろう。よく見るとi_sdsはxおよびyの両方でデータの生成に関わっている。また、yの生成過程に関してi_sdsはランダム切片に相当する値であることがわかる。実際の世界ではバックに潜んでいるこの値を知ることはできないが、シミュレーションなので今回生成したデータに関してi_sdsの記録をとっておいてある。それは下記のような値である。この値とxとyの関係を確認してみよう。

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

i_sds <- c(-3.3, -3.3, -3.3, -3.3, -3.3, -3.3, -3.3, -3.3, -3.3, -3.3, 

           1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 2.4, 2.4,

           2.4, 2.4, 2.4, 2.4, 2.4, 2.4, 2.4, 2.4, -3.8, -3.8, -3.8, 

           -3.8, -3.8, -3.8, -3.8, -3.8, -3.8, -3.8, 0.1, 0.1, 0.1, 0.1, 

           0.1, 0.1, 0.1, 0.1, 0.1, 0.1, -2.4, -2.4, -2.4, -2.4, -2.4, 

           -2.4, -2.4, -2.4, -2.4, -2.4, -3.3, -3.3, -3.3, -3.3, -3.3, 

           -3.3, -3.3, -3.3, -3.3, -3.3, -0.8, -0.8, -0.8, -0.8, -0.8, 

           -0.8, -0.8, -0.8, -0.8, -0.8, -1.4, -1.4, -1.4, -1.4, -1.4, 

           -1.4, -1.4, -1.4, -1.4, -1.4, -3, -3, -3, -3, -3, -3, 

           -3, -3, -3, -3)

plotn(x ~ i_sds)

plotn(y ~ i_sds)

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

切片ランダム効果は当然、値が大きければyも大きくなり、正の相関が出る。一方で、切片ランダム効果がxにも負の影響を与えていた。このような状況を交絡coufoundingと呼び、警戒するべき事項の一つである（交絡は別にGLMMに限らず、ANOVAだろうが、GLMだろうが出てくる問題である）。

例えば、日の当たり具合（照度）と植物の大きさの関係をモデル化することを考えよう。今回はいくつかの集団で計測したことにする。直感的予測では、普通の植物なら照度が強ければ成長も早くて大きい植物になるはずだ。ところが、照度が強いのに植物が小さく、逆に弱いのに植物が大きい事例が出てくる。

実はこの例も裏で見えない因子が照度と植物の大きさに影響しているのである。「攪乱からの時間」がいい例となる。草刈りなどの攪乱は、当然だがそこに生えている植生をリセットする行為である。当然、照度は強くなるだろう。そして、リセットされた植生から植物は種子などから成長を開始するから小さい植物ばかりが集まっているはずである。攪乱から時間がたてば、植生も回復し、うっそうとしてきて照度が弱まる。植物は生育が進み大きなものばかりとなっていることだろう。

このような交絡因子の存在によって、データの解釈は捻じ曲げられてしまう。しかし、厄介なことに、このメカニズムを見抜くのは至難の業だ。上記の例だって、攪乱からの時間を、初めから観察していたならまだしも、終わってしまった後から観測するのは困難である。ゆえに、具体的なメカニズムがわからなくても、ランダム効果として観測した集団を組み込むことが、正しい解釈をする上で重要になりえるのである。