目次

・階層的なランダム効果

・Crossed design vs. nested design

・具体的なデータの解析：Crossed design

・具体的なデータの解析：nested design

　　・group2が共通の平均を持つ場合

　　・group2が異なる平均を持つ場合

階層的なランダム効果

　様々なデータを解析していて、頭を悩ます問題の一つがランダム効果に階層性が認められるときである。例えば、ランダム効果の一つとして、栽培実験を行ったときにいくつかの区画（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の描画

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

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

}

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

さて、もしこれを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の描画

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

データの形としては、このようなランダム効果の関係となっていれば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の描画

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

危険率などに関して、どちらの解析でも大差がないことがわかる。次は係数の推定値を確認してみる。

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

#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の描画

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

係数推定値も大差はない。ランダム効果の分散も確認してみる。

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

#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の描画

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

こちらも推定に大きな差はなさそうだ。最後、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の描画

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

AICを見ると、crossed designのほうが予測性が良いことが変わる。その意味ではより適切なモデルはデータ生成過程通り、crossed designである。

具体的なデータの解析：nested design

・group2が共通の平均を持つ場合

　次に、以下のようなデータが得られているとき、どうやって解析していくかを考える。今回のデータは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の描画

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

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

}

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

さて、もしこれを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の描画

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

データの形としては、このようなランダム効果の関係となっていれば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の描画

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

危険率などに関して、どちらの解析でも大差がないことがわかる。次は係数の推定値を確認してみる。

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

#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の描画

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

係数推定値も大差はない。ランダム効果の分散も確認してみる。

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

#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の描画

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

ここで初めて、差が生じた。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の描画

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

AICを見ると、nested designのほうが予測性が良いことが変わる。その意味ではより適切なモデルはデータ生成過程通り、nested designである。

・group2が異なる平均を持つ場合

　最後に、もう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)