No.21265 Rを使ってGLMで温度と生残数(死亡数)の検定  【フラ】 2014/08/20(Wed) 12:03

Windows7,64bit,R3.1.1を使っています。

温度と生残数(死亡数)に有意な関係があるか調べています。
以下のようなデータ(Levelのアルファベットは温度の違いだと思って下さい)を用意し,test.txtという名前で保存しました。
Level	dead	live
A 50 0
B 40 10
C 35 15
D 22 28
E 11 39
F 4 46
G 1 49
H 9 41
I 20 30
J 30 20
K 50 0
そして,Rのコンソール画面に,以下のようにコマンドを入力しました。
> test1<-read.table("test.txt",header=T)
> result<-glm(cbind(dead,live)~Level,family=binomial,data=test1)
> summary(result)
すると以下の結果が得られました。
Call:
glm(formula = cbind(dead, live) ~ Level, family = binomial, data = test1)

Deviance Residuals:
[1] 0 0 0 0 0 0 0 0 0 0 0

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.663e+01 5.202e+04 0.001 1
LevelB -2.524e+01 5.202e+04 0.000 1
LevelC -2.578e+01 5.202e+04 0.000 1
(中略)
LevelK 7.732e-08 7.356e+04 0.000 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 3.1056e+02 on 10 degrees of freedom
Residual deviance: 5.4385e-10 on 0 degrees of freedom
AIC: 56.184

Number of Fisher Scoring iterations: 22
Pr(>|z|)が全て1になるのはおかしいと思いましたが,multcompパッケージを使って以下の様に多重比較もしてみました。
> mult.test<-glht(result,linfct=mcp(Level="Tukey"))
> summary(mult.test)
結果の一部を下に貼り付けましたが,AとGですら有意差が無いと言う結果になりました。
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 -2.524e+01 5.202e+04 0.000 1.0000
C - A == 0 -2.578e+01 5.202e+04 0.000 1.0000
D - A == 0 -2.687e+01 5.202e+04 -0.001 1.0000
E - A == 0 -2.790e+01 5.202e+04 -0.001 1.0000
F - A == 0 -2.907e+01 5.202e+04 -0.001 1.0000
G - A == 0 -3.052e+01 5.202e+04 -0.001 1.0000
もしデータにゼロが入っていない場合は(例えば↓のようなデータ),おかしな結果にはならなかった(と思う)のですが(結果割愛),このような検定(glmで生残比較→glhtで多重比較)では,ゼロが入っていたらダメというようなことはあるのでしょうか。
それとも,私の手順に根本的な間違いがあるのでしょうか。
もしくは,素直にカイ二乗検定を使うべきなのでしょうか。
Level	dead	live
A 49 1
B 40 10
C 35 15
D 22 28
E 11 39
F 4 46
G 1 49
H 9 41
I 20 30
J 30 20
K 42 8
よろしくお願いいたします。

No.21266 Re: Rを使ってGLMで温度と生残数(死亡数)の検定  【青木繁伸】 2014/08/20(Wed) 13:44

> もしデータにゼロが入っていない場合は(例えば↓のようなデータ),おかしな結果にはならなかった(と思う)のですが

そのようなデータの場合は,確かに,見た目が明らかに変という結果ではないですが,根本的な問題は,
Residual deviance: 5.4385e-10 on 0 degrees of freedom
の行に示されています。ゼロが入っていようがいまいが,"0 degrees of freedom" になります。これは,モデルが完全にあてはまっているということを示しています。

あなたがこの解析を行うときに参考にしたかもしれない例は,
glm(cbind(dead, live) ~ DEGREE, family=binomial, data=test1)
のようになっていませんでしたか。DEGREE は数値データだったはずです。
そのようにした場合は,
Residual deviance: xx.xxx on 9 degrees of freedom
となるはずです。自由度 9 は,データの個数から,「予測に使った変数の個数+1」を引いたものになるのです。9 = 11 - (1 + 1)
さて,あなたは,数値変数ではなく,数値変数に対応する記号を使いました。R はそのような変数は Factor 変数として扱います。そして,Factor 変数を lm や glm などに使うと,R はそれを「Factor 変数が表すカテゴリー数 - 1」個のダミー変数として使います。どのようなダミー変数になったかは分析結果に示されていますが,LevelB 〜 LevelK の 10 変数です。なので,あなたの分析結果では Residual deviance の自由度は 11 - (10 + 1) = 0 になってしまったのです。
では,なぜ自由度が 0 ではおかしいのか,「モデルが完全にあてはまっている」のがなぜおかしいのかを示しましょう。二次元座標に 2 つの点 (x1, y1) と (x2, y2) があります。データが 2 個あるということですね。このデータに y = a x + b という直線回帰をしたとします。a, b は推定値ではなく確定値です... 2 点を通る直線は一意に決まる。
> d <- data.frame(x=c(1,5), y=c(3,7))
> ans <- lm(y ~ x, d)
> summary(ans)

Call:
lm(formula = y ~ x, data = d)

Residuals:
ALL 2 residuals are 0: no residual degrees of freedom!

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2 NA NA NA
x 1 NA NA NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: NaN
F-statistic: NaN on 1 and 0 DF, p-value: NA
NA not available,NaN は not a number です。Residual standard error(Residual deviance に対応)の自由度が 0,F-statistic の第 2 自由度も 0 になっています。
あなたの glm による分析も同じ問題を持ってしまうのです。

じゃあ,どうすればよいのか?
あなたが「Levelのアルファベットは温度の違い」といっているのだから,「温度」を使って分析すればよいでしょうということ。「温度」が正確に数値データでなくても,それは「低温」,「中程度」,「高温」などのように Factor で表してもよいです。しかし,その Factor 変数が取る値は,「データの個数 - 1個」未満でなければいけません。

No.21268 Re: Rを使ってGLMで温度と生残数(死亡数)の検定  【フラ】 2014/08/20(Wed) 22:35

青木先生 様

ご回答いただき,ありがとうございます。
「ゼロが入ってなかったら検定できるのに・・・」と思っていたことが,そもそもの間違いだったことがよく分かりました。
具体例を示していただいたことで,さらによく理解出来ました。
素直に,温度を数値データとして検定します。

この度はお忙しいところ,まことにありがとうございました。

● 「統計学関連なんでもあり」の過去ログ--- 047 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る