嘿，我来帮你搞定移除异常值后重新运行多元线性回归的需求，同时帮你对比原结果，明确这类观测值对分析的影响～

首先得找出数据里的异常值/影响性观测，常用的统计指标比如Cook距离（衡量单个观测对整体模型的影响程度），一般默认用4/n作为阈值（n是样本量），超过这个值的观测就可以视为影响性观测。我们可以在原回归流程里加入这一步：

library(tidyverse)
library(broom) # 用来提取模型结果的工具包

# 先给变量重命名，让后续代码更清晰
df <- data.frame(replicate(30, sample(0:1000, 1000, rep = TRUE)))
colnames(df)[1:5] <- paste0("predictor_", 1:5)  # 5个预测变量
colnames(df)[6:30] <- paste0("outcome_", 1:25)  # 25个因变量

# 原回归+异常值识别流程
reg_base <- df %>%
  # 把25个因变量转成长格式
  gather(outcome_name, outcome_value, -starts_with("predictor_")) %>%
  group_by(outcome_name) %>%
  nest() %>%
  mutate(
    # 拟合每个因变量的原多元回归模型
    model_original = map(data, ~ lm(outcome_value ~ ., data = .x)),
    # 提取模型诊断统计量，标记异常值
    influence_data = map(model_original, ~ augment(.x) %>% 
                           mutate(is_influential = .cooksd > 4/nrow(.x))),
    # 生成移除异常值后的数据集
    data_no_outliers = map(influence_data, ~ filter(.x, !is_influential) %>% 
                             # 去掉broom生成的辅助列，保留原始变量
                             select(outcome_value, starts_with("predictor_")))
  )

有了移除异常值后的数据集，我们就可以针对每个因变量重新跑回归，同时把原模型和新模型的结果放在一起方便对比：

# 拟合无异常值的模型，并提取关键统计量
reg_comparison <- reg_base %>%
  mutate(
    # 拟合移除异常值后的模型
    model_no_outliers = map(data_no_outliers, ~ lm(outcome_value ~ ., data = .x)),
    # 提取原模型的核心结果：R²、调整R²、系数
    stats_original = map(model_original, ~ glance(.x) %>% 
                           select(r.squared, adj.r.squared) %>% 
                           bind_cols(tidy(.x) %>% select(term, estimate))),
    # 提取无异常值模型的核心结果
    stats_no_outliers = map(model_no_outliers, ~ glance(.x) %>% 
                              select(r.squared, adj.r.squared) %>% 
                              bind_cols(tidy(.x) %>% select(term, estimate)))
  )

# 把嵌套的结果展开成表格，方便直接对比
comparison_table <- reg_comparison %>%
  select(outcome_name, stats_original, stats_no_outliers) %>%
  unnest(c(stats_original, stats_no_outliers), names_sep = "_")

现在你可以通过comparison_table来分析异常值的影响：

• 观察R²和调整R²的变化：如果某个因变量的这两个指标变化很大，说明异常值严重影响了模型的解释力
• 对比预测变量的系数：重点看系数的方向是否反转，或者系数大小的波动幅度，这能帮你判断异常值是否扭曲了变量的真实影响
• 如果需要更直观的结果，还可以统计有多少个因变量的回归结果在移除异常值后发生了实质性变化（比如系数方向改变、R²变化超过5%等）

• 异常值的阈值不是固定的：除了Cook距离的4/n，也可以用1作为阈值；或者用学生化残差（绝对值>2或3）来判断，你可以根据业务场景调整
• 别光看统计指标：如果你的数据有业务含义，建议先结合领域知识判断异常值是否合理，比如某些极端值可能是真实的业务场景，而不是错误数据
• 可以可视化Cook距离分布，直观识别异常值：

reg_base %>%
  unnest(influence_data) %>%
  ggplot(aes(x = .cooksd)) +
  geom_histogram(bins = 30, fill = "#2c3e50", alpha = 0.7) +
  facet_wrap(~outcome_name, ncol = 5) +
  geom_vline(xintercept = 4/nrow(df), color = "#e74c3c", linetype = "dashed") +
  labs(title = "Cook's Distance Distribution by Outcome Variable", 
       x = "Cook's Distance", y = "Count") +
  theme_minimal()

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