「統計学関連なんでもあり」の過去ログ--- 043

No.13407　Rでの信頼区間の算出について。　　【松本】　2010/09/14(Tue) 16:30

大変お世話になっております。

Rで多重ロジスティック回帰分析を行い，オッズ比の95%信頼区間を算出する際，confint関数を用いた場合と，exp(回帰係数 +- 1.96 * 標準誤差)という式で計算した場合とで結果が微妙に異なってしまいます。
結果が異なる理由について，全く見当がつきませんので，申し訳ありませんがご教示いただけますと幸いです。
よろしくお願い申し上げます。

以下は例です。

# infertという組み込みデータを用いて，多重ロジスティック回帰を行う。
result <- glm(case ~ age + induced + spontaneous, data = infert, family = binomial(link = "logit"))

# オッズ比
exp(coef(result))

# confintを用いた場合の95%信頼区間
exp(confint(result))

# exp(回帰係数 +- 1.96 * 標準誤差)という式から算出した95%信頼区間
se <- summary(result)$coefficients[, "Std. Error"]
upper <- exp(coef(result) + 1.96 * se)
lower <- exp(coef(result) - 1.96 * se)
cbind(lower, upper)

No.13408　Re: Rでの信頼区間の算出について。　　【青木繁伸】　2010/09/14(Tue) 17:11

自分で「exp(回帰係数 +- 1.96 * 標準誤差)という式から算出した95%信頼区間」とした場合は，exp(stats:::confint.default(result)) とするのと同じです。

exp(confint(result)) で実際に使われるのは MASS:::confint.glm です。その中では，profile 関数が result を引数として呼ばれ，その結果が result に置き換えられ，その後，stats:::confint.default が呼ばれます。なぜ，profile 関数を呼んで，その結果を confint に使うかは，? MASS:::confint.glm で出てくるオンラインヘルプに若干ですが書いてあります。
Details

confint is a generic function in package stats.

These confint methods calls the appropriate profile method, then
finds the confidence intervals by interpolation in the profile
traces. If the profile object is already available it should be
used as the main argument rather than the fitted model object
itself.
以下は，実際にどの関数が使われるか，フルネームの関数指定で書いてみたもの
> exp(confint(result))
Waiting for profiling to be done...
                 2.5 %    97.5 %
(Intercept) 0.01308397 0.5797226
age         0.96636758 1.0806991
induced     1.03050638 2.3240414
spontaneous 2.24380840 5.1933476

> exp(MASS:::confint.glm(result)) # MASS:::confint.glm が使われている
Waiting for profiling to be done...
                 2.5 %    97.5 %
(Intercept) 0.01308397 0.5797226
age         0.96636758 1.0806991
induced     1.03050638 2.3240414
spontaneous 2.24380840 5.1933476

> # exp(回帰係数 +- 1.96 * 標準誤差)という式から算出した95%信頼区間
> se <- summary(result)$coefficients[, "Std. Error"]
> z <- qnorm(0.975)                                  # 1.96 を正確に
> upper <- exp(coef(result) + z * se)
> lower <- exp(coef(result) - z * se)
> cbind(lower, upper)
                 lower     upper
(Intercept) 0.01365093 0.5969427
age         0.96641450 1.0803132
induced     1.02973800 2.3147010
spontaneous 2.21749757 5.1168090

> exp(stats:::confint.default(result))    # default を使ったのと同じ
                 2.5 %    97.5 %
(Intercept) 0.01365093 0.5969427
age         0.96641450 1.0803132
induced     1.02973800 2.3147010
spontaneous 2.21749757 5.1168090

No.13412　Re: Rでの信頼区間の算出について。　　【松本】　2010/09/15(Wed) 09:47

早速のご教示ありがとうございます。
大変，すっきり致しました。
今回の場合，confintではprofile関数の結果が引数として計算されていたのですね。

以下，自分なりにまとめてみました。

glmクラスのオブジェクトをconfintに与えると，MASSパッケージのconfint.glmが使われる。
confint.glmではまずprofile関数が実行され，その結果を引数として，statsパッケージのconfint.defaultが計算される。
profile関数を用いる理由は，非線形回帰において標準誤差はパラメータの推定における不確実性を要約するものとしては不十分で，非対称な信頼区間の方がより適切なことがあるためである（参考: 「S-PLUSによる統計解析第2版」 p259-261）。

大変，勉強になりました。
ありがとうございました。
今後とも，どうぞよろしくお願い申し上げます。