36
予定
•第1回： Rの基礎と仮説検定•第2回： 分散分析と回帰•第3回： 一般線形モデル・交互作用•第4.1回： 一般化線形モデル•第4.2回： モデル選択（11/29?）•第5回： 一般化線形混合モデル•第6回： 多変量解析（12/5?）
2
今日やること
• 統計編
• 変量効果と混合モデル
• 実験デザイン
• 一般化線形混合モデル
3
一般化線形混合モデル
4
データの独立性
•今まで扱ってきたモデルはデータの独立性を仮定例）餌種が昆虫の成長量に与える影響を見たい
野外で虫を20匹採ってきて、10匹に餌Aを、10匹に餌Bを与えて成長量を見る
→餌種と成長量の関係に対する20個の独立なデータとみなせる
N=20×10
×10
0.0
0.5
1.0
1.5
2.0
餌A 餌B
5
偽の繰り返し（pseudoreplication）
• 独立性が満たされない状況
野外で虫を4匹採ってきて、それぞれ5匹の子供が産まれた（子の合計20）親A・Bから産まれた10匹に餌Aを、親C・Dから産まれた10匹に餌Bを与えた
→同親由来の子どもは成長も似てきそう
20個の独立なデータとは言えない！（真の繰り返しは4）
×10
×10×4
A B C D
0.0
0.5
1.0
1.5
2.0
親
餌A
餌B
6
どのように解析するか
・成長＝処理＋誤差
��� = �� + ������� =
MStrt
MSresid 処理内の分散（残差）
処理間の分散
���� = �� + ��(�) + ���� ���� =MStrt
MSparent
処理間の分散
処理iの中のj番目の親の効果
処理iの効果
・成長＝処理＋親間の差＋誤差
親間（処理内）の分散
7
どのように解析するか
���� = �� + ��(�) + ���� ���� =MStrt
MSparent 親間（処理内）の分散
処理間の分散
処理iの中のj番目の親の効果処理iの効果
・成長＝処理＋親間の差＋誤差
Totalばらつき ＝ 処理間ばらつき （処理内）親間ばらつき （個体間）誤差＋ ＋
処理間ばらつきが親間ばらつきに比べて十分大きければ処理の効果がある
8
変量効果と混合モデル
9
・処理の効果はまさに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))
10
出力結果
11
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）ことがある
• 変量効果を入れるかどうか・どのように入れるかは、データの取り方（実験デザイン・調査デザイン）による
12
実験デザインCompletely randomized Nested Randomized complete block
（unreplicated）
Completely randomized factorial
Split-plot
実験・調査デザインで
解析も変わってくる
13
Nested design• 要因Aの中に要因Bが入れ子、繰り返しあり
��
��(�)
��(�)
��(�)
��
��(�)
��(�)
��(�)
��(��)�
��(��)�
��(��)�
��(��)�
…
��(��)�
��(��)�
��(��)�
��(��)�…
���� = �� + ��(�) + ����
A:fixed B:random
�� =MSA
MS�(�) 地点間の分散
ヒシvs開放間の分散
ヒシの効果 地点間誤差 地点内場所間誤差
例：
ヒシ帯2地点・開放水面2地点があり、それぞれの地点で2箇所ずつ測定
14
Trapa
����(�����)
(観測数) Site1 Site2 Site3 Site4
Open 2 2
Trapa 2 2
Nested design
地点間の分散
ヒシvs開放水面間の分散
例：
ヒシ帯2地点・開放水面2地点があり、それぞれの地点で2箇所ずつ測定
Open Trapa
Site2 Site3 Site4Site1
15
Randomized block design(乱塊法)
• 要因Aと要因Bは直交、繰り返しなし
��� = �� + �� + ���
A:fixed B:random
�� =MSA
MSresid
��
��
��
��
���
���
���
��
��
��
��
���
���
���
ヒシの効果 ブロック間誤差 いずれでも説明されない残差
例：
4ブロックがあり、それぞれのブロックでヒシ刈り取り・刈り残し1地点ずつ測定
処理・ブロックで説明されない分散（残差）
刈り取りvs刈り残し間の分散
16
(観測数) Block1 Block2 Block3 Block4
刈取（O) 1 1 1 1
刈残（T) 1 1 1 1
Randomized block design
Trapa
���
説明できない誤差
刈り残しvs刈り取り間の分散
例：
4ブロックがあり、それぞれのブロックで刈り取り・刈り残し1地点ずつ測定
Blk2 Blk3 Blk4Blk1
TO TO TO TO
17
Factorial design• 要因Aと要因Bは直交、各組み合わせに繰り返しあり
��
��
��
��
��
��
��
��
����
����
����
����
…
����
����
����
����…
���� = �� + �� + ���� + ����
A:fixed B:random
�� =MSA
MSAB
例：
2ブロックがあり、それぞれのブロックで刈り取り・刈り残し2地点ずつ測定
ヒシの効果 ブロック間誤差
いずれでも説明されない残差
ヒシの効果のブロック間誤差
刈り取りの効果のブロック間でのばらつき
刈り取りvs刈り残し間の分散
18
Factorial design
例：
2ブロックがあり、それぞれのブロックで刈り取り・刈り残し2地点ずつ測定
�����
�����:�� ��
(観測数) Block1 Block2
刈取（O) 2 2
刈残（T) 2 2
Block2Block1
TO TO
刈り残しvs刈り取り間の分散
刈り取り効果のブロック間での分散
19
実験デザインによるANOVAの違い
20
Completely randomized（独立）
Nested（入れ子）
Randomized block（乱塊）
Factorial（直交）
���� = �� + ����
���� = �� + ��(�) + ����
���� = �� + �� + ���� + ����
��� = �� + �� + ���
�� =MSA
MS�(�)
�� =MSA
MSresid
�� =MSA
MSresid
�� =MSA
MSAB
処理の効果 / 処理で説明できない残差
処理の効果 / 処理とブロックの効果を除いた残差
（ブロック共通の）処理の効果 / ブロックごとの処理の効果のばらつき
処理 残差
処理の効果 / 処理内のブロック間のばらつき処理内ブロック
ブロック
処理：ブロック交互作用
練習問題2：シカの効果と調査デザインA) シカのいる3地域のシカ柵内と柵外で2個体の樹高測定B) シカのいる6地域のシカ柵内と柵外で1個体の樹高測定C) シカのいる3地域、いない3地域で2個体ずつ樹高測定D) シカのいる6地域、いない6地域で1個体ずつ樹高測定
×3地域 ×3地域 ×3地域
×6地域
A）C）
B）
×6地域 ×6地域
D）
シカ6・コントロール6のデザインでも全部違う！
21
deer site 1 11 10 10 11 21 20 20 21 31 30 30 3
×3
A）
deer site1 10 11 20 21 30 31 40 41 50 51 60 6
×6
B）
deer site1 11 11 21 21 31 30 40 40 50 50 60 6
×3
×3
C）
deer site1 11 21 31 41 51 60 70 80 90 100 110 12
×6
×6
D）
R上でのデータ記述法の違い22
解析前の注意点
• シカの在不在、地点idは数字で入っているので、これをカテゴリカル変数として認識させる必要がある
23
連続変数のまま解析すると異なる解析を行うことに（エラーは出ない）
#読み込み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
×3
A）直交
summary(aov(height~deer+Error(site/deer),dat5.2a))
������� = �� + ��� + � · ���� + ����
固定効果 変量効果 残差
siteで説明されるばらつき
deerで説明されるばらつき
deer:siteで説明されるばらつき
残り（個体差）
Error: siteDf Sum Sq Mean Sq F value Pr(>F)
Residuals 2 0.8606 0.4303
Error: site:deerDf 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: WithinDf Sum Sq Mean Sq F value Pr(>F)
Residuals 6 0.0569 0.009483
変量効果
（交互作用）
24
Randomized complete block designB）乱塊
������ = �� + ��� + ���
固定効果 変量効果 残差
siteで説明されるばらつき
deerで説明されるばらつき
残り（個体差かつ効果の地域差）
×6
Error: siteDf Sum Sq Mean Sq F value Pr(>F)
Residuals 5 0.975 0.195
Error: WithinDf 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
summary(aov(height~deer+Error(site),dat5.2b))
25
Nested designC）入れ子
������ = �� + ��� + ���
固定効果 変量効果 残差
siteで説明されるばらつき（かつ効果の地域差）
deerで説明されるばらつき
残り（個体差）
Error: siteDf Sum Sq Mean Sq F value Pr(>F)
deer 1 0.029 0.02902 0.123 0.744Residuals 4 0.944 0.23601
Error: WithinDf Sum Sq Mean Sq F value Pr(>F)
Residuals 6 0.0569 0.009483
summary(aov(height~deer+Error(site),dat5.2c))
×3
×3
�����(�)
とも表記（RCBと区別）
26
Completely randomized designD）独立
������ = �� + ���
固定効果 残差
deerで説明されるばらつき残り（個体差かつ地域差かつ効果の地域差）
Df Sum Sq Mean Sq F value Pr(>F)deer 1 0.1951 0.1951 1.698 0.222Residuals 10 1.1492 0.1149
summary(aov(height~deer,dat5.2d))
×6
×6
27
何も考えずaov(height~deer*site)で解析すると・・・
> 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.13577site 5 0.9750 0.19501deer: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.1951site 10 1.1492 0.1149
F1,2＝0.1315/0.0027P＝0.020
1-pf(0.1315/0.0027,1,2)左の表からも正しいF,Pを計算可能
1-pf(0.13577/0.01128,1,5)
F1,5＝0.13577/0.01128P＝0.018
1-pf(0.02902/0.23601,1,4)
F1,4＝0.02902/0.23601P＝0.744
1-pf(0.1951/0.1149,1,10)
F1,10＝0.1951/0.1149P＝0.222
28
×3
シカの効果
地域の効果
シカの効果の地域差
個体差
×6
シカの効果
地域の効果
個体差＋シカの効果の地域差
×3
×3
シカの効果
地域の効果＋シカの効果の地域差
個体差
×6
×6
シカの効果
地域の効果＋シカの効果の地域差＋個体差
Randomized complete block乱塊
Completely randomized factorial直交
Nested入れ子
Completely randomized独立
29
一般化線形混合モデル（GLMM）
• glmmML() in library(glmmML)glmと似た感覚で使える。binomialとpoissonのみ、ランダム効果はひとつ。
• glmer() in library(lme4)様々な分布、複数のランダム効果も扱える。
• glmmadmb() in library(glmmADMB)更に様々な分布でGLM・GLMMの解析可能。計算時間やや長い。
30
• 今までのは正規分布を仮定した、一般線形混合モデル（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”)
31
例題：シカがいるとアオキは減るか
• シカの分布（3地域）、非分布（2地域）で各地域5株に関して周囲の♀株の数を測定（Nested design）
非負整数値
family=poisson
32
atago fudago godai hiratuka kimikame
05
1015
site
fem
ale
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のほうが良い
33
推定
シカ在期待密度：λ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.3Models: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×対数尤度の差（尤度比）
「見たい効果を抜いたモデルの推定値の元でデータを発生させ、実データと同様の解析を行い、尤度比を計算させる」を繰り返し行い（1000回とか）、近似によらない尤度比の分布を得る方法。計算に時間がかかる。（詳しい方法はスクリプト参照）
シカの効果ありそう
34
実はPB法ではP=0.08～0.09になり、有意ではない
まとめ
• 個々の効果の強さに興味が有る→固定効果
• ばらつきを考慮したいだけ→変量効果
• 地域など空間的なものだけでなく、個体差、時間なども同様にモデリング可（同一個体の反復測定など）
• lmer()は正規分布も使えるが、近似による尤度比検定のP値はaov()と異なる場合がある。素直にaov()で示したほうが分散分析表が得られ、結果がわかりやすい
• ややこしい階層構造を持つモデルは、無理にGLMMでやるより、階層ベイズの方が推定しやすいかも
35
参考文献（GLM・GLMM）• Statistics: An introduction using R
• 日本語訳あり（共立出版）
• aov()レベルまでなら• ランダム要因の解説とかわかりやすい
• Extending the linear model with R：Generalized Linear, Mixed Effects and Nonparametric Regression Models
• glm()、lmer()あたりの解説詳しい
• 初学者には難しいかも
• 実例多く、使う人にはおすすめ
36
