r - 更改分散参数后计算置信区间的快速方法

Question

我正在用 R 教授建模类(class)。学生都是 SAS 用户，我必须创建与(如果可能)SAS 输出完全匹配的类(class) Material 。我正在处理泊松回归部分，并尝试匹配 PROC GENMOD，使用修改色散指数的“dscale”选项，使偏差/df==1。

这很容易做到，但我需要置信区间。我想向学生们展示如何在没有手计算的情况下做到这一点。类似于 confint_default() 或 confint()

数据

skin_cancer <- data.frame(CASES=c(1,16,30,71,102,130,133,40,4,38,
                              119,221,259,310,226,65),
                      CITY=c(rep(0,8),rep(1,8)),
                      N=c(172875, 123065,96216,92051,72159,54722,
                          32185,8328,181343,146207,121374,111353,
                          83004,55932,29007,7583),
                      agegp=c(1:8,1:8))
skin_cancer$ln_n = log(skin_cancer$N)

模型

fit <- glm(CASES ~ CITY, family="poisson", offset=ln_n, data=skin_cancer)

改变色散指数

summary(fit, dispersion= deviance(fit) / df.residual(fit)))

这让我得到了“正确的”标准错误(根据 SAS 是正确的)。但显然我不能在 summary() 对象上运行 confint()。

有什么想法吗？如果您能告诉我如何更改模型中的色散指数，那么我就不必在 summary() 调用中执行此操作，即可加分。

谢谢。

Answer 1

这是一个有趣的问题，比看起来稍微深一点。

最简单的可能答案是使用 family="quasipoisson" 而不是泊松:

fitQ <- update(fit, family="quasipoisson")
confint(fitQ)

但是，这不会让您随心所欲地调整色散；它特别改变了 R 在 summary.glm 中计算的估计值的分散度，它基于 Pearson 卡方(Pearson 残差平方和)而不是偏差，即

sum((object$weights * object$residuals^2)[object$weights > 0])/df.r

您应该知道 stats:::confint.glm()(实际上使用 MASS:::confint.glm)计算剖面置信区间而不是 Wald置信区间(即，这不仅仅是调整标准差的问题)。

如果您对 Wald 置信区间(通常不太准确)感到满意，您可以如下修改 stats::confint.default()(请注意 dispersion 标题有点误导，因为这个函数基本上假设模型的原始离散度固定为 1:如果您使用估计离​​散度的模型，这将不会按预期工作。)

confint_wald_glm <- function(object, parm, level=0.95, dispersion=NULL) {
    cf <- coef(object)
    pnames <- names(cf)
    if (missing(parm)) 
      parm <- pnames
    else if (is.numeric(parm)) 
      parm <- pnames[parm]
    a <- (1 - level)/2
    a <- c(a, 1 - a)
    pct <- stats:::format.perc(a, 3)
    fac <- qnorm(a)
    ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, 
                                                               pct))
    ses <- sqrt(diag(vcov(object)))[parm]
    if (!is.null(dispersion)) ses <- sqrt(dispersion)*ses
    ci[] <- cf[parm] + ses %o% fac
    ci
}

confint_wald_glm(fit)
confint_wald_glm(fit,dispersion=2)