rstudent() возвращает неверный результат для "mlm" (линейные модели, оснащенные несколькими LHS)
Я знаю, что поддержка линейных моделей с несколькими 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)))
2 ответа
Спасибо @李哲源; см. также https://www.r-project.org/bugs.html о том, как сообщать об ошибках... где они лучше замечены основной командой R. Кроме того, там, мы могли бы дать лучший кредит для патчей..
Также отметим, что исходный код R, в частности, его версия для разработки, всегда доступен через svn ("subversion") или через веб-браузер по адресу https://svn.r-project.org/R/trunk/
И то и другое, cooks.distance.lm()и rstandard.lm()Исходный код находится в https://svn.r-project.org/R/trunk/src/library/stats/R/lm.influence.R.... все в следующий раз. Предлагаемый вами небольшой код изменяется через outer() достаточно хороши, конечно.
Большое спасибо за тщательный анализ и хорошо предложенное исправление ошибки!
Настроить
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 терпит неудачу на "млм"
Дело в том, что не существует метода "MLM" для rstudent,
methods(rstudent)
#[1] rstudent.glm* rstudent.lm*
Когда вы звоните rstudent(fit), механизм диспетчеризации метода S3 находит rstudent.lm вместо этого, потому что inherits(fit, "lm") является 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.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 будет правильно работать на "млм".
Функции Охтера могут быть легко исправлены на уровне 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