用于 Beta 回归的 DFFIT答案

【问题标题】：DFFITs for Beta Regression用于 Beta 回归的 DFFIT
【发布时间】：2020-12-14 21:33:24
【问题描述】：

我正在尝试计算 GLM 的 DFFITS，其中响应遵循 Beta 分布。通过使用betareg R 包。但我认为这个包不支持influence.measures()，因为使用dffits() 代码

require(betareg)
df<-data("ReadingSkills")
y<-ReadingSkills$accuracy
n<-length(y)

bfit<-betareg(accuracy ~ dyslexia + iq, data = ReadingSkills)
DFFITS<-dffits(bfit, infl=influence(bfit, do.coef = FALSE))
DFFITS

结果

if (model$rank == 0) { 中的错误：参数长度为零

我是 R 的新手。我不知道如何解决这个问题。请帮助解决这个问题，或者通过 R 代码给我一些关于如何手动计算 DFFIT 的提示。问候

【问题讨论】：

见How to make a great R reproducible example。提供您迄今为止尝试过的数据和代码。
@BappaDas 我可以在没有 infl= 参数的情况下重现它。
@jay.sf 评论后，OP已更新帖子。当我发表评论时，数据集和代码不可用。
@BappaDas 啊，我现在明白了:)

标签： r outliers

【解决方案1】：

dffits 未针对 "betareg" 对象实现，但您可以尝试手动计算。

根据this Stack Overflow Q/A我们可以写这个函数：

dffits1 <- function(x1, bres.type="response") {
  stopifnot(class(x1) %in% c("lm", "betareg"))
  sapply(1:length(x1$fitted.values), function(i) {
    x2 <- update(x1, data=x1$model[-i, ]) # leave one out
    h <- hatvalues(x1)
    nm <- rownames(x1$model[i, ])
    num_dffits <- suppressWarnings(predict(x1, x1$model[i, ]) - 
                                     predict(x2, x1$model[i, ]))
    residx <- if (class(x1) == "betareg") {
      betareg:::residuals.betareg(x2, type=bres.type)
    } else {
      x2$residuals
    }
    denom_dffits <- sqrt(c(crossprod(residx)) / x2$df.residual*h[i])
    return(num_dffits / denom_dffits)
  })
}

它适用于lm：

fit <- lm(mpg ~ hp, mtcars)
dffits1(fit)
stopifnot(all.equal(dffits1(fit), dffits(fit)))

现在让我们试试betareg：

library(betareg)
data("ReadingSkills")

bfit <- betareg(accuracy ~ dyslexia + iq, data=ReadingSkills)
dffits1(bfit)
#           1           2           3           4           5           6           7 
# -0.07590185 -0.21862047 -0.03620530  0.07349169 -0.11344968 -0.39255172 -0.25739032 
#           8           9          10          11          12          13          14 
#  0.33722706  0.16606198  0.10427684  0.11949807  0.09932991  0.11545263  0.09889406 
#          15          16          17          18          19          20          21 
#  0.21732090  0.11545263 -0.34296030  0.09850239 -0.36810187  0.09824013  0.01513643 
#          22          23          24          25          26          27          28 
#  0.18635669 -0.31192106 -0.39038732  0.09862045 -0.10859676  0.04362528 -0.28811277 
#          29          30          31          32          33          34          35 
#  0.07951977  0.02734462 -0.08419156 -0.38471945 -0.43879762  0.28583882 -0.12650591 
#          36          37          38          39          40          41          42 
# -0.12072976 -0.01701615  0.38653773 -0.06440176  0.15768684  0.05629040  0.12134228 
#          43          44 
#  0.13347935  0.19670715

看起来还不错。

注意事项：

即使这在代码中有效，您也应该检查它是否符合您的统计要求！
我在5:6 的dffits1 行中使用了suppressWarnings。 predict(bfit, ReadingSkills) 以某种方式删除了 contrasts，而 predict(bfit) 没有（实际上应该是相同的）。但是结果是相同的：all.equal(predict(bfit, ReadingSkills), predict(bfit))，因此忽略警告是安全的。

【讨论】：

非常感谢。你能再解释一件事吗？在denom_dffits 中使用了哪种类型的残差？通常使用响应残差，如果我想更改其类型，那么将在哪个步骤更改？
@JAK 我已经更改了代码，但我认为这是多余的，因为您提到了残差，您认为我应该回滚我的答案吗？ ?lm 状态 "the residuals, that is response minus fitted values." 和 ?betareg 状态 "residuals: a vector of raw residuals (observed - fitted)"。
您提供了更改预测类型的机会，但我想要残差的变化。检查这个?betareg:::residuals.betareg()
非常好。 JazakAllah 上帝保佑你。