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

Question

我正在R教一个建模课。学生都是SAS用户，我必须创建与SAS输出完全匹配（如果可能）的课程材料。我正在研究Poisson回归部分并尝试匹配PROC GENMOD，并使用“dscale”选项修改色散指数，使deviance / 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的说法......）。但显然我不能在confint()对象上运行summary()。

有任何想法吗？如果您可以告诉我如何更改模型中的色散指数，那么我不需要在summary()调用中执行此操作。

谢谢。

Answer 1

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

最简单的可能答案是使用family="quasipoisson"而不是poisson：

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

但是，这不会让你将色散调整为你想要的任何东西;它特别改变了在summary.glm中估计R计算的离散度，这是基于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)