階層的なランダム効果が
仮定されるときの考え方
階層的なランダム効果
様々なデータを解析していて、頭を悩ます問題の一つがランダム効果に階層性が認められるときである。例えば、ランダム効果の一つとして、栽培実験を行ったときにいくつかの区画(A, B……)に分けて管理した場合を考えよう。各栽培区画には、一つの鉢に何個体かの植物を混植し、これをいくつも用意した(鉢1、2……)とする。ここで植え付けた植物の結実率を調査するとしよう。
すると区画以外にも、各鉢内に複数の植物が存在し、疑似反復になっている。なので、鉢の効果もランダム効果として取り込むことが可能である。
さてこの時、区画と鉢の関係性に目を向けると、A区の鉢1は、当たり前だがB区に存在しないし、その逆もしかりである。便宜上鉢に1、2……と番号を振っているが、A区の鉢1とB区の鉢1は別物で、それぞれが個別の効果を持っているはずである。見方を変えれば、区画という上位階層に対して、鉢という下位階層が存在するということである。このような状態を鉢が区画にnestされていると表現する。このnest状態は注意をしないと、想定した解析にならないことがある。
Crossed design vs. nested design
通常、2つ以上の要因(固定効果)やランダム効果を考えるとき、例えば下記のような状況が考えられる。
このようなA区とB区の1区が同じ意味を持ち、逆に1区と2区でA区が同じ意味をもつ状況をcrossed designと呼ぶ。
一方、改めて階層的なランダム効果があるときを考えると、A区の鉢1はB区の鉢1と同じ意味を持たない。
ゆえにA区内の鉢にはA区だけの効果、B区の鉢にはB区だけの効果が働くはずだ。このような状況をnested designと呼ぶ。
そもそもランダム効果に限らず、固定効果にもcrossed design/nested designが存在し、今までは暗にcrossed designを仮定して扱ってきた。これは統計そのものというより実験計画法であるので、別の機会に学んでほしいが、階層的なランダム効果は実際のところかなりの頻度で直面する問題なので、ここで取り扱う。
具体的なデータの解析:Crossed design
では、以下のようなデータが得られているとき、どうやって解析していくかを考える。このデータをまず図示してみる。まず、今回のデータはランダム効果のgroup1と2に関してcrossed designとなっている。後述するように、group1と2を定義しているので、各自、データフレーム化した後に中身を確認するとよい。模式図的には下記の通り。
------------------------------------------------------
library(plotn)
library(glmmTMB)
x <- c(7.6, -3.2, -4.7, 12.9, 13.1, -8.9, -7.2, 5.1, -13.8, -2.3, 0.6, 11.2, 11.7, -2.7, 9.8, 1.1, 9.2, 8, 14.4, 5.7, -12.5, -12.7, 7.8, 13.3, 6.9, 6.9, -7.9, -12.3, -7.4, -2.7, -6.4, 13.7, 13.7, 11.3, 14.8, 5.1, -3.2, -2, -8.7, -10.5, 8.9, 13.6, -3.3, -7.8, -1.5, 0, 4.3, -5.7, -2, 13.1, -11.6, -6, -10.9, -9, -11.8, -13.5, 7, -7.4, 0.4, -5.9, 9.2, -5.2, 3.9, 9.2, -2.5, 5.7, 7.5, -5.1, -11.1, -8.2, 1.3, -1.5, -10.9, 7, 0.3, -3.5, 5.4, -10.9, -14.4, 9.1, 8.3, 14.2, -1.8, 1.3, -10.6, -6.6, -13.3, 13.2, -1.2, 9.8, 8.6, -13.2, 10, -0.4, 0.3, -2.4, -7.6, 0.8, -3.7, -12.5, -2, 8, -7.6, 3.5, 10.7, 12.2, 13, 9.3, 5.3, -11.9, 6.3, -7.4, 4.2, -5.3, -15, -12.3, 1.5, 4.8, 6.7, -4.3, -13.3, -5.2, 9.8, -11.4, -4, -9.6, 8, 0.4, 4.5, 2.3, 0.4, -1.3, -6.6, -5.4, 5.8, -13, -9.7, 6.2, 3.4, 11.4, -3.8, -13, -14.1, -7, 4.4, 14.7, -11.7, -1.6, 3.5, -4.9, 3, 6.4, 3.5, 2.6, 7.7, -4.9, -14.1, -3.8, 3.2, -6.1, -8, 14, 10, 3, -12.1, 9.5, 9.5, -3.1, -10.1, 5.5, 3.3, 9.4, -2.4, 13, -3.2, -6.4, -10.9, -4.9, -14.4, 7.8, 10.4, 7.6, 0.4, 10.8, -12.5, -2.7, 11.8, 9.7, 3.6, 2.9, -8.4, 5.4, -4.7, 7.5, -2.5, -12.4, -12.1, 2.6, -13, -14.9, 1.7, 7.1, 14, -8.6, -13.2, -10.8, -9, -9.1, -9, 3.1, -2, 8.1, -11.3, 0, 0.6, 13.3, -2.2, 4.4, 2, -14, -3.3, 13.6, -4.4, -10.6, -13.3, -7.6, -5.5, -4.9, 2.7, -2.4, 12.2, 8.7, 1.7, 14.9, -4.4, 1.6, -10.1, 10.8, 13.6, 6.4, -9.6, -3.4, -13, 8.7, 4.7, -4.8, 10.5, -3.7, 1.1, -13.7)
y <- c(0, 2.9, 5.2, 2, 4.4, 3, 0.6, 4.8, 6.7, 0.5, 7.7, 0.9, 5.9, 5.9, 5.2, 3.7, 4.9, 6.4, 3.6, 7.1, 4.7, 2.7, 0.2, 8.8, 4, 5.6, 5, 3.7, 1, 4.9, 9.7, 0.7, 1.7, -1.3, 4.5, 7.1, 8.2, 6.5, 10.9, 3.6, 4.1, 8.8, 7.5, 1.2, 4.3, 9.4, 7.9, 1.3, 5.3, 1.4, 1.8, 4.2, 0.8, 4.9, 1.2, -0.8, -2.9, 1.3, -3.7, -3.9, -1.2, 0, -5, -0.2, 4.3, -5.1, -3.7, 4.9, 3.8, -2, -2.7, 3.3, -2.1, 2.4, 2.2, 2.9, 0, -5.3, -5.4, 1.8, 3.7, -1.7, -1.1, 1.5, 8.5, 0, 0.7, 0.1, 5.6, 0.9, -4.7, 0.7, 4, -4.7, 0.9, 3.2, 7, -1, -1.8, 2.2, 0.7, -2.7, -4.2, 2.1, -2.4, -2.4, 3.7, -3.3, -3.2, -2, 4.6, -5, -1.2, -1.2, 1.2, 3.1, 0.4, -4.5, 2.4, 2.5, -1.7, -1.1, -2.6, -3.8, 0.2, -2.7, -2.6, 0.5, -0.7, 2.9, 3.5, 5, 8, -1.3, 3.6, 0.9, -0.8, 5.9, 2.7, -2.8, -6.7, -0.7, -3.7, -3.3, -1, -1.6, 1, 2.4, -2.3, 1.8, -2.3, -1, -5.8, 1, -3.2, -0.9, -3.8, -1.4, -1.3, -4.9, 1.6, -6.5, -0.2, -2.9, -2, -5.3, -0.5, -7.3, 1.7, 6.7, -7.8, -4.2, 1.3, -3.5, -1.3, -5.1, -1.4, 1.7, -4.1, -0.1, 1.6, 4.6, 3.2, -0.7, -2.4, -1.3, 0.8, 0.3, -5.6, -0.3, 2.9, -5.7, 4, 0.4, -4.2, -2.4, 2.9, 0.6, -6.8, 2.5, -4.9, -0.4, -3, 0.1, -5.2, 2.2, 2.4, 1.4, -1.6, 0.7, -0.5, -2.3, -4.2, 0.2, -1.4, 0.2, -1.7, -0.9, 1.5, 5.4, 1.6, -2.6, -2.1, 0.8, -6.5, -3, -2.9, -0.6, -0.3, -2.2, -2.7, 0.6, -2.2, 1.7, 5.1, 2, 3.3, -2.8, 4.4, -0.2, 2.3, -2, -2.4, -4.7, -3, 0, -1.4, 0.7, -1, 0.5)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
plotn(y ~ x, data = d, mode = "m", group = "group1")#図1の描画
plotn(y ~ x, data = d, mode = "m", group = "group2")#図2の描画
------------------------------------------------------
図1 データの図示。group1ごとに色を分けた。
図2 データの図示。group2ごとに色を分けた。
xとyはあまり強い関係がなさそうだ。では、このデータに関して、crossed designを想定したランダム効果の組み込み方を紹介する。また、切片項だけにランダム効果を仮定する。下記のように行う。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1) + (1 | group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1283.7 1301.3 -636.8 1273.7 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group1 (Intercept) 4.3913 2.0955
## group2 (Intercept) 0.8647 0.9299
## Residual 8.6890 2.9477
## Number of obs: 250, groups: group1, 5; group2, 5
##
## Dispersion estimate for gaussian family (sigma^2): 8.69
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51571 1.04212 0.495 0.621
## x -0.04154 0.02276 -1.825 0.068 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
2つのランダム効果を(1|group1)+(1|group2)という和の形で入力した。現在のデータフレーム内のgroup1と2の入力状態で、上記のような記述をすればcrossedを想定することになる。group1、2、それぞれに関して標準偏差2.10、0.93の正規分布から生成されていると推定している。
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"))
plotn(y ~ x, data = d, mode = "m", group = "group1")#図4の描画
for(i in 1:length(unique(d$group1))){
tmp <- d[d$group1 == unique(d$group1)[i],]
xx <- seq(min(tmp$x), max(tmp$x), length = 200)
yy <- (b0e + ranef(fit1)$cond$group1[i,1] +
mean(unlist(ranef(fit1)$cond$group2))) + b1e * xx
overdraw(points(xx, yy, type = "l", col = .default_col[i]))
}
plotn(y ~ x, data = d, mode = "m", group = "group2")#図5の描画
for(i in 1:length(unique(d$group2))){
tmp <- d[d$group2 == unique(d$group2)[i],]
xx <- seq(min(tmp$x), max(tmp$x), length = 200)
yy <- (b0e + ranef(fit1)$cond$group2[i,1] +
mean(unlist(ranef(fit1)$cond$group1))) + b1e * xx
overdraw(points(xx, yy, type = "l", col = .default_col[i]))
}
------------------------------------------------------
図3 データの図示。回帰曲線も描いた。
図4 データの図示。group1ごとに回帰曲線も描いた。
図5 データの図示。group2ごとに回帰曲線も描いた。
さて、もしこれをnestedな想定で解析するとどうなるだろうか。下記のようにすればnestedな解析になる。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1/group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1293.6 1311.2 -641.8 1283.6 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group2:group1 (Intercept) 0.5254 0.7248
## group1 (Intercept) 4.0917 2.0228
## Residual 8.9869 2.9978
## Number of obs: 250, groups: group2:group1, 25; group1, 5
##
## Dispersion estimate for gaussian family (sigma^2): 8.99
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51922 0.93561 0.555 0.579
## x -0.03036 0.02340 -1.298 0.194
------------------------------------------------------
今回は、(1|group1/group2)という項を足している。中身は上の階層/下の階層、となっており、この表記ならgroup1にgroup2がnestされていることを想定している。解析結果も変わったことがわかるだろう。注目するべきは、Random effectsにgroup2の代わりにgroup2:group1という項が登場していることであり、これがgroup1ごとのgroup2の効果を意味している。それを反映して、group2:group1のグループ数も5×5 = 25になっている。
ちなみにこれらの表記法は重要でgroup1/group2 = group1 + group2:group1であることを意味している。固定効果として組み込むことを考えるときも、formulaをx + group1/group2やx + group1 + group2:group1とすれば、nested ANOVA、nested glmなどを実行可能である。なお、このデータであればAICを見ると、crossed想定のほうが予測性が高いことがわかる。
今回のデータはcrossedを想定して、下記のように生成した。
------------------------------------------------------
b0 <- 2
b1 <- 0
s <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rep(rnorm(g2, mean = 0, sd = is2), each = n), times = g1)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s)
------------------------------------------------------
その時の、ランダム切片効果を記録してある。下記のとおりである。この時のgroup1、2のランダム効果の関係性を確認してみよう。
------------------------------------------------------
i_sds1 <- c(3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, 3.6, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.2, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4, -2.4)
i_sds2 <- c(-0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, -1.1, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2, -1.2)
plotn(i_sds2 ~ i_sds1)#図6の描画
------------------------------------------------------
図6 横軸=group1、縦軸=group2のランダム効果
データの形としては、このようなランダム効果の関係となっていればcrossed designである。
では、crossed designのデータの時、crossedあるいはnested designを想定した解析を行うとどうなるかシミュレーションしてみよう。
------------------------------------------------------
b0 <- 2
b1 <- 0
s <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
p1_c <- NULL
p2_c <- NULL
b0e_c <- NULL
b1e_c <- NULL
is_e1_c <- NULL
is_e2_c <- NULL
p1_n <- NULL
p2_n <- NULL
b0e_n <- NULL
b1e_n <- NULL
is_e1_n <- NULL
is_e2_n <- NULL
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
for (i in 1:1000){
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rep(rnorm(g2, mean = 0, sd = is2), each = n), times = g1)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity")) #crossed design
fit2 <- summary(fit1)
fit3 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity")) #nested design
fit4 <- summary(fit3)
p1_c <- c(p1_c, coef(fit2)$cond[1,4])
p2_c <- c(p2_c, coef(fit2)$cond[2,4])
b0e_c <- c(b0e_c, coef(fit2)$cond[1,1])
b1e_c <- c(b1e_c, coef(fit2)$cond[2,1])
is_e1_c <- c(is_e1_c, as.numeric(fit2$varcor$cond$group1[1,1]))
is_e2_c <- c(is_e2_c, as.numeric(fit2$varcor$cond$group2[1,1]))
p1_n <- c(p1_n, coef(fit4)$cond[1,4])
p2_n <- c(p2_n, coef(fit4)$cond[2,4])
b0e_n <- c(b0e_n, coef(fit4)$cond[1,1])
b1e_n <- c(b1e_n, coef(fit4)$cond[2,1])
is_e1_n <- c(is_e1_n, as.numeric(fit4$varcor$cond$group1[1,1]))
is_e2_n <- c(is_e2_n, as.numeric(fit4$varcor$cond$`group2:group1`[1,1]))
}
------------------------------------------------------
では、この時の危険率、検出力などを確認してみよう。
------------------------------------------------------
#crossed design
sum(p1_c < 0.05, na.rm = T)
## [1] 396
sum(p2_c < 0.05, na.rm = T)
## [1] 47
sum(!is.na(p1_c))
## [1] 1000
sum(!is.na(p2_c))
## [1] 1000
#nested design
sum(p1_n < 0.05, na.rm = T)
## [1] 447
sum(p2_n < 0.05, na.rm = T)
## [1] 45
sum(!is.na(p1_n))
## [1] 1000
sum(!is.na(p2_n))
## [1] 1000
#crossed design
histn(p1_c[!is.na(p1_c)], xlab = "P value", ylab = "頻度")#図7の描画
histn(p2_c[!is.na(p2_c)], xlab = "P value", ylab = "頻度")#図8の描画
#nested design
histn(p1_n[!is.na(p1_n)], xlab = "P value", ylab = "頻度")#図9の描画
histn(p2_n[!is.na(p2_n)], xlab = "P value", ylab = "頻度")#図10の描画
------------------------------------------------------
図7 crossed designの切片の係数のP値のヒストグラム
図8 crossed designのxの係数のP値のヒストグラム
図9 nested designの切片の係数のP値のヒストグラム
図10 nested designのxの係数のP値のヒストグラム
危険率などに関して、どちらの解析でも大差がないことがわかる。次は係数の推定値を確認してみる。
------------------------------------------------------
#crossed design
mean(b0e_c[!is.na(p1_c)])
## [1] 2.028986
mean(b1e_c[!is.na(p2_c)])
## [1] -0.0007117002
var(b0e_c[!is.na(p1_c)])
## [1] 2.039897
var(b1e_c[!is.na(p2_c)])
## [1] 0.000468113
#nested design
mean(b0e_n[!is.na(p1_n)])
## [1] 2.028979
mean(b1e_n[!is.na(p2_n)])
## [1] -0.0006816258
var(b0e_n[!is.na(p1_n)])
## [1] 2.04002
var(b1e_n[!is.na(p2_n)])
## [1] 0.000479986
#crossed design
histn(b0e_c[!is.na(p1_c)], xlab = "b0", ylab = "頻度")#図11の描画
histn(b1e_c[!is.na(p2_c)], xlab = "b1", ylab = "頻度")#図12の描画
#nested design
histn(b0e_n[!is.na(p1_n)], xlab = "b0", ylab = "頻度")#図13の描画
histn(b1e_n[!is.na(p2_n)], xlab = "b1", ylab = "頻度")#図14の描画
------------------------------------------------------
図11 crossed designの切片の係数のヒストグラム
図12 crossed designのxの係数のヒストグラム
図13 nested designの切片の係数のヒストグラム
図14 nested designのxの係数のヒストグラム
係数推定値も大差はない。ランダム効果の分散も確認してみる。
------------------------------------------------------
#crossed design
mean(is_e1_c[!is.na(p1_c)])
## [1] 7.320218
mean(is_e2_c[!is.na(p2_c)])
## [1] 0.9382674
#nested design
mean(is_e1_n[!is.na(p1_n)])
## [1] 6.967927
mean(is_e2_n[!is.na(p2_n)])
## [1] 1.013542
#crossed design
histn(is_e1_c[!is.na(p1_c)], ylab = "頻度")#図15の描画
histn(is_e2_c[!is.na(p2_c)], ylab = "頻度")#図16の描画
#nested design
histn(is_e1_n[!is.na(p1_n)], ylab = "頻度")#図17の描画
histn(is_e2_n[!is.na(p2_n)], ylab = "頻度")#図18の描画
------------------------------------------------------
図15 crossed designの切片のランダム効果の分散のヒストグラム
図16 crossed designのxの係数のランダム効果の分散のヒストグラム
図17 nested designの切片のランダム効果の分散のヒストグラム
図18 nested designのxの係数のランダム効果の分散のヒストグラム
こちらも推定に大きな差はなさそうだ。最後、AICを確認しておこう。
------------------------------------------------------
#crossed design
mean(aic_c[!is.na(p1_c)], na.rm = T)
## [1] 1290.751
#nested design
mean(aic_n[!is.na(p1_n)], na.rm = T)
## [1] 1297.111
#crossed design
histn(aic_c[!is.na(p1_c)], ylab = "頻度")#図19の描画
#nested design
histn(aic_n[!is.na(p1_n)], ylab = "頻度")#図20の描画
------------------------------------------------------
図19 crossed designのAICのヒストグラム
図20 nested designのAICのヒストグラム
AICを見ると、crossed designのほうが予測性が良いことが変わる。その意味ではより適切なモデルはデータ生成過程通り、crossed designである。
具体的なデータの解析:nested design
次に、以下のようなデータが得られているとき、どうやって解析していくかを考える。今回のデータはgroup2がgroup1にnestされていることを想定する。つまり、group2でaというラベルになっていても、group1におけるラベルによってその意味は変わる。ということは単純にgroup2のラベルだけでグループを振り分けた図にしても意味はない。ここでは、group1についてだけ振り分けた図を描画する。group2に関しては各自考えてみよう(ヒント:group1と2を組み合わせた識別因子があれば実現可能)。
------------------------------------------------------
x <- c(10.1, 4.2, 10.4, -11.4, -7.3, 10.2, 3.6, -10.2, 12.8, 3, -4.7, -13.7, -11.1, 5.4, 11, -0.2, 14.7, 5.6, -10.5, 14.5, 9.5, 0.4, 8.7, -12.8, 12.8, 3.2, -13.9, 8.5, -4.9, -9.5, -10.3, 2.1, -8.7, -11.2, 8.5, -6.8, -2.4, -2.3, 9.5, 1.5, -0.8, 2.6, -2.3, -11.1, -7.7, -8.2, -1.1, 12.6, 11.3, -6.8, -4.8, 2.6, -13.2, -5.6, -4.9, -12.5, 12.4, -4.8, 3.2, -7, -11.3, -9.1, 13.1, 5.8, 6.8, -4.7, 9.6, 13.6, -11.1, -0.6, 14.6, -6.2, 5.8, -4.5, 10.6, -5.2, -8.9, -5.2, -12.6, 7, -9.2, 0.9, -1.3, -1.6, -0.8, 1.6, -7.9, 11, 8.3, 3, 1, 1.9, 0.7, -4, -11, 12.1, -3.1, 14.3, 11.7, -1.6, -5, 1.2, -6.2, 9.8, 9.7, -5.1, -11.8, -0.6, -12, 5.4, 11.6, 9.6, 14.1, -4.9, 4.6, 13.1, -13.2, -4.9, -1.7, -12.4, -12.3, 3.2, -6.6, -6.9, -11.3, 4.6, -6.4, 12, 6.3, 7.6, 2.1, -10.7, -6.7, 3.4, -0.4, -2.5, -3.1, 3.3, 0, 3.7, -8.7, 1.6, 5.6, 7.1, 13.7, 7.6, -6.2, -14.1, -9.2, 14.5, 12.8, -0.5, -3.1, -3.2, -4.2, -6.5, 12.1, 0.2, 9.9, 13.7, -0.3, 5.1, -6.7, 4, -13.6, 10.7, 5.7, -11.9, -0.4, 14.2, 2.8, -0.5, 5.2, 7.4, -13.4, 7.5, 8.8, 11.5, -10.6, -5, 1.5, -14.1, 1.6, -14.5, -3.1, -14.9, 8.4, 11, 3.6, -5.1, 2.7, 2.1, 6.2, -13.6, 4.4, -14, -12.8, -8, 0, 6.4, -10.3, 12.4, 12, 3.7, 12.7, 7.9, -8.3, -10, 3.3, -3.3, -9.9, -2.3, -11.9, 4, 0.2, 14.8, 8.7, -10.5, 1.9, -2.5, 4.3, -14.5, -7.7, -14.7, -2.2, -3.7, 6, 13, 1.2, 7.7, 11.5, 6.4, 4.2, -2.6, -8.9, 11.2, -14.9, -8.6, 6.9, -7.4, 1.3, -8.8, -11.5, 13.1, -1, 13.7, -13.5, 0.1, -12.1, 7.7)
y <- c(-4.6, -0.1, -0.6, 1.4, 2.4, -3.8, 1.7, -3.5, -5.1, -1.5, -2.3, -7.6, -6.6, 0.2, -4.2, 0.2, -3.9, -1.3, 1, -2.5, -11.5, -5.3, -9.6, -6.2, -11, 0.6, -4.1, -5.5, -4.7, -6.5, -7.1, 0.1, 0.7, -1.3, -1.6, 1.4, -2.1, 0.2, -4.2, 0.9, -6.9, -0.9, 1.9, -4.9, -1.4, -4, 0.6, -6.9, 1.2, -5.3, -2.1, 5.5, -2.1, 1.5, 3.7, 1.6, -1.8, 0.5, -0.3, 0.8, 4.9, 8.8, 3, 2.2, 2.3, 6.7, 1.8, -2.4, 3.7, 3.7, 4.2, 6.6, 3.1, 3.8, 2.5, 1.8, 3.6, 2.9, 1.3, 0.6, 3.8, 7.7, 3.6, 6.4, 3.2, 8.5, 1.8, 5.2, 0.5, 7.6, 3.8, 7.8, 2.1, 6.9, 5.7, 3.8, 6.2, 6.2, 0.6, 1.7, -1.6, -3.9, -5.4, -2.3, -2.6, -5.6, -6.1, -7.4, -2.8, -6.6, -3.9, -9.4, -2.9, -1.1, -12.8, -7.9, -4.7, 2.6, -6.2, -0.2, 2.3, -6.3, -5.8, -3.8, -0.1, -1.8, -0.8, -0.3, 2.6, 1.1, -2.2, -8.6, 1.1, -2.9, -2.7, -3.1, -3.5, -2.4, -2.1, -0.1, -1, -7.4, -1.2, 0, 0.7, 4.7, -0.7, -6.2, -6.7, -3.6, 0.1, 5.6, 2.2, 0.2, 3, -1, 1.4, -0.9, 4.3, 3, -4, -4.9, 1.2, -2.9, 3.1, 0.6, -2, 0.2, 7.9, -5.2, 0.1, -5, 0.6, -2.1, -4.4, -0.1, 0.2, -3.6, 0.5, -3.4, -1, -1.7, -1.9, -4, 2, 0.3, -2.3, -7.6, 0.7, -0.9, 2.5, 0.3, 1, -1.3, 1.7, 0, -2.8, 6, -5.8, -2.6, 2.3, 1.7, -5.2, 0.2, -4.9, -1.6, 1, -2.7, 1.1, -2.1, 1.5, 5.4, -4.3, 0.3, -5.7, -2.7, 0.2, -6.1, -6.1, -5.6, -1.9, -0.5, -1, -7.7, 1.8, 0.6, -7.3, 1.1, -2, -5.6, -6.2, -3.2, -0.8, -3.7, 1.5, -3.9, -1.1, -2.3, -1.1, -1.6, -5.5, 1.1, -2.4, -3.5, -5, -4.5, 2.4, -1.1, 4.2, -2.5)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
plotn(y ~ x, data = d, mode = "m", group = "group1")#図21の描画
------------------------------------------------------
図19 データの図示。group1ごとに色を分けた。
xとyはあまり強い関係がなさそうだ。まあ、お察しの通り、そういうデータを生成した。では、このデータに関して、nested designを想定したランダム効果の組み込み方を紹介する。また、切片項だけにランダム効果を仮定する。下記のように行う。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1/group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1304.9 1322.5 -647.4 1294.9 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group2:group1 (Intercept) 1.473 1.214
## group1 (Intercept) 4.999 2.236
## Residual 8.986 2.998
## Number of obs: 250, groups: group2:group1, 25; group1, 5
##
## Dispersion estimate for gaussian family (sigma^2): 8.99
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.99397 1.04630 -0.950 0.3421
## x -0.04387 0.02247 -1.953 0.0508 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
上記で紹介したように、group1に2がnestされていることを想定した、(1|group1/group2)という項を足している。group1は標準偏差2.24、group1と2の組み合わせ効果は標準偏差1.21の正規分布から生成されていると推定している。
GLMMの結果を図示すると以下のようになる。今回もgroup1のみ図示することとする。
------------------------------------------------------
b0e <- coef(fit2)$cond[1,1]
b1e <- coef(fit2)$cond[2,1]
plotn(y ~ x, data = d)#図22の描画
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 = "group1")#図23の描画
for(i in 1:length(unique(d$group1))){
tmp <- d[d$group1 == unique(d$group1)[i],]
xx <- seq(min(tmp$x), max(tmp$x), length = 200)
yy <- (b0e + ranef(fit1)$cond$group1[i,1] +
mean(unlist(ranef(fit1)$cond$group2))) + b1e * xx
overdraw(points(xx, yy, type = "l", col = .default_col[i]))
}
------------------------------------------------------
図22 データの図示。回帰曲線も描いた。
図23 データの図示。group1ごとに回帰曲線も描いた。
さて、もしこれをcrossedな想定で解析するとどうなるだろうか。下記のようにすればcrossedな解析になる。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1) + (1 | group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1315.9 1333.5 -653.0 1305.9 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group1 (Intercept) 5.28518 2.2990
## group2 (Intercept) 0.07338 0.2709
## Residual 10.12274 3.1816
## Number of obs: 250, groups: group1, 5; group2, 5
##
## Dispersion estimate for gaussian family (sigma^2): 10.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.99352 1.05461 -0.942 0.346
## x -0.03550 0.02352 -1.509 0.131
------------------------------------------------------
結果が変わったことがわかるだろう。注目するべきは、group2の標準偏差がnested過少推定となっていることである。その詳細な理由は後述する。
さて、一つテクニックとして、上記のようなcrossed想定の表記でもnestedな解析をすることができる。それは、nestされているグループを上位の階層ごとに違う識別子を用いることだ。例えば、今回のデータではgroup2がnestされており、その識別子はa ~ eの5つである。そうではなく、group1のAのgroup2にはa ~ e、Bのgroup2にはf ~ jみたいに振り分ける。
つまり、下記のようなgroup2を用意して解析すると、nestedな解析と同じ結果になる。
------------------------------------------------------
group2 <- rep(letters[1:(g1*g2)], each = n)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1) + (1 | group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1304.9 1322.5 -647.4 1294.9 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group1 (Intercept) 4.999 2.236
## group2 (Intercept) 1.473 1.214
## Residual 8.986 2.998
## Number of obs: 250, groups: group1, 5; group2, 25
##
## Dispersion estimate for gaussian family (sigma^2): 8.99
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.99397 1.04630 -0.950 0.3421
## x -0.04387 0.02247 -1.953 0.0508 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
今回のデータはnestedを想定して、下記のように生成した。group2のランダム効果は平均0、標準偏差1の共通の正規分布から生成されるが、crossedの時と異なり、異なるgroup1ごとに異なる値が生成されるようにしている。
------------------------------------------------------
b0 <- 2
b1 <- 0
s <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rnorm(g1*g2, mean = 0, sd = is2), each = n)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s)
------------------------------------------------------
その時の、ランダム切片効果を記録してある。下記のとおりである。この時のgroup1、2のランダム効果の関係性を確認してみよう。
------------------------------------------------------
i_sds1 <- c(-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4)
i_sds2 <- c(0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3)
plotn(i_sds2 ~ i_sds1)#図24の描画
------------------------------------------------------
図24 横軸=group1、縦軸=group2のランダム効果
データの形としては、このようなランダム効果の関係となっていればnested designである。今回の場合、すべてのgroup1におけるgroup2のランダム効果の平均は0となっている。後述するようにnestな状況では、これ以外にもgroup1ごとにgroup2の平均が違う場合も考えられる。
では、nested designのデータの時、crossedあるいはnested designを想定した解析を行うとどうなるかシミュレーションしてみよう。
------------------------------------------------------
b0 <- 2
b1 <- 0
s <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
p1_c <- NULL
p2_c <- NULL
b0e_c <- NULL
b1e_c <- NULL
is_e1_c <- NULL
is_e2_c <- NULL
p1_n <- NULL
p2_n <- NULL
b0e_n <- NULL
b1e_n <- NULL
is_e1_n <- NULL
is_e2_n <- NULL
aic_c <- NULL
aic_n <- NULL
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
for (i in 1:1000){
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rnorm(g1*g2, mean = 0, sd = is2), each = n)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル
fit2 <- summary(fit1)
fit3 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル
fit4 <- summary(fit3)
p1_c <- c(p1_c, coef(fit2)$cond[1,4])
p2_c <- c(p2_c, coef(fit2)$cond[2,4])
b0e_c <- c(b0e_c, coef(fit2)$cond[1,1])
b1e_c <- c(b1e_c, coef(fit2)$cond[2,1])
is_e1_c <- c(is_e1_c, as.numeric(fit2$varcor$cond$group1[1,1]))
is_e2_c <- c(is_e2_c, as.numeric(fit2$varcor$cond$group2[1,1]))
p1_n <- c(p1_n, coef(fit4)$cond[1,4])
p2_n <- c(p2_n, coef(fit4)$cond[2,4])
b0e_n <- c(b0e_n, coef(fit4)$cond[1,1])
b1e_n <- c(b1e_n, coef(fit4)$cond[2,1])
is_e1_n <- c(is_e1_n, as.numeric(fit4$varcor$cond$group1[1,1]))
is_e2_n <- c(is_e2_n, as.numeric(fit4$varcor$cond$`group2:group1`[1,1]))
aic_c <- c(aic_c, fit2$AICtab[1])
aic_n <- c(aic_n, fit4$AICtab[1])
}
------------------------------------------------------
では、この時の危険率、検出力などを確認してみよう。
------------------------------------------------------
#crossed design
sum(p1_c < 0.05, na.rm = T)
## [1] 419
sum(p2_c < 0.05, na.rm = T)
## [1] 55
sum(!is.na(p1_c))
## [1] 1000
sum(!is.na(p2_c))
## [1] 1000
#nested design
sum(p1_n < 0.05, na.rm = T)
## [1] 433
sum(p2_n < 0.05, na.rm = T)
## [1] 51
sum(!is.na(p1_n))
## [1] 1000
sum(!is.na(p2_n))
## [1] 1000
#crossed design
histn(p1_c[!is.na(p1_c)], xlab = "P value", ylab = "頻度")#図25の描画
histn(p2_c[!is.na(p2_c)], xlab = "P value", ylab = "頻度")#図26の描画
#nested design
histn(p1_n[!is.na(p1_n)], xlab = "P value", ylab = "頻度")#図27の描画
histn(p2_n[!is.na(p2_n)], xlab = "P value", ylab = "頻度")#図28の描画
------------------------------------------------------
図25 crossed designの切片の係数のP値のヒストグラム
図26 crossed designのxの係数のP値のヒストグラム
図27 nested designの切片の係数のP値のヒストグラム
図28 nested designのxの係数のP値のヒストグラム
危険率などに関して、どちらの解析でも大差がないことがわかる。次は係数の推定値を確認してみる。
------------------------------------------------------
#crossed design
mean(b0e_c[!is.na(p1_c)])
## [1] 1.932659
mean(b1e_c[!is.na(p2_c)])
## [1] -0.0002694504
var(b0e_c[!is.na(p1_c)])
## [1] 1.764986
var(b1e_c[!is.na(p2_c)])
## [1] 0.0005595486
#nested design
mean(b0e_n[!is.na(p1_n)])
## [1] 1.932627
mean(b1e_n[!is.na(p2_n)])
## [1] -3.659782e-05
var(b0e_n[!is.na(p1_n)])
## [1] 1.76492
var(b1e_n[!is.na(p2_n)])
## [1] 0.0005318871
#crossed design
histn(b0e_c[!is.na(p1_c)], xlab = "b0", ylab = "頻度")#図29の描画
histn(b1e_c[!is.na(p2_c)], xlab = "b1", ylab = "頻度")#図30の描画
#nested design
histn(b0e_n[!is.na(p1_n)], xlab = "b0", ylab = "頻度")#図31の描画
histn(b1e_n[!is.na(p2_n)], xlab = "b1", ylab = "頻度")#図32の描画
------------------------------------------------------
図29 crossed designの切片の係数のヒストグラム
図30 crossed designのxの係数のヒストグラム
図31 nested designの切片の係数のヒストグラム
図32 nested designのxの係数のヒストグラム
係数推定値も大差はない。ランダム効果の分散も確認してみる。
------------------------------------------------------
#crossed design
mean(is_e1_c[!is.na(p1_c)])
## [1] 7.507709
mean(is_e2_c[!is.na(p2_c)])
## [1] 0.1987615
#nested design
mean(is_e1_n[!is.na(p1_n)])
## [1] 7.278962
mean(is_e2_n[!is.na(p2_n)])
## [1] 1.031972
#crossed design
histn(is_e1_c[!is.na(p1_c)], ylab = "頻度")#図33の描画
histn(is_e2_c[!is.na(p2_c)], ylab = "頻度")#図34の描画
#nested design
histn(is_e1_n[!is.na(p1_n)], ylab = "頻度")#図35の描画
histn(is_e2_n[!is.na(p2_n)], ylab = "頻度")#図36の描画
------------------------------------------------------
図33 crossed designの切片のランダム効果の分散のヒストグラム
図34 crossed designのxの係数のランダム効果の分散のヒストグラム
図35 nested designの切片のランダム効果の分散のヒストグラム
図36 nested designのxの係数のランダム効果の分散のヒストグラム
ここで初めて、差が生じた。nestedなデータをcrossedな解析を行うと、group2の分散が過少推定となる。
原因の一端を解説しよう。crossedな解析を行うということは、各group1に存在するgroup2を区別しないということになる。例えば、group2のaに関しては(group1, group2) = (A, a), (B, a), (C, a), (D, a), (E, a)はランダム切片を共有するということだから、推定値としてはこれらの平均となる。真のデータ生成過程ではgroup2のaは平均0、分散1^2の正規分布から、group1のA ~ Eで独立に生じている。つまり、group2のa ~ eのランダム切片の推定値はそれぞれ平均0、分散1の正規分布からサンプルサイズ5のデータを生成して平均をとった時の分布に従うと考えてよいから、中心極限定理よりその平均値が従う正規分布の分散は1/5 = 0.2程度になると予測できる。なお、実際の中身はこんなシンプルではないようで、今回は予測からよく一致した結果になっているが、group1とgroup2の推定は互いに影響を与え合うため、片方の分散の効果をもう片方に押し付けてしまうみたいなことも生ずる。
ちなみにnestedな解析の場合のgroup2の分散推定値は不偏分散の性質から1^2 × 24/25 = 0.96となることが期待され、その期待通りとなっている(というか、サンプルサイズが大きいので区別つかないだろう)。
最後、AICを比較してみよう。
------------------------------------------------------
#crossed design
mean(aic_c[!is.na(p1_c)], na.rm = T)
## [1] 1305.463
#nested design
mean(aic_n[!is.na(p1_n)], na.rm = T)
## [1] 1299.941
#crossed design
histn(aic_c[!is.na(p1_c)], ylab = "頻度")#図37の描画
#nested design
histn(aic_n[!is.na(p1_n)], ylab = "頻度")#図38の描画
------------------------------------------------------
図37 crossed designのAICのヒストグラム
図38 nested designのAICのヒストグラム
AICを見ると、nested designのほうが予測性が良いことが変わる。その意味ではより適切なモデルはデータ生成過程通り、nested designである。
最後に、もう1パターン、nestedなデータを扱ってみよう。今回も、group1にgroup2がnestされていることを想定する。今回はデータの図示を飛ばして、淡々と解析する。
------------------------------------------------------
x <- c(-10.1, 7.7, -13, -6.6, 11.2, -14, -6.4, 14, -3.3, 3.6, -14.1, -7.8, -0.2, 5.7, -1.5, -14.8, 6, -13.5, 5.7, 3.6, -9.2, 8, 3.2, -9.5, -7.5, -8.8, -5.3, 1.3, 11.4, -3.3, -13, -3.8, 14.5, 9.7, 10, -0.9, -10.4, -14.4, 9.6, 9.7, 6.6, 14.3, 1.6, 0.1, -7.8, -8.3, 3.9, 1.8, -11.8, 0, 3, 2.6, 6.1, 8.3, -5.8, 0.7, 1.3, -4.7, 11.4, -9, 12, -0.4, -14.2, 5.2, 10.9, -7.3, -9.3, 13.8, 6.5, 5.6, -2.8, -1.5, 14.5, 11.2, -0.1, -6.4, 3.4, 9.1, 9.5, 14.1, 10.3, 10.2, -4.5, 6.9, -7.8, -1.5, -6.2, -6.4, 6.8, 1.9, -2.4, 0.9, -6.6, 5, -2.8, -3.7, 6.2, 11.8, -1.1, -6.6, -5.9, 13, 9.7, -10.5, -2.1, -1.8, 4.7, -3.3, -3.3, -0.6, 3.9, -0.2, 12.4, -6.6, 1.8, 12.8, 10.3, 12.5, 7.8, -13.4, -13.6, -6, -8.7, -7.4, -4.3, -11.7, -7.5, 0.4, 1.3, 8.7, -5.3, -4.6, -7.3, -4, 13.7, 5.3, -5.5, -2.2, -7.9, -9.5, -10.3, -4.7, 12.5, 14.8, 2.5, 4.8, -13.1, -13.7, 6.4, 12.3, -0.2, 1.4, 9.4, 3.3, -3, -5.6, -3.4, -8, -3, 5.2, -12.8, 5.7, -3, 1.3, 12.9, -14.8, -14.3, 11.6, 4.7, 14.7, 14.2, 7.8, 10, -5.6, -4.3, 5.3, -6.9, 4.6, 1.7, 5.3, 4.5, -8.1, -8, 2.3, 9.5, 8.5, -11.7, -7.4, 4.5, 0.9, -2.1, 6.1, 2.3, -10.7, -1, -1, -4.8, 5.9, 11.2, -5, 5.5, 9.3, 8.1, 7.1, 10.4, -6.7, -10.9, -3.6, -12.5, 0.3, -0.8, -10, -9, 11.5, -13.6, -5.2, 0.7, -10.3, -7.1, -14, -4.9, -4.3, 10.3, -12.9, -1, -2.8, 11.1, -6.8, 11.6, -3.8, -3.1, -7.4, -9.1, 5.5, -5.2, 8, 0.3, -13.3, 5.4, 8.6, -6.5, -5.4, -8.4, -9.7, 6.8, -1.6, 14.2, -2.3, 1.2, -12.3)
y <- c(6.5, 1.5, -2.5, 0.6, 1.8, 5, 1.7, 5, 5.1, 4, 0.5, -1.9, 5.4, -5.7, 2.1, 1.2, -2.8, -1.2, 3.1, 3.1, 1.8, 0.8, -0.8, 1.7, 0.7, -1.9, -6.3, -6.6, -5.1, 0.2, -1.8, -4, 2.1, 2.1, -5, -6.9, 0.2, -0.9, -4.7, -5.9, 1.8, 1.1, -3, -2.9, -3.7, -0.1, -1, 2.1, -0.3, -3.5, -5.4, -10.6, -8.8, -3.6, -1.3, -3.8, 1.6, -6, -1.7, -2.6, 0.7, 2.2, -2.8, 0.3, 3.1, -0.2, -4.6, 4.1, 0.2, -3, -1.5, -2.2, -10.6, -2.5, -5.9, -5.4, -0.5, -1.6, -2.8, -1.3, 2.3, -2.4, -1, -6, -7, -4.7, 0.7, -1.9, -2.7, -6.2, -4.5, -4.4, -4.7, -11.5, 1, -2.5, -5.4, -4.2, -3.3, -5.6, 0.4, -2.4, 2.7, -3.5, -2, 3.5, -1.2, 3, 5, -1.8, 2.1, -0.8, 5.3, -0.3, -4.7, 4.7, 1.9, 2.2, 1.4, 1.3, -1.5, -3.7, 0, 2.9, -3.7, -7.9, 2.3, 1, -0.6, 4.7, -4.5, 1.8, 3.2, -2, -0.5, -4.6, 0.6, -3.2, 2.3, -4.5, 3.4, 1.3, 2.8, -2.7, -3, 3.4, -2, 3.3, -5, 0.8, 3.2, 2.5, 2.5, 11.9, 5.7, 3.4, 9.2, 6.1, 8.4, 6.6, 7.6, 3, 8.3, 12.7, 4.5, 6.2, 5.1, 4, 2.7, 10.9, 9.4, 8.3, 7, 5.4, 1.6, 0.6, 13, 12.3, 7.5, 7.6, 6.3, 7.4, 5.2, 6.3, 9.2, 4.4, 6, 9.2, 12.4, 12.1, 6, 6.1, 3.8, 9.5, 4.2, 7.4, 6.4, 8.2, 8.9, 9.6, 2, 7.3, 5.9, 9.8, -1.7, 7, 4.3, 1.9, 1.7, 4.5, -0.1, 6.5, 4.8, 5.4, 1.8, 10.2, 2.7, -0.7, 7.2, -1.4, 0.8, 7.5, 5.8, 3.8, 4.2, 1.1, 5.2, 6.6, 1.3, 5.9, 7.6, 4.2, 4.1, 9, 8.5, 6, 4.9, 6.8, 9.3, 9.3, 7.1, 7.2, 6.4, 10.5, 4.2, 2, 9.8, 6.6, 6.9, 5.2)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
------------------------------------------------------
では、このデータをnested designで解析しよう。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1/group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1318.1 1335.7 -654.0 1308.1 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group2:group1 (Intercept) 1.299 1.140
## group1 (Intercept) 13.286 3.645
## Residual 9.385 3.064
## Number of obs: 250, groups: group2:group1, 25; group1, 5
##
## Dispersion estimate for gaussian family (sigma^2): 9.39
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.74849 1.65729 1.055 0.291
## x 0.02553 0.02459 1.038 0.299
------------------------------------------------------
group1は標準偏差3.65、group1と2の組み合わせ効果は標準偏差1.14の正規分布から生成されていると推定している。
さて、もしこれをcrossedな想定で解析するとどうなるだろうか。下記のようにすればcrossedな解析になる。
------------------------------------------------------
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1) + (1 | group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1326.9 1344.5 -658.5 1316.9 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group1 (Intercept) 1.354e+01 3.6791067
## group2 (Intercept) 7.391e-08 0.0002719
## Residual 1.044e+01 3.2318491
## Number of obs: 250, groups: group1, 5; group2, 5
##
## Dispersion estimate for gaussian family (sigma^2): 10.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.74901 1.65800 1.055 0.291
## x 0.02850 0.02542 1.121 0.262
------------------------------------------------------
やはり結果は変わることがわかる。そして、group2の標準偏差はやはり過少推定されている。
下記のようなgroup2を用意して解析すると、nestedな解析と同じ結果になる。
------------------------------------------------------
group2 <- rep(letters[1:(g1*g2)], each = n)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity"))
fit2 <- summary(fit1)
fit2
## Family: gaussian ( identity )
## Formula: y ~ x + (1 | group1) + (1 | group2)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 1318.1 1335.7 -654.0 1308.1 245
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group1 (Intercept) 13.286 3.645
## group2 (Intercept) 1.299 1.140
## Residual 9.385 3.064
## Number of obs: 250, groups: group1, 5; group2, 25
##
## Dispersion estimate for gaussian family (sigma^2): 9.39
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.74849 1.65729 1.055 0.291
## x 0.02553 0.02459 1.038 0.299
------------------------------------------------------
今回のデータはnestedを想定して、下記のように生成した。今回は、上記のnestedよりももっとあり得る状況として、異なるgroup1について、異なる平均の正規分布からgroup2のランダム効果が生成されるようにしている。
------------------------------------------------------
b0 <- 2
b1 <- 0
s1 <- 3
s2 <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rnorm(g1*g2, mean = rep(rnorm(g1, mean = 0, sd = s1),
each = g2), sd = is2), each = n)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s2)
------------------------------------------------------
その時の、ランダム切片効果を記録してある。下記のとおりである。この時のgroup1、2のランダム効果の関係性を確認してみよう。
------------------------------------------------------
i_sds1 <- c(-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, 2.3, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -4.7, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4, -3.4)
i_sds2 <- c(0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, -0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.6, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, -1.9, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.8, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3, -1.3)
plotn(i_sds2 ~ i_sds1)#図39の描画
------------------------------------------------------
図39 横軸=group1、縦軸=group2のランダム効果
では、このタイプのnested designのデータの時、crossedあるいはnested designを想定した解析を行うとどうなるかシミュレーションしてみよう。
------------------------------------------------------
b0 <- 2
b1 <- 0
s1 <- 3
s2 <- 3
is1 <- 3
is2 <- 1
g1 <- 5
g2 <- 5
n <- 10
p1_c <- NULL
p2_c <- NULL
b0e_c <- NULL
b1e_c <- NULL
is_e1_c <- NULL
is_e2_c <- NULL
p1_n <- NULL
p2_n <- NULL
b0e_n <- NULL
b1e_n <- NULL
is_e1_n <- NULL
is_e2_n <- NULL
aic_c <- NULL
aic_n <- NULL
group1 <- rep(LETTERS[1:g1], each = n*g2)
group2 <- rep(rep(letters[1:g2], each = n), times = g1)
for (i in 1:1000){
x <- round(runif(n*g1*g2, -15, 15), digits = 1)
i_sds1 <- rep(rnorm(g1, mean = 0, sd = is1), each = n*g2)
i_sds2 <- rep(rnorm(g1*g2, mean = rep(rnorm(g1, mean = 0, sd = s1),
each = g2), sd = is2), each = n)
y <- rnorm(n = n*g1*g2, mean = (b0+i_sds1+i_sds2) + b1*x, sd = s2)
d <- data.frame(x = x, group1 = group1, group2 = group2, y = y)
fit1 <- glmmTMB(y ~ x + (1|group1) + (1|group2), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル
fit2 <- summary(fit1)
fit3 <- glmmTMB(y ~ x + (1|group1/group2), data = d, family = gaussian(link = "identity")) #一般化線形混合モデル
fit4 <- summary(fit3)
p1_c <- c(p1_c, coef(fit2)$cond[1,4])
p2_c <- c(p2_c, coef(fit2)$cond[2,4])
b0e_c <- c(b0e_c, coef(fit2)$cond[1,1])
b1e_c <- c(b1e_c, coef(fit2)$cond[2,1])
is_e1_c <- c(is_e1_c, as.numeric(fit2$varcor$cond$group1[1,1]))
is_e2_c <- c(is_e2_c, as.numeric(fit2$varcor$cond$group2[1,1]))
p1_n <- c(p1_n, coef(fit4)$cond[1,4])
p2_n <- c(p2_n, coef(fit4)$cond[2,4])
b0e_n <- c(b0e_n, coef(fit4)$cond[1,1])
b1e_n <- c(b1e_n, coef(fit4)$cond[2,1])
is_e1_n <- c(is_e1_n, as.numeric(fit4$varcor$cond$group1[1,1]))
is_e2_n <- c(is_e2_n, as.numeric(fit4$varcor$cond$`group2:group1`[1,1]))
aic_c <- c(aic_c, fit2$AICtab[1])
aic_n <- c(aic_n, fit4$AICtab[1])
}
------------------------------------------------------
では、この時の危険率、検出力などを確認してみよう。
------------------------------------------------------
#crossed design
sum(p1_c < 0.05, na.rm = T)
## [1] 326
sum(p2_c < 0.05, na.rm = T)
## [1] 49
sum(!is.na(p1_c))
## [1] 1000
sum(!is.na(p2_c))
## [1] 1000
#nested design
sum(p1_n < 0.05, na.rm = T)
## [1] 330
sum(p2_n < 0.05, na.rm = T)
## [1] 55
sum(!is.na(p1_n))
## [1] 1000
sum(!is.na(p2_n))
## [1] 1000
#crossed design
histn(p1_c[!is.na(p1_c)], xlab = "P value", ylab = "頻度")#図40の描画
histn(p2_c[!is.na(p2_c)], xlab = "P value", ylab = "頻度")#図41の描画
#nested design
histn(p1_n[!is.na(p1_n)], xlab = "P value", ylab = "頻度")#図42の描画
histn(p2_n[!is.na(p2_n)], xlab = "P value", ylab = "頻度")#図43の描画
------------------------------------------------------
図40 crossed designの切片の係数のP値のヒストグラム
図41 crossed designのxの係数のP値のヒストグラム
図42 nested designの切片の係数のP値のヒストグラム
図43 nested designのxの係数のP値のヒストグラム
危険率などに関して、どちらの解析でも大差がないことがわかる。次は係数の推定値を確認してみる。
------------------------------------------------------
#crossed design
mean(b0e_c[!is.na(p1_c)])
## [1] 2.063145
mean(b1e_c[!is.na(p2_c)])
## [1] 2.523602e-05
var(b0e_c[!is.na(p1_c)])
## [1] 3.791186
var(b1e_c[!is.na(p2_c)])
## [1] 0.0005421164
#nested design
mean(b0e_n[!is.na(p1_n)])
## [1] 2.063094
mean(b1e_n[!is.na(p2_n)])
## [1] -4.706159e-05
var(b0e_n[!is.na(p1_n)])
## [1] 3.791442
var(b1e_n[!is.na(p2_n)])
## [1] 0.0005223886
#crossed design
histn(b0e_c[!is.na(p1_c)], xlab = "b0", ylab = "頻度")#図44の描画
histn(b1e_c[!is.na(p2_c)], xlab = "b1", ylab = "頻度")#図45の描画
#nested design
histn(b0e_n[!is.na(p1_n)], xlab = "b0", ylab = "頻度")#図46の描画
histn(b1e_n[!is.na(p2_n)], xlab = "b1", ylab = "頻度")#図47の描画
------------------------------------------------------
図44 crossed designの切片の係数のヒストグラム
図45 crossed designのxの係数のヒストグラム
図46 nested designの切片の係数のヒストグラム
図47 nested designのxの係数のヒストグラム
係数推定値も大差はない。ランダム効果の分散も確認してみる。
------------------------------------------------------
#crossed design
mean(is_e1_c[!is.na(p1_c)])
## [1] 14.5204
mean(is_e2_c[!is.na(p2_c)])
## [1] 0.2012162
#nested design
mean(is_e1_n[!is.na(p1_n)])
## [1] 14.29
mean(is_e2_n[!is.na(p2_n)])
## [1] 1.023415
#crossed design
histn(is_e1_c[!is.na(p1_c)], ylab = "頻度")#図48の描画
histn(is_e2_c[!is.na(p2_c)], ylab = "頻度")#図49の描画
#nested design
histn(is_e1_n[!is.na(p1_n)], ylab = "頻度")#図50の描画
histn(is_e2_n[!is.na(p2_n)], ylab = "頻度")#図51の描画
------------------------------------------------------
図48 crossed designの切片のランダム効果の分散のヒストグラム
図49 crossed designのxの係数のランダム効果の分散のヒストグラム
図50 nested designの切片のランダム効果の分散のヒストグラム
図51 nested designのxの係数のランダム効果の分散のヒストグラム
さて、一つ前のnestedな解析の時と同様に、crossedな解析ではgroup2のランダム効果の分散が過少推定になる。もう一点、興味深いのが、group2ではなくgroup1の分散が2つの解析で今までの解析よりも大きくなっていることだ。私はgroup2の分散が大きくなると思っていた。
なぜこんなことを言うかというと、データの生成過程を確認してみると、group1は今までと同じように生成し、group2の生成の仕方を変更している。一つ前のnestedなデータと比較して、group2のランダム効果は、今回の場合、平均が「平均0、分散3^2の正規分布に従うとき」の分散1^2の正規分布から生成されていることになる。正規-正規混合分布の性質から、この時の混合分布は平均0、分散9 + 1 = 10の正規分布に従う。ゆえに、group2の分散が10程度になると期待した。
しかし、実際の推定では、平均の従う正規分布の分散は影響がなかった。つまり、group2のランダム効果の分散の推定はデータ全体ではなく、各group1で行われているのだと推定できる。では、group2の平均のばらつきの影響はどこに行ったかというと、group1に押し付けられたといえよう。ここで、天下り的に考えると、group1が平均が「分散3^2の正規分布に従うとき」の分散3^2の分布の混合分布は、その分散が9 + 9 = 18である。これの4/5倍が大体14である。不偏分散と母分散の関係から、サンプルサイズが5なら、シミュレーションで得た分散=不偏分散=4/5×混合分布の分散=4/5×母分散である。おそらく、group1の分散の推定は、group2の平均が従う正規分布の分散とgroup1の分散の和となっている。私の手元でいくつか計算してみると、おおむねこの予想通りとなっている。
最後、AICを比較してみよう。
------------------------------------------------------
#crossed design
mean(aic_c[!is.na(p1_c)], na.rm = T)
## [1] 1306.751
#nested design
mean(aic_n[!is.na(p1_n)], na.rm = T)
## [1] 1301.178
#crossed design
histn(aic_c[!is.na(p1_c)], ylab = "頻度")#図52の描画
#nested design
histn(aic_n[!is.na(p1_n)], ylab = "頻度")#図53の描画
------------------------------------------------------
図52 crossed designのAICのヒストグラム
図53 nested designのAICのヒストグラム
AICを見ると、nested designのほうが予測性が良いことが変わる。その意味ではより適切なモデルはデータ生成過程通り、nested designである。
以上、いくつかのパターンから考えられることは、「nestedな関係があるなら考慮したほうが予測性は高まるが、推定値には大きな影響がないだろう」ということになる。そういう意味では、nest関係が複雑なデータにおいて、過剰に神経質になる必要はないかもしれない。