No.21454　ロジスティック回帰モデル(説明変数が質的変数)　　【tani】　2014/11/29(Sat) 20:59

お世話になります。R初心者でございます。
以下のようなデータから，ランダム効果を考慮したロジスティック回帰モデルを作成するコードを書いているのですが，期待通りの結果にならず，悩んでおります。

> X
center group n sideeffect
1 1 A 32 14
2 1 B 33 18
3 2 A 30 4
4 2 B 28 8
5 3 A 23 14
6 3 B 24 9
7 4 A 22 7
8 4 B 22 10
9 5 A 20 6
10 5 B 21 12
11 6 A 19 1
12 6 B 20 3
13 7 A 17 2
14 7 B 17 6
15 8 A 16 7
16 8 B 15 9
17 9 A 13 1
18 9 B 14 5
19 10 A 13 3
20 10 B 13 1
21 11 A 11 1
22 11 B 12 2
23 12 A 10 1
24 12 B 9 0
25 13 A 9 2
26 13 B 9 6
27 14 A 8 1
28 14 B 8 1
29 15 A 7 1
30 15 B 8 0

library(glmmML)
m <- glmmML(cbind(sideeffect,n-sideeffect)~group,family=binomial,data=X,cluster=center)

summary(m)

Call: glmmML(formula = cbind(sideeffect, n - sideeffect) ~ group, family = binomial, data = X, cluster = center)

coef se(coef) z
(Intercept) -1.3370 0.2667 -5.014
groupB 0.4995 0.2057 2.428
Pr(>|z|)
(Intercept) 5.34e-07
groupB 1.52e-02

Scale parameter in mixing distribution: 0.7826 gaussian
Std. Error: 0.2011

LR p-value for H_0: sigma = 0: 3.759e-08

Residual deviance: 56.61 on 27 degrees of freedom AIC: 62.61
>

ここで，要因groupの水準Aにおける係数の推定値が計算されないのは何故でしょうか？
コードが間違っているでしょうか？

よろしくお願いいたします。
No.21455　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【青木繁伸】　2014/11/29(Sat) 21:05

> 要因groupの水準Aにおける係数の推定値が計算されないのは何故でしょうか？

要因groupの水準Aを基準とするためです。つまり，係数の推定値は0です。
要因groupの水準Aを0としたときに，要因groupの水準Bはいくつかということですね。
No.21456　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【tani】　2014/11/29(Sat) 21:28

青木先生

早速回答頂き，ありがとうございます。
では，水準Bを基準とする場合はどのように書けばいいのでしょうか？
No.21457　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【波音】　2014/11/29(Sat) 22:47

本質的な回答ではないと思いますが，，

例えばlm()の場合，

lm(Y ~ Group, contrasts=list(Group="contr.SAS"))

といったように引数に　contrasts=list(Group="contr.SAS")　という指定をしてあげることで，最後の水準がベースラインとして計算されます。

glmMML()でも同じようにしたらどうでしょうか？

ただし，これはあくまでも「最後の水準をベースラインにする」ということなので，今回はAとBの2水準しかないので結果として，Bがベースラインにおかれるという発想にすぎません。

なおcontr.SASとはSASで採用されている計画行列を用いるという意味です。
（難しい話なのでおまじないとでも思っておけばOKかと）
No.21458　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【tani】　2014/11/29(Sat) 23:55

波音様

回答ありがとうございます。そうです！SASの解析結果との整合性を取りたいのです。

glmmML関数の場合，引数に指定するのではなく，下記のようにすることで上手く動くようです。

options("contrasts"=c("contr.SAS","contr.SAS"))
library(glmmML)
m <- glmmML(cbind(sideeffect,n-sideeffect)~group,family=binomial,data=X,cluster=center)

頂いた回答をヒントに検索をかけて，解決致しました。
ありがとうございました。
No.21459　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【青木繁伸】　2014/11/29(Sat) 23:57

ダミー変数を使う分析は R では，最初のカテゴリーを基準とします。では，何が最初のカテゴリーになるかといえば，factor 関数で明示的に指定すればそこで指定された最初のカテゴリーが基準になります。では，R は特に指定しないときにはカテゴリーの順序をどのように解釈するかということですが，それは，単にカテゴリーに与えられた名前の辞書順ということです。カテゴリーの名前が辞書でもっとも早く出現する順ということ。アルファベットならば単純ですが，そうでない場合は納得のいかないこともあるかもしれないですが，とにかくコンピュータの内部ではコンピュータが準拠する辞書順です。ということで，納得のいかない順序づけをされたくないならば，factor 関数でカテゴリーの順序づけをするのがお勧めです。
No.21460　Re: ロジスティック回帰モデル(説明変数が質的変数)　　【tani】　2014/12/01(Mon) 00:11

青木先生

X$group <- factor(X$group, levels=c("B", "A"))
と書けば，最初の水準Bが基準となる，ということですね。

ありがとうございます。助かりました。

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