Triad sou.

glm() でオッズ比と信頼区間を求める

R

基礎的すぎるから知らなかったのか confint(ModelName, level = 1-alpha) という関数があったらしい、プロファイル尤度に基づく信頼区間を計算できる (2012/6/19 記述を修正)。

d <- c(
  1,1,1,1,0,0,0,0,0,0, 1,1,1,1,1,1,1,1,0,0,
  1,0,0,0,0,0,0,0,0,0, 1,1,1,1,1,0,0,0,0,0
)
trt <- gl(2, 20, 40)
foo <- gl(2, 10, 40)
GLM.1 <- glm(d ~ trt + foo, family=binomial())
summary(GLM.1)

ModelName <- GLM.1
alpha <- 0.05
x <- summary(ModelName)
y <- confint(ModelName, level=1-alpha)
OR <- exp(x$coefficients[,1])
LowerCL <- exp(y[,1])
UpperCL <- exp(y[,2])
z <- cbind(OR, LowerCL, UpperCL)
z

                   OR    LowerCL    UpperCL
(Intercept) 0.6228675 0.18065550  1.9492149
trt[T.2]    0.2110707 0.03920214  0.8931083
foo[T.2]    7.1291202 1.70056921 38.8209458

追記 (2012/6/19)

epicalc パッケージの logistic.display 関数も便利です。

require(epicalc)
d <- c(
  1,1,1,1,0,0,0,0,0,0, 1,1,1,1,1,1,1,1,0,0,
  1,0,0,0,0,0,0,0,0,0, 1,1,1,1,1,0,0,0,0,0
)
trt <- gl(2, 20, 40)
foo <- gl(2, 10, 40)
GLM.1 <- glm(d ~ trt + foo, family=binomial())
logistic.display(GLM.1)
logistic.display(GLM.1, simplified = TRUE)

> logistic.display(GLM.1)
Logistic regression predicting d 
 
            crude OR(95%CI)    adj. OR(95%CI)     P(Wald's test) P(LR-test)
trt: 2 vs 1 0.29 (0.08,1.06)   0.21 (0.05,0.97)   0.045          0.034     
                                                                           
foo: 2 vs 1 5.57 (1.42,21.86)  7.13 (1.55,32.85)  0.012          0.006     
                                                                           
Log-likelihood = -21.9491
No. of observations = 40
AIC value = 49.8982

調整あり・なしのオッズ比、Wald 型の信頼区間、Wald 検定と尤度比検定および AIC などを出力してくれます。

> logistic.display(GLM.1, simplified = TRUE)

            OR  lower95ci  upper95ci   Pr(>|Z|)
trt2 0.2110707 0.04596388  0.9692577 0.04548816
foo2 7.1291202 1.54702330 32.8529986 0.01174547

さらに、simplified = TRUE オプションを指定すると数値を出力してくれる。

補足

あとヘルプを読んでいて、gl() という関数を見つけた。

gl(3, 1, 9)
gl(3, 2, 9)
gl(3, 3, 9)

[1] 1 2 3 1 2 3 1 2 3
Levels: 1 2 3
[1] 1 1 2 2 3 3 1 1 2
Levels: 1 2 3
[1] 1 1 1 2 2 2 3 3 3
Levels: 1 2 3

第 1 引数で Levels の数を、第 2 引数で繰り返し数、第 3 引数で長さを指定できる。
返り値は "factor" のベクトルになるので、規則的な因子データを生成するときに便利かもしれない。