2

我知道对具有多个 LHS 的线性模型的支持是有限的。但是，当可以在“mlm”对象上运行函数时，我希望结果是可靠的。使用时rstudent，会产生奇怪的结果。这是一个错误还是有其他解释？

在下面的示例中fittedA和fittedB是相同的，但在rstudent第 2 列的情况下不同。

y <- matrix(rnorm(20), 10, 2)
x <- 1:10
fittedA <- fitted(lm(y ~ x))
fittedB <- cbind(fitted(lm(y[, 1] ~ x)), fitted(lm(y[, 2] ~ x)))
rstudentA <- rstudent(lm(y ~ x))
rstudentB <- cbind(rstudent(lm(y[, 1] ~ x)), rstudent(lm(y[, 2] ~ x)))
4

2 回答 2

3

谢谢@李哲源；另请参阅https://www.r-project.org/bugs.html以了解如何报告错误……R 核心团队会更好地注意到这些错误。此外，在那里，我们可以对补丁给予更好的评价。
另外请注意，R 的源代码 - 特别是它的开发版本 - 始终可以通过 svn（“subversion”）或https://svn.r-project 上的网络浏览器获得。 org/R/主干/

cooks.distance.lm()'s 和rstandard.lm()'s 的源代码都在https://svn.r-project.org/R/trunk/src/library/stats/R/lm.influence.R ....下一次。当然，您建议的小代码更改outer()已经足够了。

非常感谢您的彻底分析和您提出的错误修复！

于 2018-09-16T16:59:06.060 回答
3

设置

set.seed(0)
y <- matrix(rnorm(20), 10, 2)
x <- 1:10
fit <- lm(y ~ x)           ## class: "mlm", "lm"
fit1 <- lm(y[, 1] ~ x)     ## class: "lm"
fit2 <- lm(y[, 2] ~ x)     ## class: "lm"

rstudent(fit)
#          [,1]        [,2]
#1   0.74417620  0.89121744
#2  -0.67506054 -0.50275275
#3   0.76297805 -0.74363941
#4   0.71164461  0.01075898
#5   0.03337192  0.03355209
#6  -1.75099724 -0.02701558
#7  -1.05594284  0.56993056
#8  -0.48486883 -0.35286612
#9  -0.23468552  0.79610101
#10  2.90701182 -0.93665406

cbind(rstudent(fit1), rstudent(fit2))
#          [,1]        [,2]
#1   0.74417620  1.90280959
#2  -0.67506054 -0.92973971
#3   0.76297805 -1.47237918
#4   0.71164461  0.01870820
#5   0.03337192  0.06042497
#6  -1.75099724 -0.04056992
#7  -1.05594284  1.02171222
#8  -0.48486883 -0.64316472
#9  -0.23468552  1.69605079
#10  2.90701182 -1.25676088

正如您所观察到的，只有第一个响应的结果由rstandard(fit).

为什么rstudent“传销”失败

问题是，没有“传销”方法rstudent。

methods(rstudent)
#[1] rstudent.glm* rstudent.lm*

当您调用rstudent(fit)时，S3 方法调度机制会找到rstudent.lm，因为inherits(fit, "lm")is TRUE。不幸的是，stats:::rstudent.lm没有对“mlm”模型进行正确的计算。

stats:::rstudent.lm
#function (model, infl = lm.influence(model, do.coef = FALSE), 
#    res = infl$wt.res, ...) 
#{
#    res <- res/(infl$sigma * sqrt(1 - infl$hat))
#    res[is.infinite(res)] <- NaN
#    res
#}

lm.influence没有给出正确sigma的“传销”。底层 C 例程C_influence仅计算sigma“lm”。如果您给出lm.influence“mlm”，则仅返回第一个响应变量的结果。

## pass in "mlm"
.Call(stats:::C_influence, fit$qr, FALSE, residuals(fit), 10 * .Machine$double.eps)$sigma
# [1] 1.3130265 1.3216357 1.3105706 1.3171621 1.3638689 1.1374385 1.2668101
# [8] 1.3416338 1.3586428 0.9180828

## pass in "lm"
.Call(stats:::C_influence, fit1$qr, FALSE, residuals(fit1), 10 * .Machine$double.eps)$sigma
# [1] 1.3130265 1.3216357 1.3105706 1.3171621 1.3638689 1.1374385 1.2668101
# [8] 1.3416338 1.3586428 0.9180828

显然，对于“传销”，sigma应该是一个矩阵。现在考虑到这个不正确sigma的，“回收规则”应用"/"在下一行的后面stats:::rstudent.lm，因为res它左边的（残差）是一个矩阵，但它右边的东西是一个向量。

res <- res / (infl$sigma * sqrt(1 - infl$hat))

实际上，计算结果只对第一个响应变量是正确的；所有其余的响应变量都会使用错误的sigma.

R核心团队需要修补一些诊断功能

请注意，文档页面中列出的几乎所有功能?influence.measures对于“mlm”都是错误的。他们应该发出警告说“传销”方法尚未实现。

lm.influnce需要在 C 级别进行修补。一旦完成，rstudent.lm将在“mlm”上正常工作。

其他功能可以在 R 级别轻松修补，例如stats:::cooks.distance.lm，stats:::rstandard. 目前（R 3.5.1）它们是：

stats:::cooks.distance.lm
#function (model, infl = lm.influence(model, do.coef = FALSE), 
#    res = weighted.residuals(model),
#    sd = sqrt(deviance(model)/df.residual(model)), 
#    hat = infl$hat, ...) 
#{
#    p <- model$rank
#    res <- ((res/(sd * (1 - hat)))^2 * hat)/p
#    res[is.infinite(res)] <- NaN
#    res
#}

stats:::rstandard.lm
#function (model, infl = lm.influence(model, do.coef = FALSE), 
#    sd = sqrt(deviance(model)/df.residual(model)), type = c("sd.1", 
#        "predictive"), ...) 
#{
#    type <- match.arg(type)
#    res <- infl$wt.res/switch(type, sd.1 = sd * sqrt(1 - infl$hat), 
#        predictive = 1 - infl$hat)
#    res[is.infinite(res)] <- NaN
#    res
#}

并且可以使用（通过使用outer）修补它们

patched_cooks.distance.lm <-
function (model, infl = lm.influence(model, do.coef = FALSE), 
    res = weighted.residuals(model),
    sd = sqrt(deviance(model)/df.residual(model)), 
    hat = infl$hat, ...) 
{
    p <- model$rank
    res <- ((res / c(outer(1 - hat, sd))) ^ 2 * hat) / p
    res[is.infinite(res)] <- NaN
    res
}

patched_rstandard.lm <- 
function (model, infl = lm.influence(model, do.coef = FALSE), 
    sd = sqrt(deviance(model)/df.residual(model)), type = c("sd.1", 
        "predictive"), ...) 
{
    type <- match.arg(type)
    res <- infl$wt.res/switch(type, sd.1 = c(outer(sqrt(1 - infl$hat), sd)), 
        predictive = 1 - infl$hat)
    res[is.infinite(res)] <- NaN
    res
}

快速测试：

oo <- cbind(cooks.distance(fit1), cooks.distance(fit2))  ## correct result
all.equal(cooks.distance(fit), oo)
#[1] "Mean relative difference: 0.8211302"
all.equal(patched_cooks.distance.lm(fit), oo)
#[1] TRUE

rr <- cbind(rstandard(fit1), rstandard(fit2))  ## correct result
all.equal(rstandard(fit), rr)
#[1] "Mean relative difference: 0.5290036"
all.equal(patched_rstandard.lm(fit), rr)
#[1] TRUE
于 2018-09-15T14:25:43.687 回答