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 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る