基礎的すぎるから知らなかったのか 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" のベクトルになるので、規則的な因子データを生成するときに便利かもしれない。