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

Size: px

Start display at page:

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

なごみとどろき
5 years ago
Views:

1 1 環境統計学ぷらす第 5 回一般 ( 化 ) 線形混合モデル高木俊 shun.takagi@sci.toho-u.ac.jp 2013/11/21

2 2 予定第 1 回 : Rの基礎と仮説検定第 2 回 : 分散分析と回帰第 3 回 : 一般線形モデル交互作用第 4.1 回 : 一般化線形モデル第 4.2 回 : モデル選択 (11/29?) 第 5 回 : 一般化線形混合モデル第 6 回 : 多変量解析 (12/5?)

3 3 今日やること統計編変量効果と混合モデル実験デザイン一般化線形混合モデル

4 4 一般化線形混合モデル

5 5 データの独立性今まで扱ってきたモデルはデータの独立性を仮定例 ) 餌種が昆虫の成長量に与える影響を見たい野外で虫を20 匹採ってきて 10 匹に餌 Aを 10 匹に餌 Bを与えて成長量を見る餌種と成長量の関係に対する 20 個の独立なデータとみなせる N= 餌 A 餌 B

6 6 偽の繰り返し (pseudoreplication) 独立性が満たされない状況野外で虫を 4 匹採ってきてそれぞれ 5 匹の子供が産まれた ( 子の合計 20) 親 A B から産まれた 10 匹に餌 A を親 C D から産まれた 10 匹に餌 B を与えた同親由来の子どもは成長も似てきそう 20 個の独立なデータとは言えない!( 真の繰り返しは 4) A B C D 親餌 A 餌 B

7 7 どのように解析するか成長 = 処理 + 誤差 = + = MS trt 処理間の分散 MS resid 処理内の分散 ( 残差 ) 成長 = 処理 + 親間の差 + 誤差 = + ( ) + = MS trt 処理 i の効果処理 i の中の j 番目の親の効果 MS parent 処理間の分散親間 ( 処理内 ) の分散

8 8 どのように解析するか成長 = 処理 + 親間の差 + 誤差 = + ( ) + = MS trt 処理 iの効果処理 iの中のj 番目の親の効果 MS parent 処理間の分散親間 ( 処理内 ) の分散 Totalばらつき = 処理間ばらつき + ( 処理内 ) 親間ばらつき +( 個体間 ) 誤差処理間ばらつきが親間ばらつきに比べて十分大きければ処理の効果がある

9 9 変量効果と混合モデル処理の効果はまさに i 番目の処理がどのように影響するかに注目固定効果 (fixed effect) 親 j の効果は個々の効果そのものではなくそのバラつきに注目変量効果 ( ランダム効果 random effect) 両者を含むモデルが混合モデル (Mixed effect model) 一般線形混合モデル一般化線形混合モデル

10 10 練習問題 :aov による解析 data5.1<- read.csv("data5.1.csv",t) data5.1$parent<- as.factor(data5.1$parent) # カテゴリカル変数として認識させる 1. 親の効果を無視して解析 summary(aov(growth~trt,data5.1)) 2. 親の効果を変量効果として解析 summary(aov(growth~trt+error(parent),data5.1)) 3. 親の効果を固定効果として解析するのと何が違う? summary(aov(growth~trt+parent,data5.1)) 4. 親ごとに平均値を計算して最初からデータ数 4 で解析するのと何が違う? par.growth<-tapply(data5.1$growth,data5.1$parent,mean) par.trt<- factor(c("c","c","t","t")) summary(aov(par.growth~par.trt))

11 11 出力結果 > summary(aov(growth~trt+error(parent),data5.1)) Error: parent Df Sum Sq Mean Sq F value Pr(>F) trt Residuals 親を誤差とした分散分析表この Residual は ( 処理内 ) 親間の誤差親間の誤差と比較すると処理の効果は有意ではない (P=0.1) Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals こっち Residual は ( 親内 ) 個体間の誤差個体を誤差とした分散分析表

12 12 ここまでのまとめデータが独立ではなく同一親同一地点といったデータのばらつきに影響する要因が入っている場合それらを変量効果として考慮する必要がある変量効果を無視すると個々のデータは偽の繰り返しとなり差がないものが有意になる (Type I Error) 差があるものが検出できない (Type II Error) ことがある変量効果を入れるかどうかどのように入れるかはデータの取り方 ( 実験デザイン調査デザイン ) による

13 13 実験デザイン Completely randomized Nested Randomized complete block (unreplicated) Completely randomized factorial Split-plot 実験調査デザインで解析も変わってくる

14 14 Nested design 要因 A の中に要因 B が入れ子繰り返しあり ( ) ( ) ( ) ( ) 例 : ヒシ帯 2 地点開放水面 2 地点がありそれぞれの地点で 2 箇所ずつ測定 ( ) ( ) ( ) = + ( ) + ( ) ( ) ( ) ( ) ヒシの効果地点間誤差地点内場所間誤差 A:fixed B:random ( ) ( ) ( ) = MS A ヒシvs 開放間の分散 MS ( ) 地点間の分散

15 15 Nested design 例 : ヒシ帯 2 地点開放水面 2 地点がありそれぞれの地点で 2 箇所ずつ測定 ( 観測数 ) Site1 Site2 Site3 Site4 Open 2 2 Trapa 2 2 ヒシ vs 開放水面間の分散 Trapa ( ) 地点間の分散 Site1 Site2 Site3 Site4 Open Trapa

16 16 Randomized block design( 乱塊法 ) 要因 A と要因 B は直交繰り返しなし例 : 4 ブロックがありそれぞれのブロックでヒシ刈り取り刈り残し 1 地点ずつ測定 = + + ヒシの効果ブロック間誤差いずれでも説明されない残差 A:fixed B:random = MS A MS resid 刈り取り vs 刈り残し間の分散処理ブロックで説明されない分散 ( 残差 )

17 17 Randomized block design 例 : 4 ブロックがありそれぞれのブロックで刈り取り刈り残し 1 地点ずつ測定 ( 観測数 ) Block1 Block2 Block3 Block4 刈取 (O) 刈残 (T) 刈り残しvs 刈り取り間の分散 Trapa 説明できない誤差 Blk1 O T Blk2 Blk3 Blk4 O T O T O T

18 18 Factorial design 要因 A と要因 B は直交各組み合わせに繰り返しあり例 : 2 ブロックがありそれぞれのブロックで刈り取り刈り残し 2 地点ずつ測定いずれでも説明されない残差 = ヒシの効果ブロック間誤差 A:fixed B:random ヒシの効果のブロック間誤差 = MS A MS AB 刈り取り vs 刈り残し間の分散刈り取りの効果のブロック間でのばらつき

19 19 Factorial design 例 : 2 ブロックがありそれぞれのブロックで刈り取り刈り残し 2 地点ずつ測定 ( 観測数 ) Block1 Block2 刈取 (O) 2 2 刈残 (T) 2 2 刈り残し vs 刈り取り間の分散 : 刈り取り効果のブロック間での分散 O Block1 T O Block2 T

20 20 実験デザインによる ANOVA の違い Completely randomized ( 独立 ) Nested ( 入れ子 ) = + 処理残差 = + ( ) + 処理内ブロック = = MS A MS resid 処理の効果 / 処理で説明できない残差 MS A MS ( ) 処理の効果 / 処理内のブロック間のばらつき Randomized block ( 乱塊 ) = + + ブロック = MS A MS resid 処理の効果 / 処理とブロックの効果を除いた残差 Factorial ( 直交 ) = 処理 : ブロック交互作用 = MS A MS AB ( ブロック共通の ) 処理の効果 / ブロックごとの処理の効果のばらつき

21 21 練習問題 2: シカの効果と調査デザイン A) シカのいる3 地域のシカ柵内と柵外で2 個体の樹高測定 B) シカのいる6 地域のシカ柵内と柵外で1 個体の樹高測定 C) シカのいる3 地域いない3 地域で2 個体ずつ樹高測定 D) シカのいる6 地域いない6 地域で1 個体ずつ樹高測定 A) シカ 6 コントロール 6 のデザインでも全部違う! C) B) 3 地域 3 地域 3 地域 D) 6 地域 6 地域 6 地域

22 22 R 上でのデータ記述法の違い A) B) C) D) 3 6 deer site deer site deer site deer site

23 23 解析前の注意点シカの在不在地点 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)

A) 直交 Completely randomized factorial design = + + + 3 固定効果変量効果変量効果残差 ( 交互作用 ) summary(aov(height~deer+error(site/deer),dat5.2a)) 24 Error: site Df Sum Sq Mean Sq F value Pr(>F) Residuals 2 0.

24 A) 直交 Completely randomized factorial design = 固定効果変量効果変量効果残差 ( 交互作用 ) summary(aov(height~deer+error(site/deer),dat5.2a)) 24 Error: site Df Sum Sq Mean Sq F value Pr(>F) Residuals site で説明されるばらつき Error: site:deer Df Sum Sq Mean Sq F value Pr(>F) deer * Residuals Signif. codes: 0 *** ** 0.01 * deer で説明されるばらつき deer:site で説明されるばらつき Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 残り ( 個体差 )

25 B) 乱塊 Randomized complete block design 25 = 固定効果変量効果残差 summary(aov(height~deer+error(site),dat5.2b)) Error: site Df Sum Sq Mean Sq F value Pr(>F) Residuals site で説明されるばらつき Error: Within Df Sum Sq Mean Sq F value Pr(>F) deer * Residuals Signif. codes: 0 *** ** 0.01 * deer で説明されるばらつき残り ( 個体差かつ効果の地域差 )

26 C) 入れ子 Nested design 3 = ( ) とも表記 (RCB と区別 ) 3 固定効果変量効果残差 summary(aov(height~deer+error(site),dat5.2c)) Error: site Df Sum Sq Mean Sq F value Pr(>F) deer Residuals Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals deer で説明されるばらつき site で説明されるばらつき ( かつ効果の地域差 ) 残り ( 個体差 )

D) 独立 Completely randomized design 27 6 6 = + 固定効果残差 summary(aov(height~deer,dat5.

27 D) 独立 Completely randomized design = + 固定効果残差 summary(aov(height~deer,dat5.2d)) Df Sum Sq Mean Sq F value Pr(>F) deer Residuals deer で説明されるばらつき残り ( 個体差かつ地域差かつ効果の地域差 )

28 28 何も考えず aov(height~deer*site) で解析すると > summary(aov(height~deer*site,data5.2a)) Df Sum Sq Mean Sq F value Pr(>F) deer ** site *** deer:site Residuals 左の表からも正しいF,Pを計算可能 1-pf(0.1315/0.0027,1,2) F1,2=0.1315/ P=0.020 > summary(aov(height~deer*site,dat5.2b)) Df Sum Sq Mean Sq deer site deer:site pf( / ,1,5) F1,5= / P=0.018 > summary(aov(height~deer*site,dat5.2c)) Df Sum Sq Mean Sq F value Pr(>F) deer site *** Residuals > summary(aov(height~deer*site,dat5.2d)) Df Sum Sq Mean Sq deer site pf( / ,1,4) F1,4= / P= pf(0.1951/0.1149,1,10) F1,10=0.1951/ P=0.222

29 29 Completely randomized factorial 直交 Randomized complete block 乱塊 3 6 シカの効果地域の効果シカの効果の地域差個体差シカの効果地域の効果個体差 +シカの効果の地域差 Nested 入れ子 3 3 シカの効果地域の効果 + シカの効果の地域差個体差 Completely randomized 独立 6 6 シカの効果地域の効果 + シカの効果の地域差 + 個体差

30 30 一般化線形混合モデル (GLMM) 今までのは正規分布を仮定した一般線形混合モデル (LMM と書くことも ) 一般線形混合モデル (LMM) で正規分布以外を仮定する一般化線形混合モデル (GLMM) も存在 R において GLMM の解析が可能な関数 glmmml() in library(glmmml) glm と似た感覚で使える binomial と poisson のみランダム効果はひとつ glmer() in library(lme4) 様々な分布複数のランダム効果も扱える glmmadmb() in library(glmmadmb) 更に様々な分布で GLM GLMM の解析可能計算時間やや長い #LMM は lmer で解析 #R のバージョンによっては website からのインストール

31 31 モデル式はほぼ同じただ最尤法を使っての推定なので計算は複雑一般線形混合モデル aov(y~trt+error(site)) = + 処理の効果一般化線形混合モデル log λ = + 場所間ばらつき ~ (0, 2 ) ~ (, 2 ) glmmml(y~trt,cluster=site,family=poisson) glmer(y~trt+(1 site),family=poisson) glmmadmb(y~trt+(1 site),family= poisson ) ~ (0, 2 ) ~ (λ )

32 32 例題 : シカがいるとアオキは減るかシカの分布 (3 地域 ) 非分布 (2 地域 ) で各地域 5 株に関して周囲の株の数を測定 (Nested design) female 非負整数値 family=poisson atago fudago godai hiratuka kimikame site

33 33 推定 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 Random effects: Groups Name Variance Std.Dev. site (Intercept) Number of obs: 25, groups: site, 5 log λ = + ~ (0, 2 ) Fixed effects: Estimate Std. Error z value Pr(> z ) (Intercept) *** deernd ** --- Signif. codes: 0 *** ** 0.01 * シカ在期待密度 : λd = exp(1.0085) シカ不在期待密度 : λnd = exp( ) Correlation of Fixed Effects: (Intr) deernd 係数の推定値に対する Wald 検定による P 値変数の効果をみるなら尤度比検定 Parametric Bootstrap のほうが良い

34 34 検定カイ二乗近似を用いた尤度比検定 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) 見たい効果抜き (null) glmer.model5.3: female ~ deer + (1 site) 見たい効果入り (deer) Df AIC BIC loglik deviance Chisq Chi Df Pr(>Chisq) glmer.model glmer.model * --- Signif. codes: 0 *** ** 0.01 * 対数尤度の差 ( 尤度比 ) シカの効果ありそう Parametric bootstrap 法による近似によらない検定見たい効果を抜いたモデルの推定値の元でデータを発生させ実データと同様の解析を行い尤度比を計算させるを繰り返し行い (1000 回とか ) 近似によらない尤度比の分布を得る方法計算に時間がかかる ( 詳しい方法はスクリプト参照 ) 実は PB 法では P=0.08~0.09 になり有意ではない

35 35 まとめ個々の効果の強さに興味が有る固定効果ばらつきを考慮したいだけ変量効果地域など空間的なものだけでなく個体差時間なども同様にモデリング可 ( 同一個体の反復測定など ) lmer() は正規分布も使えるが近似による尤度比検定の P 値は aov() と異なる場合がある素直に aov() で示したほうが分散分析表が得られ結果がわかりやすいややこしい階層構造を持つモデルは無理に GLMM でやるより階層ベイズの方が推定しやすいかも

36 36 参考文献 (GLM GLMM) Statistics: An introduction using R 日本語訳あり ( 共立出版 ) aov() レベルまでならランダム要因の解説とかわかりやすい Extending the linear model with R:Generalized Linear, Mixed Effects and Nonparametric Regression Models glm() lmer() あたりの解説詳しい初学者には難しいかも実例多く使う人にはおすすめ

Microsoft PowerPoint - GLMMexample_ver pptx

Microsoft PowerPoint - GLMMexample_ver pptx Linear Mixed Model ( 以下混合モデル ) の短い解説この解説のPDFは http://www.lowtem.hokudai.ac.jp/plantecol/akihiro/sumida-index.html のお勉強のページにあります. ver 20121121 ととの間に次のような関係が見つかったとしよう全体的な傾向に対する回帰直線を点線で示したところがこれらのデータは実は異なる

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

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