階層的なランダム効果が

仮定されるときの考え方

階層的なランダム効果

 様々なデータを解析していて、頭を悩ます問題の一つがランダム効果に階層性が認められるときである。例えば、ランダム効果の一つとして、栽培実験を行ったときにいくつかの区画(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関係が複雑なデータにおいて、過剰に神経質になる必要はないかもしれない。