set.seed(0)
## 2 response of 10 observations each
response <- matrix(rnorm(20), 10, 2)
## 3 covariates with 10 observations each
predictors <- matrix(rnorm(30), 10, 3)
fit <- lm(response ~ predictors)
我一直在使用以下方法为整个模型生成残差图：
plot(fitted(fit),residuals(fit))
但是，我想为每个预测变量协变量制作单独的图。我可以一次做一个：
f <- fitted(fit)
r <- residual(fit)
plot(f[,1],r[,1])
然而，这种方法的问题在于它需要对具有更多预测协变量的数据集进行泛化。有没有一种方法可以在迭代（f）和（r）的每一列时使用绘图？或者有没有一种方法plot()可以按颜色对每个协变量进行分组？
确保您使用的是标准化残差而不是原始残差
我经常看到plot(fitted(fit), residuals(fit))，但它在统计上是错误的。我们用于plot(fit)生成诊断图，因为我们需要标准化残差而不是原始残差。阅读?plot.lm更多。但是plot“mlm”的方法得到的支持很差：
plot(fit)
# Error: 'plot.mlm' is not implemented yet
为“mlm”定义“rstandard”S3 方法
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)。
现在，让我们总结一下这些想法，为“mlm”获得我们自己的“rstandard”方法：
## 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
}
注意使用.mlmin 函数名来告诉 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)分布的。
获取“mlm”的残差与拟合图
您最初的尝试：
f <- fitted(fit); r <- rstandard(fit); plot(f, r)
这不是一个坏主意，前提是可以相互识别不同模型的点。所以我们可以尝试对不同的模型使用不同的点颜色：
plot(f, r, col = as.numeric(col(f)), pch = 19)
col像,pch和的图形参数cex可以接受向量输入。我要求plot使用col = jfor r[,j] ~ f[,j], where 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])
}
脑洞大开。