Получите стандартизованные остатки и Остаточные данные по сравнению с другими. Построенный график для объекта mlm из `lm ()`

Убедитесь, что вы используете стандартизированные остатки, а не необработанные остатки

Я часто вижу plot(fitted(fit), residuals(fit)), но это статистически неверно. Мы используем plot(fit) для создания диагностической диаграммы, потому что нам нужны стандартизированные остатки, а не необработанные. Прочтите ?plot.lm, чтобы узнать больше. Но plot метод "mlm" плохо поддерживается:

plot(fit)
# Error: 'plot.mlm' is not implemented yet

Определите "rstandard" S3-метод для "mlm"

plot.mlm не поддерживается по многим причинам, одна из которых - отсутствие rstandard.mlm. Для классов «lm» и «glm» существует общий метод S3 «rstandard» для получения стандартизованных остатков:

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

Нет поддержки "млм". Так что сначала восполним этот пробел.

Получить стандартизованные остатки несложно. Пусть hii будет диагоналями матрицы шляпы, точечная оценочная стандартная ошибка для остатков будет sqrt(1 - hii) * sigma, где sigma = sqrt(RSS / df.residual) - оцененная стандартная остаточная ошибка. RSS - остаточная сумма квадратов; df.residual - остаточная степень свободы.

hii можно вычислить из матричного коэффициента Q QR-факторизации матрицы модели: rowSums(Q ^ 2). Для «mlm» существует только одно QR-разложение, поскольку матрица модели одинакова для всех ответов, поэтому нам нужно вычислить hii только один раз.

У разных ответов разные sigma, но они элегантные colSums(residuals(fit) ^ 2) / df.residual(fit).

Теперь давайте завершим эти идеи, чтобы получить наш собственный "rstandard" метод для "mlm":

## define our own "rstandard" method for "mlm" class
rstandard.mlm <- function (model) {
  Q <- with(model, qr.qy(qr, diag(1, nrow = nrow(qr$qr), ncol = qr$rank)))  ## Q matrix
  hii <- rowSums(Q ^ 2)  ## diagonal of hat matrix QQ'
  RSS <- colSums(model$residuals ^ 2)  ## residual sums of squares (for each model)
  sigma <- sqrt(RSS / model$df.residual)  ##  ## Pearson estimate of residuals (for each model)
  pointwise_sd <- outer(sqrt(1 - hii), sigma)  ## point-wise residual standard error (for each model)
  model$residuals / pointwise_sd  ## standardised residuals
  }

Обратите внимание на использование .mlm в имени функции, чтобы сообщить R, что это метод S3. После того, как мы определили эту функцию, мы можем увидеть ее в методе "rstandard":

## now there are method for "mlm"
methods(rstandard)
# [1] rstandard.glm* rstandard.lm*  rstandard.mlm

Чтобы вызвать эту функцию, нам не нужно явно вызывать rstandard.mlm; достаточно позвонить rstandard:

## test with your fitted model `fit`
rstandard(fit)
#          [,1]       [,2]
#1   1.56221865  2.6593505
#2  -0.98791320 -1.9344546
#3   0.06042529 -0.4858276
#4   0.18713629  2.9814135
#5   0.11277397  1.4336484
#6  -0.74289985 -2.4452868
#7   0.03690363  0.7015916
#8  -1.58940448 -1.2850961
#9   0.38504435  1.3907223
#10  1.34618139 -1.5900891

Стандартизированные остатки N(0, 1) распределяются.

Получение остатков по сравнению с подогнанный сюжет для "млм"

Ваша первая попытка:

f <- fitted(fit); r <- rstandard(fit); plot(f, r)

неплохая идея, при условии, что точки для разных моделей можно отличить друг от друга. Итак, мы можем попробовать использовать разные цвета точек для разных моделей:

plot(f, r, col = as.numeric(col(f)), pch = 19)

Графические аргументы, такие как col, pch и cex, могут принимать векторные входные данные. Я прошу plot использовать col = j вместо r[,j] ~ f[,j], где j = 1, 2,..., ncol(f). Прочтите «Спецификацию цвета» ?par, чтобы узнать, что означает col = j. pch = 19 указывает plot рисовать сплошные точки. Прочтите основные графические параметры для различных вариантов.

Наконец, вам может понадобиться легенда. Ты можешь сделать

plot(f, r, col = as.numeric(col(f)), pch = 19, ylim = c(-3, 4))
legend("topleft", legend = paste0("response ", 1:ncol(f)), pch = 19,
       col = 1:ncol(f), text.col = 1:ncol(f))

Чтобы оставить место для поля легенды, мы немного расширим ylim. Поскольку стандартизированные остатки равны N(0,1), ylim = c(-3, 3) - хороший диапазон. Если мы хотим разместить поле легенды вверху слева, мы расширяем ylim до c(-3, 4). Вы можете еще больше настроить легенду с помощью ncol, title и т. Д.

Сколько у вас ответов?

Если у вас есть не более нескольких ответов, приведенное выше предложение работает нормально. Если у вас их много, предлагается вывести их на отдельный участок. Цикл for, как вы выяснили, неплох, за исключением того, что вам нужно разделить область построения на разные подзаголовки, возможно, используя par(mfrow = c(?, ?)). Также установите внутреннее поле mar и внешнее поле oma, если вы воспользуетесь этим подходом. Вы можете прочитать Как создать более красивый график для моих категориальных данных временных рядов в матрице? в качестве одного из примеров этого.

Если у вас еще больше ответов, возможно, вы захотите сочетать и то, и другое? Скажем, если у вас есть 42 ответа, вы можете сделать par(mfrow = c(2, 3)), а затем нанесите 7 ответов на каждую подфигурку. Теперь решение больше основано на мнениях.

Вот как я это решил.

for(i in 1:ncol(f)) {
    plot(f[,i],r[,i])
}

Разум взорван.