Быстрый способ расчета доверительного интервала после изменения параметра дисперсии

Question

Я преподаю класс моделирования в R. Все студенты являются пользователями SAS, и я должен создать материалы курса, которые точно соответствуют (когда это возможно) результатам SAS. Я работаю над разделом регрессии Пуассона и пытаюсь сопоставить 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 ...). Но, очевидно, я не могу запустить confint() на summary() объекте.

Есть идеи? Бонусные баллы, если вы можете сказать мне, как изменить индекс дисперсии в модели, чтобы мне не приходилось делать это при вызове summary().

Спасибо.

Answer 1

Это интересный вопрос, и он немного глубже, чем кажется.

Самый простой потенциальный ответ - использовать family="quasipoisson" вместо пуассона:

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

Однако это победило 't позволяет вам настроить дисперсию так, как вы хотите;он конкретно изменяет дисперсию до оценки R, вычисленной в summary.glm, которая основана на хи-квадрат Пирсона (сумма квадратов остатков Пирсона), а не на отклонении, то есть

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

Вы должны знатьчто stats:::confint.glm() (который фактически использует MASS:::confint.glm) вычисляет доверительные интервалы профиля, а не доверительные интервалы Вальда (т. е. это не просто вопрос корректировки стандартных отклонений).

Если вы удовлетворены Уолдомдоверительные интервалы (которые обычно менее точны) можно взломать 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)