在R中拟合数据后识别异常值——特殊环境高频次样本排查

Hey there! 针对你的问题，咱们可以把整个分析拆成验证线性关系假设和识别显著高于预期的特殊环境频次样本两个核心部分，用R就能一步步搞定，下面我给你详细拆解：

第一步：数据准备与初步探索 #

首先假设你的数据是一个包含3列的data.frame：sample_id（样本类型标识）、general（普通环境频次）、special（特殊环境频次）。先做基础的数据检查：

1. 加载必备工具包：

   library(tidyverse)   # 数据处理+可视化
   library(lmtest)      # 统计检验辅助
   library(MASS)        # 可选：处理计数数据的负二项回归
   

2. 检查数据结构与缺失值：

   # 查看数据基本信息
   str(sample_data)
   # 检查缺失值总数
   sum(is.na(sample_data))
   # 如果有缺失值，可直接删除（确保不影响结果的前提下）
   sample_data <- sample_data %>% drop_na()
   

第二步：验证「普通环境频次越高，特殊环境频次越高」的线性关系 #

咱们先通过线性回归验证你的假设，同时做模型诊断确保结果可靠：

1. 拟合线性回归模型：

   lm_model <- lm(special ~ general, data = sample_data)
   

2. 查看模型统计结果：

   summary(lm_model)
   

   重点关注这几个指标：
   • general变量的系数：如果系数为正且p值<0.05，就支持你的线性关系假设；
   • 调整后R-squared：数值越接近1，说明普通环境频次对特殊环境频次的解释力越强；
3. 模型诊断（确保线性回归假设成立）：
   • 残差正态性检验：

     # Shapiro-Wilk检验，p>0.05说明残差近似正态
     shapiro.test(residuals(lm_model))
     # 可视化QQ图，点越贴近直线越符合正态
     ggplot(data.frame(resid = residuals(lm_model)), aes(sample = resid)) +
       geom_qq() + geom_qq_line(color = "red")
     

   • 方差齐性检验：

     # Breusch-Pagan检验，p>0.05说明方差齐
     bptest(lm_model)
     # 可视化残差-拟合值图，点随机分布在0线附近就没问题
     plot(lm_model, which = 1)
     

注意：如果你的special是计数型数据（非负整数），线性回归可能不是最优选择——这时候可以用泊松回归（数据无过度离散）或负二项回归（数据过度离散）：

# 泊松回归
poisson_model <- glm(special ~ general, data = sample_data, family = poisson)
# 检查过度离散：比值>1说明存在过度离散
sum(residuals(poisson_model, type = "pearson")^2) / df.residual(poisson_model)
# 负二项回归（如果过度离散）
nb_model <- glm.nb(special ~ general, data = sample_data)

第三步：识别特殊环境频次显著高于预期的样本 #

核心思路是找到那些实际special值远高于模型预测值的样本，用标准化残差结合统计显著性来判断：

1. 计算预测值、残差及显著性指标：

   # 以线性回归为例，泊松/负二项回归只需替换model为对应模型即可
   sample_data <- sample_data %>%
     mutate(
       pred_special = predict(lm_model),          # 模型预测的special值
       residual = special - pred_special,         # 原始残差（正残差=实际>预期）
       std_resid = rstandard(lm_model),           # 标准化残差
       p_value = 2 * (1 - pnorm(abs(std_resid)))  # 双侧检验p值
     )
   

2. 筛选显著高于预期的样本（这里用p<0.05作为显著性阈值）：

   high_outliers <- sample_data %>%
     filter(residual > 0 & p_value < 0.05) %>%
     arrange(desc(std_resid))  # 按异常程度从高到低排序
   

3. 可视化标注异常样本：

   ggplot(sample_data, aes(x = general, y = special)) +
     geom_point(aes(color = ifelse(residual > 0 & p_value < 0.05, "异常高值样本", "普通样本")), alpha = 0.7) +
     geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.8) +
     scale_color_manual(values = c("异常高值样本" = "#e74c3c", "普通样本" = "#95a5a6")) +
     labs(title = "普通环境频次 vs 特殊环境频次", x = "普通环境出现频次", y = "特殊环境出现频次") +
     theme_minimal() +
     theme(legend.title = element_blank())
   

   图里的红色点就是你要找的、特殊环境频次显著高于正常水平的样本。

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