R中的问题 - GLM残差 revoScaleR::rxGlm()

Question

我可能在这里找不到答案，因为我认为revoScaleR包并没有被广泛使用。

如果我使用rxGlm()创建GLM，它可以正常工作。但是，通过rxPredict()获得的模型残差似乎只是“原始”残差，即观察值减去拟合值。各种转换版本（偏差残差、Pearson残差等）似乎不可用。

有人知道是否有办法实现这一点吗？我可以通过再次使用glm()（使用相同的公式、数据、误差结构、链接函数、权重）运行模型，并使用residuals(glm_object, type = "deviance")来获得模型的偏差残差（例如），但这很麻烦，因为glm()运行非常缓慢（大数据集，许多模型参数）。

谢谢。

编辑：包括我正在尝试遵循的文献中的此指导。

- Alan

---

你能具体说明你想要哪种残差类型（默认为deviance，还有pearson、working、response和partial）吗？我认为revoScaleR并没有提供这些选项，但可以编写代码来计算残差。 - broti

我认为理想情况下，我希望首选标准化偏差残差。谢谢。 - Alan

你的意思是指这篇文章中所定义的内容吗：https://stats.stackexchange.com/a/99723/240805？即“标准化”或“内部学生化”残差？ - broti

也许...我担心我的数学有点生疏。我编辑了我的原始帖子，包括了一个描述我正在尝试创建的东西的 .pdf 文件的截图...那是同一件事吗？ - Alan

谢谢提供截图，我开始看明白你的意思了。为了完全理解，能否说明一下书中phi、zeta和omega分别代表什么？此外，您的数据中是否涉及二元结果？在这种情况下残差会更加复杂... - broti

Answer 1

根据您的问题，很难完全理解RevoScaleR包在残差方面提供的内容以及您需要哪些残差。此外，在残差术语方面存在相当多的混淆，例如这里和这里提到的那样。

以下是一些可能会帮助您的要点/观察结果。

在线性回归中，“原始”残差与“偏差”残差相同

至少从使用glm运行玩具回归并预测结果的角度来看是这样的：

df <- mtcars
modl <- glm(formula = mpg ~ wt + qsec + am, data = mtcars)
y_hat <- predict(modl)

接下来，计算“原始”残差（预测的结果减去实际结果）以及偏差残差：

y <- as.vector(df[["mpg"]])
res_raw <- y - y_hat
res_dev <- residuals(modl, type = "deviance")

这两者是相同的:

identical(res_raw, res_dev)
[1] TRUE

我猜一旦涉及到二元结果等，就会更加复杂。

计算标准差残差的公式

使用 rstandard 方法，从 glm 中计算标准化差异残差。

res_std <- rstandard(modl)

查看getAnywhere(rstandard.glm)可以了解如何通过偏差残差手动计算标准化残差：

function (model, infl = influence(model, do.coef = FALSE), type = c("deviance", 
    "pearson"), ...) 
{
    type <- match.arg(type)
    res <- switch(type, pearson = infl$pear.res, infl$dev.res)
    res <- res/sqrt(summary(model)$dispersion * (1 - infl$hat)) # this is the key line
    res[is.infinite(res)] <- NaN
    res
}

在我的例子中，您需要手动计算标准化残差，方法是运行res/sqrt(summary(modl)$dispersion * (1 - influence(modl)$hat))。因此，您需要两个东西：hat和dispersion。我假设RevoScaleR提供了离散参数。如果在RevoScaleR中没有像influence(modl)$hat这样获取帽值的函数，那么您就必须从头开始做：

X <- as.matrix(df[, c("wt", "qsec", "am")]) # Gets the X variables
X <- cbind(rep(1, nrow(df)), X) # adds column for the constant
hat <- diag(X %*% solve(t(X) %*% X) %*% t(X)) # formula for hat values

现在计算您的标准偏差残差：

res_man <- res_raw/sqrt(summary(modl)$dispersion * (1 - hat))

这些与使用rstandard得出的相同：

head(res_man)
        Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive Hornet Sportabout           Valiant 
       -0.6254171        -0.4941877        -1.4885771         0.2297471         0.7217423        -1.1790097 
head(res_std)
        Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive Hornet Sportabout           Valiant 
       -0.6254171        -0.4941877        -1.4885771         0.2297471         0.7217423        -1.1790097 

- broti