GLM <- glm(cbind(Dis,Health) ~ TRT + REP , family=binomial(logit), data=Dataset)使用したデータは以下のとおりです。
Anova(GLM, type="II")
Response: cbind(Dis, Health)
LR Chisq Df Pr(>Chisq)
TRT 6.7399 1 0.009428 **
REP 3.1048 1 0.078064 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
試験場所 管理方法 ブロック 発病 健全
PL TRT REP Dis Health
1 1 1 1 34 12
2 1 1 2 30 12
3 1 1 3 28 9
4 1 2 1 14 28
5 1 2 2 13 36
6 1 2 3 24 22
7 1 3 1 21 24
8 1 3 2 28 16
9 1 3 3 25 15
No.08731 Re: 反復がある多群の比率の差の検定について 【知ったかぶり】 2008/12/24(Wed) 17:12
管理方法(TRT)とブロック(REP)が数値データのままです.まずは,因子型に変換しましょう.
多重比較に関しては,管理方法の全てのペア毎にロジスティック回帰を行って,Bonferroniで有意水準を調整するのが一番簡単ではないかと思います.検定の回数は全部で3ですから,検出力の低下をあまり気にしなくてもよいのでは.
地域による違いを検討したいのであれば,試験場所(PL)を説明変数に加えてみましょう(その際,ブロックの水準のつけかたに注意.試験場所毎にブロック番号は変えなければいけない).
No.08733 Re: 反復がある多群の比率の差の検定について 【青木繁伸】 2008/12/24(Wed) 22:56
一部結果が再現できないのですが，お使いいただいているRのバージョンなどを報告してください。また，タイプミス というか，実行結果をコピー・ペーストしていない場合には，間違いがないか，再確認していただけると幸いかと。。。Anova って，どこで定義されている関数でしょうかということなのですけど。。。
No.08735 Re: 反復がある多群の比率の差の検定について 【はた】 2008/12/25(Thu) 01:27
知ったかぶり様，青木繁伸様
お返事ありがとうございます。
青木 様
Rのバージョンは2.8.0です。
改めて，同一データを処理したところ，前記の出力結果と違っていました。オブジェクトのクリアができていなかったのでしょうか。下記は再実施した経過ですが，今度はいかがでしょうか？
Anovaについては，Rコマンダーメニューの「モデル」「仮説検定」「分散分析表」で検定のタイプをデフォルトの「Type�供廚鯀�択しました。
知ったかぶり 様
ロジスティック回帰とは，最初のRコマンドの関数glmによるロジスティック回帰でよいでしょうか？また，多重比較は以下のようにすると9区総当たりの対 比較となってしまいます。管理方法(TRT)による群間での検定はどのようにすれば良いでしょうか？以下のRコマンドではデフォルトでHolmの方法がと られていますが，Bonferroniの方法を選択すべきなのでしょうか？> print(pairwise.prop.test(Dis,Dis+Health))3試験場所でのブロック（REP)番号は，以下のように変更するということでしょうか？試験場所 管理方法 ブロック
PL TRT REP
1 1〜3 1〜3 → 1〜3
2 1〜3 1〜3 → 4〜6
3 1〜3 1〜3 → 7〜9
------------------------------------------------------
（再実施した出力）
> attach(Dataset)
> PL <- as.factor(PL)
> REP <- as.factor(REP)
> TRT <- as.factor(TRT)
> GLM <- glm(cbind(Dis, Health) ~ TRT + REP, family=binomial(logit),
+ data=Dataset)
> summary(GLM)
Call:
glm(formula = cbind(Dis, Health) ~ TRT + REP, family = binomial(logit), data = Dataset)
Deviance Residuals:
1 2 3 4 5 6 7 8
0.62171 -0.06081 -0.69741 0.19819 -1.14686 0.92353 -0.74881 1.19295
9
-0.46186
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.8360 0.2448 3.414 0.00064 ***
TRT[T.2] -1.5942 0.2726 -5.849 4.94e-09 ***
TRT[T.3] -0.7461 0.2718 -2.745 0.00606 **
REP[T.2] 0.1011 0.2576 0.393 0.69464
REP[T.3] 0.5728 0.2691 2.129 0.03330 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 46.5522 on 8 degrees of freedom
Residual deviance: 5.2812 on 4 degrees of freedom
AIC: 52.104
Number of Fisher Scoring iterations: 3
> Anova(GLM, type="II")
Anova Table (Type II tests)
Response: cbind(Dis, Health)
LR Chisq Df Pr(>Chisq)
TRT 37.560 2 6.98e-09 ***
REP 5.172 2 0.07531 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No.08736 Re: 反復がある多群の比率の差の検定について 【知ったかぶり】 2008/12/25(Thu) 09:31
Anova()は,library(car)の関数ですね.Rコマンダーではこれがデフォルトなんですか…なんでだろう？null modelとの比較ならanova()でもよさそうですが.
> print(pairwise.prop.test(Dis,Dis+Health))
これは,比率の差の多重比較であり,一般化線形モデル(この場合は,ロジスティック回帰）とは「枠組」が違います.両者を併用するのはマズイでしょう.
研究の目的にもよりますが,多重比較を考える前にモデル選択(変数選択)を行うべきでは？説明変数を全て投入したモデルを作成して,stepAIC()関数を使えば,AICによるモデル選択を行うことができます.
library(MASS)
GLM <- glm(cbind(Dis, Health) ~ PL + TRT + REP, family=binomial(logit), data=Dataset)
stepAIC(GLM)
多重比較は,どういうモデルが選択されるか,はたさんが何と何を比較したいか,によると思います.必ずしも9区総当たりの検定を行う必要はないでしょう.
＞3試験場所でのブロック（REP)番号は，以下のように変更するということでしょうか？
その通りです.異なる試験場所に同じブロック番号を付けると,「同じ」ブロックとして処理されてしまいます.
No.08752 Re: 反復がある多群の比率の差の検定について 【はた】 2008/12/25(Thu) 19:14
知ったかぶり 様
ご教示ありがとうございます。
この試験では，圃場でのばらつきを考慮 し，前作での管理方法（TRT：3水準)の違いにより，後作の植物の発病が少なくなるのか。また発病を少なくさせるためには，どの管理方法を選択すべきか を明らかにしたいと思っております。さらに，標高が異なっても同様なこと（安定した効果）が言えるかということです。
3試験場所での データについて説明変数（PL，TRT，REP)を全て投入したモデルにおいて，関数stepAIC()を行うと，TRT＋REPを説明変数としたモデル が選択されました。また，分散分析表でも説明変数として，下記のようにTRT，REPそれぞれが有意である結果となっています。この場合，TRTの効果に ついて多重比較を行いたいところですが方法が判りません。また，REPの影響も大きいようですのでどのように解析を進めて行ったらよいのか困っています。 よろしくご教示下さい。> GLM2 <- stepAIC(GLM)
Start: AIC=139.35
cbind(Dis, Health) ~ PL + TRT + REP
Step: AIC=139.35
cbind(Dis, Health) ~ TRT + REP
Df Deviance AIC48.52 139.35
- TRT 2 118.97 205.80
- REP 8 505.35 580.19
> summary(GLM2)
省略
> Anova(GLM2, type="II")
Anova Table (Type II tests)
Response: cbind(Dis, Health)
LR Chisq Df Pr(>Chisq)
TRT 70.45 2 5.551e-16 ***
REP 456.84 8 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No.08759 Re: 反復がある多群の比率の差の検定について 【知ったかぶり】 2008/12/25(Thu) 23:59
＞TRT＋REPを説明変数としたモデルが選択されました。
ということは，PLを考慮しなくて良いわけですから，全てのデータをプールしてTRTについて多重比較を行えば良いと思います．
多重比較の具体的な手順ですが，TRTの3水準について2水準ごとを組み合わせたデータセット（3組）を作成して，そのそれぞれについて，ロジスティック回帰を行えば良いのでは．Anova()は良くわからないので，anova()を使うと，
GLM<-glm(cbind(Dis,Health)~TRT+REP,family=binomial(logit),data=dataset)
GLM.REP<-glm(cbind(Dis,Health)~REP,family=binomial(logit),data=dataset)
anova(GLM,GLM.REP,test="Chisq")
と すれば，REPのみのモデルとの比較において，当該のモデルが有意に当てはまりが良いか否かの検定が行われます（尤度比検定）．TRTを含めたモデルの方 が当てはまりが良ければ，TRTの水準間には差があるということです．これを上記の3組のデータセットそれぞれについて行い，Bonferroniで有意 水準を調整すれば良いと思います（検出力を高めたいのであれば，Holmの方法などを使う．まあ，大差はないと思いますが）．
No.08767 Re: 反復がある多群の比率の差の検定について 【はた】 2008/12/26(Fri) 17:56
知ったかぶり 様
丁寧にご説明いただきありがとうございます。
ご教示いただいた方法をやってみました。また，尤度比検定という用語が出てきましたので，それを足がかりに調べますと中澤先生の著書「Rによる保健医療データ解析演習」の中にも尤度比検定として以下の方法2が示されていましたので試しに行ってみました。
ともに有意水準5%より小さい値であることから，帰無仮説が棄却されてREPのみよりもTRT+REPのモデルが選択され，TRTの水準間の比較に移行できるわけですね。
この尤度比検定の2方法は同じに使用できると考えてよろしいでしょうか？一般化線形モデルは初めてなのですが，尤度比検定については複数のやり方があるの ですね。また，変数選択に用いたstepAIC（）関数は，REPのみではAICが大きくなり，選択されなかった結果と考え，今回のような場面でも使用で きると考えることは無理なのでしょうか？方法1
> GLM<-glm(cbind(Dis,Health)~TRT+REP,family=binomial(logit),data=Dataset)
> GLM.REP<-glm(cbind(Dis,Health)~REP,family=binomial(logit),data=Dataset)
> anova(GLM,GLM.REP,test="Chisq")
Analysis of Deviance Table
Model 1: cbind(Dis, Health) ~ TRT + REP
Model 2: cbind(Dis, Health) ~ REP
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 24 529.47
2 25 534.57 -1 -5.09 0.02
方法2
> GLM<-glm(cbind(Dis,Health)~TRT+REP,family=binomial(logit),data=Dataset)
> GLM.REP<-glm(cbind(Dis,Health)~REP,family=binomial(logit),data=Dataset)
> lambda <- -2*(logLik(GLM.REP)-logLik(GLM))
> 1-pchisq(lambda,1)
'log Lik.' 0.02400585 (df=2)
No.08768 Re: 反復がある多群の比率の差の検定について 【知ったかぶり】 2008/12/27(Sat) 00:39
＞この尤度比検定の2方法は同じに使用できると考えてよろしいでしょうか？
違っているように見えるかもしれませんが，どちらの方法もやっていることは同じ．anova()を使うと関数内部で何をしているか見えないだけです．
anova(GLM,GLM.REP,test="Chisq")$P[2]
とすれば，P値だけが出力されますが，方法2の値と全く同じになるはずです．
＞変数選択に用いたstepAIC（）関数は，REPのみではAICが大きくなり，選択されなかった結果と考え，今回のような場面でも使用できると考えることは無理なのでしょうか？
「検定」にこだわらないのであれば，全て「モデル選択」というテもあります．詳しくは，こちらを．
http://hosho.ees.hokudai.ac.jp/~kubo/ce/2004/kubo/kuboJES2004.pdf
http://hosho.ees.hokudai.ac.jp/~kubo/ce/2004/kubo/R.html
No.08794 Re: 反復がある多群の比率の差の検定について 【はた】 2008/12/29(Mon) 03:39
紹介いただいたサイトはたいへん参考になりました。
さて，TRT2水準の組み合わせ3通りで以下を実施してみました。> GLM<-glm(cbind(Dis,Health)~TRT+REP,family=binomial(logit),data=Dataset)ボンフェローニの方法による有意水準の補正(全体α=0.05または0.01，総検定数3回）
> GLM.REP<-glm(cbind(Dis,Health)~REP,family=binomial(logit),data=Dataset)
> anova(GLM,GLM.REP,test="Chisq")
Analysis of Deviance Table
Model 1: cbind(Dis, Health) ~ TRT + REP
Model 2: cbind(Dis, Health) ~ REP
TRT水準(1,2)の場合
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 15 336.45
2 16 379.71 -1 -43.26 4.789e-11
TRT水準(1,3)の場合
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 15 239.361
2 16 245.161 -1 -5.800 0.016
TRT水準(2,3)の場合
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 15 403.28
2 16 420.41 -1 -17.13 3.483e-05
補正した有意水準は3組共通で α=0.05/3=0.0167，α=0.01/3=0.0033TRT水準(1,2) p=4.789e-11 ** → TRT水準1と2の間に差があるといった理解でよいでしょうか？
TRT水準(1,3) p=0.016 * → TRT水準1と3の間に差がある（微妙?）
TRT水準(2,3) p=3.483e-05 ** → TRT水準2と3の間に差がある
最初の全説明変数を投入したロジスティック回帰において，試験場所PLが組み込まれませんでしたが，データを見ると3箇所のうち1箇所は他の場所と違い， 各区ともほぼ100%が発病している状態です。この1箇所を除いてTRTの効果を再検討してみる必要があるとも考えられますが，勝手に除くのではなく，何 か分離して考えるべき根拠を示して先に進める方法はあるのでしょうか？紹介いただいたサイトでは一般化線形混合モデルは場所の差をモデルに組み込んで説明 できるようなのですが，今回のような場面で適用できるでしょうか？
● 「統計学関連なんでもあり」の過去ログ--- 042 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る