1 環境統計学ぷらす 第 5 回 一般 ( 化 ) 線形混合モデル 高木俊 2013/11/21

第

₅

回

一般（化）線形混合モデル

高木 俊

[email protected]

予定

•

第

1

回：

R

の基礎と仮説検定

•

第

2

回： 分散分析と回帰

•

第

3

回： 一般線形モデル・交互作用

•

第

4.1

回： 一般化線形モデル

•

第

4.2

回： モデル選択（

11/29?

）

•

第

5

回： 一般化線形混合モデル

•

第

6

回： 多変量解析（

12/5?

）

今日やること

•

統計編

•

変量効果と混合モデル

•

実験デザイン

一般化線形

混合モデル

データの独立性

•

今まで扱ってきたモデルはデータの独立性を仮定

例）餌種が昆虫の成長量に与える影響を見たい

野外で虫を

₂₀

匹採ってきて、

₁₀

匹に餌

_A

を、

₁₀

匹に餌

_B

を与えて成長量を見る

→餌種と成長量の関係に対する

20

個の独立なデータとみなせる

N=20

_×

10

×

10

0 .0 0 .5 1 .0 1 .5 2 .0 餌A 餌B

偽の繰り返し

（

_{pseudoreplication}

）

•

独立性が満たされない状況

野外で虫を

₄

匹採ってきて、それぞれ

₅

匹の子供が産まれた（子の合計

₂₀

）

親

_A

・

_B

から産まれた

₁₀

匹に餌

_A

を、親

_C

・

_D

から産まれた

₁₀

匹に餌

_B

を与えた

→同親由来の子どもは成長も似てきそう

20

個の独立なデータとは言えない！（真の繰り返しは

4

）

×

10

×

₁₀

×

₄

A B C D 0 .0 0 .5 1 .0 1 .5 2 .0

親

餌A 餌B

どのように解析するか

・成長＝処理＋誤差

_௜௞

= 

_௜

+ 

_௜௞



୲୰୲

=

MS

_trt

MS

_resid

_{処理内の分散（残差）}

処理間の分散

_௜௝௞

= 

_௜

+



_௝(௜)

+ 

_௜௝௞

_

୲୰୲

=

MS

_trt

MS

_parent

処理間の分散

処理

i

の中の

j

番目

の親の効果

処理

_i

の効果

・成長＝処理＋親間の差＋誤差

親間（処理内）の分散

どのように解析するか

_௜௝௞

=



_௜

+



_௝

_(௜)

+



_௜௝௞

_

୲୰୲

=

MS

_trt

MS

_parent

親間（処理内）の分散

処理間の分散

処理

_i

の中の

_j

番目の親の効果

処理

_i

の効果

・成長＝処理＋親間の差＋誤差

Total

ばらつき

＝

処理間ばらつき

＋

（処理内）親間ばらつき

＋

（個体間）誤差

処理間ばらつき

が

親間ばらつき

に比べて十分大きければ処理の効果がある

変量効果と混合モデル

・処理の効果はまさに

_i

番目の処理がどのように影響するかに注目

→

固定効果（

fixed effect

）

・親

_j

の効果は個々の効果そのものではなく、その

バラつき

に注目

→

変量効果（ランダム効果・

random effect

）

両者を含むモデルが混合モデル（

Mixed effect model

）

一般線形混合モデル

一般化線形混合モデル

練習問題：

_aov

による解析

1.

親の効果を無視して解析

2.

親の効果を変量効果として解析

3.

親の効果を固定効果として解析するのと何が違う？

4.

親ごとに平均値を計算して、最初からデータ数

4

で解析する

のと何が違う？

data5.1<- read.csv("data5.1.csv",T)

data5.1$parent<- as.factor(data5.1$parent)

#

カテゴリカル変数として認識させる

summary(aov(growth~trt,data5.1))

summary(aov(growth~trt+Error(parent),data5.1))

summary(aov(growth~trt+parent,data5.1))

par.growth<-tapply(data5.1$growth,data5.1$parent,mean)

par.trt<- factor(c("c","c","t","t"))

summary(aov(par.growth~par.trt))

出力結果

Error: parent

Df Sum Sq Mean Sq F value Pr(>F)

trt 1 1.7553 1.7553 8.218 0.103

Residuals 2 0.4272 0.2136

Error: Within

Df Sum Sq Mean Sq F value Pr(>F)

Residuals 16 0.124 0.007749

> summary(aov(growth~trt+Error(parent),data5.1))

親を誤差とした

分散分析表

このResidualは（処理内）親間の誤差 こっちResidualは（親内）個体間の誤差

個体を誤差とし

た分散分析表

親間の誤差と比較すると処理の効果 は有意ではない（P=0.1）

ここまでのまとめ

•

データが独立ではなく、

同一親・同一地点

といったデータの

ばらつきに影響する要因が入っている場合、それらを

変量効

果

として考慮する必要がある

•

変量効果を無視すると個々のデータは

偽の繰り返し

となり、

差がないものが有意になる（

_{Type I Error}

）、差があるものが

検出できない（

_{Type II Error}

）ことがある

•

変量効果を入れるかどうか・どのように入れるかは、データの

取り方（実験デザイン・調査デザイン）による

実験デザイン

Completely randomized

Nested

Randomized complete block

（

unreplicated

）

Completely

randomized factorial

Split-plot

実験・調査デザインで

解析も変わってくる

Nested design

•

要因

A

の中に要因

B

が入れ子、繰り返しあり

ܣ

_ଵ

ܤ

_ଵ(ଵ)

ܤ

_ଶ(ଵ)

ܤ

_ଷ(ଵ)

ܣ

_ଶ

ܤ

_ଵ(ଶ)

ܤ

_ଶ(ଶ)

ܤ

_ଷ(ଶ)

ݕ

_{ଵ(ଵଵ)ଵ}

ݕ

_{ଵ(ଵଵ)ଶ}

ݕ

_{ଵ(ଷଵ)ଵ}

ݕ

_{ଵ(ଷଵ)ଶ}

…

ݕ

_{ଶ(ଵଶ)ଵ}

ݕ

_{ଶ(ଵଶ)ଶ}

ݕ

_{ଶ(ଷଶ)ଵ}

ݕ

_{ଶ(ଷଶ)ଶ}

…

_௜௝௞

= 

_௜

+



_௝(௜)

+ 

_௜௝௞

A:fixed B:random



_஺

=

MS

A

MS

_୆(୅)

地点間の分散

ヒシ

vs

開放間の分散

ヒシの効果 地点間誤差 地点内場所間誤差

例：

ヒシ帯

₂

地点・開放水面

₂

地点があり、

それぞれの地点で

2

箇所ずつ測定

Trapa

ୗ୧୲ୣ(୘୰ୟ୮ୟ)

(

観測数

)

Site1

Site2

Site3

Site4

Open

2

2

Trapa

2

2

Nested design

地点間の分散

ヒシ

vs

開放水面間の分散

例：

ヒシ帯

₂

地点・開放水面

₂

地点があり、

それぞれの地点で

2

箇所ずつ測定

Open

Trapa

Site2

Site3

Site4

Randomized block design

(

乱塊法

)

•

要因

A

と要因

B

は直交、繰り返しなし

_௜௝

= 

_௜

+



_௝

+ 

_௜௝

A:fixed B:random



_஺

=

MS

A

MS

_resid

ܣ

_ଵ

ܤ

_ଵ

ܤ

_ଶ

ܤ

_ଷ

ݕ

_ଵଵ

ݕ

_ଵଶ

ݕ

_ଵଷ

ܣ

_ଶ

ܤ

_ଵ

ܤ

_ଶ

ܤ

_ଷ

ݕ

_ଶଵ

ݕ

_ଶଶ

ݕ

_ଶଷ ヒシの効果 ブロック間誤差 いずれでも説明されない残差

例：

4

ブロックがあり、それぞれのブロックで

ヒシ刈り取り・刈り残し

1

地点ずつ測定

処理・ブロックで説明さ

れない分散（残差）

刈り取り

vs

刈り残し間の分散

(

観測数

)

Block1

Block2

Block3

Block4

刈取（

O)

1

1

1

1

刈残（

T)

1

1

1

1

Randomized block design

Trapa

୰ୣୱ୧ୢ

説明できない誤差

刈り残し

vs

刈り取り間の分散

例：

4

ブロックがあり、それぞれのブロックで

刈り取り・刈り残し

1

地点ずつ測定

Blk2

Blk3

Blk4

Blk1

T

O

O

T

O

T

O

T

Factorial design

•

要因

A

と要因

B

は直交、各組み合わせに繰り返しあり

ܣ

_ଵ

ܤ

_ଵ

ܤ

_ଶ

ܤ

_ଷ

ܣ

_ଶ

ܤ

_ଵ

ܤ

_ଶ

ܤ

_ଷ

ݕ

_ଵଵଵ

ݕ

_ଵଵଶ

ݕ

_ଵଷଵ

ݕ

_ଵଷଶ

…

ݕ

_ଶଵଵ

ݕ

_ଶଵଶ

ݕ

_ଶଷଵ

ݕ

_ଶଷଶ

…

_௜௝௞

= 

_௜

+



_௝

+



_௜௝

+ 

_௜௝௞

A:fixed B:random



_஺

=

MS

A

MS

_AB

例：

2

ブロックがあり、それぞれのブロックで

刈り取り・刈り残し

2

地点ずつ測定

ヒシの効果 ブロック間誤差 いずれでも説明されない残差 ヒシの効果のブロック間誤差

刈り取りの効果のブロック間

でのばらつき

刈り取り

vs

刈り残し間の分散

Factorial design

例：

2

ブロックがあり、それぞれのブロックで

刈り取り・刈り残し

2

地点ずつ測定

୘୰ୟ୮ୟ

୘୰ୟ୮ୟ:୆୪୭ୡ୩

(

観測数

)

Block1

Block2

刈取（

O)

2

2

刈残（

T)

2

2

Block2

Block1

T

O

O

T

刈り残し

_vs

刈り取り間の分散

刈り取り効果のブロック間での分散

実験デザインによる

_ANOVA

の違い

Completely

randomized

（独立）

Nested

（入れ子）

Randomized

block

（乱塊）

Factorial

（直交）

_௜௝௞

= 

_௜

+ 

_௜௝௞

_௜௝௞

= 

_௜

+



_௝

_(௜)

+ 

_௜௝௞

_௜௝௞

= 

_௜

+



_௝

+



_௜௝

+ 

_௜௝௞

_௜௝

= 

_௜

+



_௝

+ 

_௜௝



_஺

=

MS

A

MS

_୆(୅)



_஺

=

MS

A

MS

_resid



_஺

=

MS

A

MS

_resid



_஺

=

MS

A

MS

_AB

処理の効果/ 処理で説明できない残差 処理の効果/ 処理とブロックの効果を除いた残差 （ブロック共通の）処理の効果/ブロックごとの処理の効果のばらつき 処理 残差 処理の効果/ 処理内のブロック間のばらつき 処理内ブロック ブロック 処理：ブロック交互作用

練習問題

₂

：シカの効果と調査デザイン

A)

シカのいる

3

地域のシカ柵内と柵外で

2

個体の樹高測定

B)

シカのいる

6

地域のシカ柵内と柵外で

1

個体の樹高測定

C)

シカのいる

3

地域、いない

3

地域で

2

個体ずつ樹高測定

D)

シカのいる

6

地域、いない

6

地域で

1

個体ずつ樹高測定

×

₃

地域

×

3

地域

×

3

地域

×

₆

地域

A

）

C

）

B

）

×

₆

地域

_×

₆

_地域

D

）

シカ

₆

・コントロール

₆

のデザインでも全部違う！

deer site

1 1

1 1

0 1

0 1

1 2

1 2

0 2

0 2

1 3

1 3

0 3

0 3

×

3

A

）

deer site

1 1

0 1

1 2

0 2

1 3

0 3

1 4

0 4

1 5

0 5

1 6

0 6

×

₆

B

）

deer site

1 1

1 1

1 2

1 2

1 3

1 3

0 4

0 4

0 5

0 5

0 6

0 6

×

₃

×

₃

C

）

deer site

1 1

1 2

1 3

1 4

1 5

1 6

0 7

0 8

0 9

0 10

0 11

0 12

×

₆

×

₆

D

）

R

解析前の注意点

•

シカの在不在、地点

id

は数字で入っているので、これをカテゴ

リカル変数として認識させる必要がある

連続変数のまま解析すると異なる解析を行うことに（エラーは出ない）

#読み込み

data5.2a<- read.csv("data5.2a.csv",T)

data5.2b<- read.csv("data5.2b.csv",T)

data5.2c<- read.csv("data5.2c.csv",T)

data5.2d<- read.csv("data5.2d.csv",T)

#deerをカテゴリー変数として認識させる

data5.2a$deer<- as.factor(data5.2a$deer)

data5.2b$deer<- as.factor(data5.2b$deer)

data5.2c$deer<- as.factor(data5.2c$deer)

data5.2d$deer<- as.factor(data5.2d$deer)

#siteをカテゴリー変数として認識させる

data5.2a$site<- as.factor(data5.2a$site)

data5.2b$site<- as.factor(data5.2b$site)

data5.2c$site<- as.factor(data5.2c$site)

data5.2d$site<- as.factor(data5.2d$site)

Completely randomized factorial design

×

₃

summary(aov(height~deer+Error(site/deer),dat5.2a))

 

_௜௝௞

= 

_௜

+ 

_௝

+  · 

_௜௝

+ 

_௜௝௞

固定効果

変量効果

残差

site

で説明されるばらつき

deer

で説明されるばらつき

deer:site

で説明されるばらつき

残り（個体差）

Error: site

Df Sum Sq Mean Sq F value Pr(>F) Residuals 2 0.8606 0.4303 Error: site:deer

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.13152 0.13152 48.29 0.0201 * Residuals 2 0.00545 0.00272

---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Error: Within

Df Sum Sq Mean Sq F value Pr(>F) Residuals 6 0.0569 0.009483

変量効果

（交互作用）

Randomized complete block design

 

_௜௝

= 

_௜

+ 

_௝

+ 

_௜௝

固定効果

変量効果

残差

site

で説明されるばらつき

deer

で説明されるばらつき

残り（個体差かつ効果の地域差）

×

₆

Error: site

Df Sum Sq Mean Sq F value Pr(>F) Residuals 5 0.975 0.195 Error: Within

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.13577 0.13577 12.03 0.0179 * Residuals 5 0.05641 0.01128

---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nested design

 

_௜௝

= 

_௜

+ 

_௝

+ 

_௜௝

固定効果

変量効果

残差

site

で説明されるばらつき

（かつ効果の地域差）

deer

で説明されるばらつき

残り（個体差）

Error: site

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.029 0.02902 0.123 0.744 Residuals 4 0.944 0.23601 Error: Within

Df Sum Sq Mean Sq F value Pr(>F) Residuals 6 0.0569 0.009483

summary(aov(height~deer+Error(site),dat5.2c))

×

₃

×

₃

ݏ݅ݐ݁

_௝(௜)

とも表記（RCBと区別）

Completely randomized design

 

_௜௝

= 

_௜

+ 

_௜௝

固定効果

残差

deer

で説明されるばらつき

残り（個体差かつ地域差

かつ効果の地域差）

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.1951 0.1951 1.698 0.222 Residuals 10 1.1492 0.1149

summary(aov(height~deer,dat5.2d))

×

₆

×

₆

> summary(aov(height~deer*site,data5.2a))

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.1315 0.1315 13.868 0.009806 ** site 2 0.8606 0.4303 45.373 0.000239 *** deer:site 2 0.0054 0.0027 0.287 0.760136 Residuals 6 0.0569 0.0095 > summary(aov(height~deer*site,dat5.2b)) Df Sum Sq Mean Sq deer 1 0.1358 0.13577 site 5 0.9750 0.19501 deer:site 5 0.0564 0.01128 > summary(aov(height~deer*site,dat5.2c))

Df Sum Sq Mean Sq F value Pr(>F) deer 1 0.0290 0.02902 3.06 0.130838 site 4 0.9440 0.23601 24.89 0.000703 *** Residuals 6 0.0569 0.00948 > summary(aov(height~deer*site,dat5.2d)) Df Sum Sq Mean Sq deer 1 0.1951 0.1951 site 10 1.1492 0.1149

F

1,2

＝

0.1315/0.0027

P

＝

0.020

1-pf(0.1315/0.0027,1,2)

左の表からも正しい

F,P

を計算可能

1-pf(0.13577/0.01128,1,5)

F

1,5

＝

0.13577/0.01128

P

＝

0.018

1-pf(0.02902/0.23601,1,4)

F

1,4

＝

0.02902/0.23601

P

＝

0.744

1-pf(0.1951/0.1149,1,10)

F

1,10

＝

0.1951/0.1149

P

＝

0.222

×

₃

地域の効果

シカの効果の地域差

個体差

×

6

シカの効果

地域の効果

個体差＋シカの効果の地域差

×

₃

×

3

シカの効果

地域の効果＋シカの効果の地域差

個体差

×

₆

×

₆

シカの効果

地域の効果＋シカの効果の地域差＋個体差

Randomized

complete

block

乱塊

randomized

factorial

直交

Nested

入れ子

Completely

randomized

独立

一般化線形混合モデル（

_GLMM

）

• glmmML() in library(glmmML)

glm

と似た感覚で使える。

binomial

と

poisson

のみ、ランダム効果はひとつ。

• glmer() in library(lme4)

様々な分布、複数のランダム効果も扱える。

• glmmadmb() in library(glmmADMB)

更に様々な分布で

GLM

・

GLMM

の解析可能。

計算時間やや長い。

•

今までのは正規分布を仮定した、一般線形混合モデル（

LMM

と書くことも）

•

一般線形混合モデル（

LMM

）で

正規分布以外

を仮定する

一般化線形

混合モデル（

_GLMM

）

も存在

R

において

GLMM

の解析が可能な関数

#LMM

は

lmer

で解析

#R

のバージョンによっ

ては

website

からのイン

ストール



_௜௝

= 

_௜

+ 

_௝

log λ

_௜௝

= 

_௜

+ 

_௝

・一般線形混合モデル

・一般化線形混合モデル

_௜௝௞

~(λ

_௜௝

)



_௝

~(0, 

_௦௜௧௘

2

₎

_௜௝௞

~(

_௜௝

, 

_ఌ

2

₎



_௝

~(0, 

_௦௜௧௘

2

)

処理の効果

場所間ばらつき

ただ、最尤法を使っての推定なので、計算は複雑

aov(y~trt+Error(site))

glmmML(y~trt,cluster=site,family=poisson)

glmer(y~trt+(1|site),family=poisson)

glmmadmb(y~trt+(1|site),family=“poisson”)

例題：シカがいるとアオキは減るか

•

シカの分布（

3

地域）、非分布（

2

地域）で各地域

5

株に関して周

囲の♀株の数を測定（

Nested design

）

非負整数値

family=poisson

atago fudago godai hiratuka kimikame

0 5 1 0 1 5

site

fe

m

a

le

glmer.model5.3<- glmer(female~deer+(1|site),data5.3,family=poisson)

summary(glmer.model5.3)

Generalized linear mixed model fit by maximum likelihood ['glmerMod'] Family: poisson ( log )

Formula: female ~ deer + (1 | site) Data: data5.3

AIC BIC logLik deviance 125.0005 128.6572 -59.5003 119.0005 Random effects:

Groups Name Variance Std.Dev. site (Intercept) 0.1544 0.393 Number of obs: 25, groups: site, 5 Fixed effects:

Estimate Std. Error z value Pr(>|z|) (Intercept) 1.0085 0.2760 3.654 0.000258 *** deernd 1.2864 0.4043 3.182 0.001463 **

---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects:

(Intr) deernd -0.683

log λ

_௜௝

= ߙ

_௜

+ ߛ

_௝

ߛ

_௝

~ܰ݋ݎ݈݉ܽ(0, ߪ

2

₎

係数の推定値に対するWald検定によるP値 変数の効果をみるなら尤度比検定・Parametric Bootstrapのほうが良い

シカ在期待密度：

λ

d

= exp(1.0085)

シカ不在期待密度：

λ

nd

= exp(1.0085＋1.2864)

ߙ

_௜

•

カイ二乗近似を用いた尤度比検定

•

Parametric bootstrap

法による近似によらない検定

glmer.model5.3.0<- glmer(female~1+(1|site),data5.3,family=poisson)

anova(glmer.model5.3.0,glmer.model5.3)

Data: data5.3 Models:

glmer.model5.3.0: female ~ 1 + (1 | site) glmer.model5.3: female ~ deer + (1 | site)

Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) glmer.model5.3.0 2 128.35 130.79 -62.177 124.36 glmer.model5.3 3 125.00 128.66 -59.500 119.00 5.3543 1 0.02067 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ←見たい効果入り（deer） ←見たい効果抜き（null） 2×対数尤度の差（尤度比）

「見たい効果を抜いたモデルの推定値の元でデータを発生させ、

実データと同様の解析を行い、尤度比を計算させる」を繰り返し行

い（

₁₀₀₀

回とか）、近似によらない尤度比の分布を得る方法。計

算に時間がかかる。（詳しい方法はスクリプト参照）

シカの効果ありそう

実は

_PB

法では

P=0.08

～

0.09

になり、

有意ではない

まとめ

•

個々の効果の強さに興味が有る→固定効果

•

ばらつきを考慮したいだけ→変量効果

•

地域など空間的なものだけでなく、個体差、時間なども同様に

モデリング可（同一個体の反復測定など）

•

lmer()

は正規分布も使えるが、近似による尤度比検定の

P

値

は

_aov()

と異なる場合がある。素直に

_aov()

で示したほうが分

散分析表が得られ、結果がわかりやすい

•

ややこしい階層構造を持つモデルは、無理に

GLMM

でやるより、階層ベイ

ズの方が推定しやすいかも

参考文献（

_GLM

・

_GLMM

）

•

Statistics: An introduction using R

•

日本語訳あり（共立出版）

•

aov()

レベルまでなら

•

ランダム要因の解説とかわかりやすい

•

Extending the linear model with R

：

Generalized Linear, Mixed

Effects and Nonparametric Regression Models

•

glm()

、

lmer()

あたりの解説詳しい

•

初学者には難しいかも