在R中手动计算Cook's Distance时for循环遇问题求助

嘿，我来帮你解决手动计算Cook's Distance时卡在for循环的问题！先理清楚咱们的目标：Cook's Distance是用来衡量单个样本对整个回归模型拟合结果的影响，你现在选择的是逐个去掉样本重新拟合的思路，咱们一步步把循环跑通。

先梳理你已有的基础工作： #

• 你已经拿到了原始模型的拟合值：fitted <- lm10$fitted.values，长度为n
• 建好了n×n的空白矩阵lev.mat，用来存去掉每个样本后的新拟合值，还填好了第一列——这步已经很棒了！

解决for循环的具体操作： #

首先得确保你的响应变量和设计矩阵能正常调用，然后咱们把循环的逻辑理清楚：遍历每个样本i，去掉它之后重新拟合模型，再把新的拟合值塞进lev.mat对应的列里。

完整代码示例（带注释）：

# 先把原始模型的响应变量和设计矩阵提出来（假设你的lm10是用lm()拟合的）
y <- lm10$model[[1]]  # 拽出响应变量
X <- model.matrix(lm10)  # 提取带截距的设计矩阵
n <- length(y)
p <- ncol(X)  # 模型参数个数（包含截距）

# 初始化存储拟合值的矩阵（每列对应去掉第i个样本后的所有拟合值）
lev.mat <- matrix(0, nrow = n, ncol = n)

# 你已经完成的第一列填充，这里给个参考写法
lm_no1 <- lm(y[-1] ~ X[-1, ] - 1)  # 去掉第1个样本拟合，X[-1,]已有截距，所以-1避免重复加
lev.mat[, 1] <- predict(lm_no1, newdata = data.frame(X))  # 对所有样本预测拟合值

# 重点：for循环处理剩下的i=2到n
for (i in 2:n) {
  # 去掉第i个样本的响应变量和设计矩阵
  y_remove_i <- y[-i]
  X_remove_i <- X[-i, ]
  
  # 重新拟合线性模型
  lm_remove_i <- lm(y_remove_i ~ X_remove_i - 1)
  
  # 对所有n个样本预测新的拟合值，填入lev.mat的第i列
  lev.mat[, i] <- predict(lm_remove_i, newdata = data.frame(X))
}

# 现在计算Cook's Distance
# 先拿原始模型的残差和MSE
resid <- lm10$residuals
mse <- sum(resid^2) / (n - p)

# 用拟合值的变化量计算每个样本的Cook's D
cooks_d <- numeric(n)
for (i in 1:n) {
  delta_fit_sum <- sum((fitted - lev.mat[, i])^2)
  cooks_d[i] <- delta_fit_sum / (p * mse)
}

# 给你个高效验证方法：不用循环拟合，直接用杠杆值计算
hat_matrix <- X %*% solve(t(X) %*% X) %*% t(X)
h_ii <- diag(hat_matrix)
cooks_d_fast <- (resid^2 / (p * mse)) * (h_ii / (1 - h_ii)^2)

# 对比两种方法的结果，应该几乎一致（浮点误差忽略不计）
all.equal(cooks_d, cooks_d_fast)

循环里容易踩的坑，帮你提前避坑： #

• 索引搞混：确保去掉第i个样本时，y[-i]和X[-i,]是去掉对应行，别把行和列搞反
• 预测数据格式：predict函数要求newdata是数据框，所以必须把设计矩阵X转成数据框，不然会报错
• 截距重复添加：因为model.matrix出来的X已经包含截距项了，重新拟合时一定要加-1，不然模型会多拟合一个截距，结果就错了

额外小提示： #

其实Cook's Distance不用逐个拟合模型也能算，用杠杆值和残差的公式（就是上面的cooks_d_fast）效率高很多，你可以用它来验证你循环计算的结果对不对，相当于双重校验～

内容的提问来源于stack exchange，提问作者KVemuri