求R/Python实现：筛选满足回归经典假设的最优50样本组合

核心思路 #

直接遍历所有C(100,50)种组合完全不现实（该数值约为10^29），因此采用随机抽样+经典假设验证+模型性能排序的启发式循环方案：通过多次随机抽取50条样本，逐一验证回归经典假设，留存满足所有假设且拟合效果最优的样本组合。

需验证的经典假设及判定标准：

• 多重共线性：方差膨胀因子（VIF）<5（可根据需求放宽至10）
• 同方差性：Breusch-Pagan检验p值>0.05
• 无自相关：Durbin-Watson检验值接近2（范围1.5-2.5）或Ljung-Box检验p值>0.05
• 正态性：Shapiro-Wilk检验p值>0.05（样本量≤50适用）
• 线性性：Rainbow检验p值>0.05或残差与拟合值无显著趋势

R语言实现代码 #

假设你的数据框为df，因变量为y，自变量为其余列：

library(car)       # 用于VIF、Breusch-Pagan检验
library(lmtest)    # 用于Durbin-Watson检验
library(MASS)      # 用于Rainbow检验

# 初始化参数
total_rows <- nrow(df)
target_size <- 50
loop_times <- 1000  # 循环次数，可根据算力调整
best_r2 <- 0
best_sample_idx <- NULL
best_model <- NULL

set.seed(123)  # 固定随机种子，保证结果可复现

for (i in 1:loop_times) {
  # 随机抽取50条样本索引
  sample_idx <- sample(1:total_rows, target_size, replace = FALSE)
  sample_data <- df[sample_idx, ]
  
  # 拟合回归模型
  model <- lm(y ~ ., data = sample_data)
  
  # 逐一验证经典假设
  ## 1. 多重共线性检查
  if (any(vif(model) >= 5)) next
  
  ## 2. 同方差性检查
  if (bptest(model)$p.value <= 0.05) next
  
  ## 3. 自相关检查
  dw_stat <- dwtest(model)$statistic
  if (dw_stat < 1.5 || dw_stat > 2.5) next
  
  ## 4. 残差正态性检查
  if (shapiro.test(residuals(model))$p.value <= 0.05) next
  
  ## 5. 线性性检查
  if (rainbow(model)$p.value <= 0.05) next
  
  # 更新最优样本组合（以R²为判定指标）
  current_r2 <- summary(model)$r.squared
  if (current_r2 > best_r2) {
    best_r2 <- current_r2
    best_sample_idx <- sample_idx
    best_model <- model
  }
  
  # 打印循环进度（可选）
  if (i %% 100 == 0) {
    cat(paste("完成第", i, "次循环，当前最优R²：", round(best_r2, 4), "\n"))
  }
}

# 输出结果
cat("最优样本索引：", best_sample_idx, "\n")
summary(best_model)

Python语言实现代码 #

假设你的数据为pandas.DataFrame格式的df，因变量为y：

import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
from scipy.stats import shapiro

# 初始化参数
total_rows = df.shape[0]
target_size = 50
loop_times = 1000  # 循环次数，可根据算力调整
best_r2 = 0
best_sample_idx = None
best_model = None

np.random.seed(123)  # 固定随机种子

for i in range(loop_times):
  # 随机抽取50条样本索引
  sample_idx = np.random.choice(total_rows, target_size, replace=False)
  sample_data = df.iloc[sample_idx, :]
  
  # 拟合OLS模型（添加常数项）
  X = sm.add_constant(sample_data.drop('y', axis=1))
  y = sample_data['y']
  model = sm.OLS(y, X).fit()
  
  # 逐一验证经典假设
  ## 1. 多重共线性检查
  vif = [variance_inflation_factor(X.values, j) for j in range(X.shape[1])]
  if any(v >= 5 for v in vif):
    continue
  
  ## 2. 同方差性检查
  bp_pvalue = het_breuschpagan(model.resid, X)[1]
  if bp_pvalue <= 0.05:
    continue
  
  ## 3. 自相关检查（Ljung-Box检验替代Durbin-Watson）
  lb_pvalue = acorr_ljungbox(model.resid, lags=[1])['lb_pvalue'].values[0]
  if lb_pvalue <= 0.05:
    continue
  
  ## 4. 残差正态性检查
  if shapiro(model.resid).pvalue <= 0.05:
    continue
  
  ## 5. 线性性检查（残差与拟合值无显著趋势）
  residual_trend = sm.OLS(model.resid, sm.add_constant(model.fittedvalues)).fit()
  if residual_trend.pvalues[1] <= 0.05:
    continue
  
  # 更新最优样本组合
  current_r2 = model.rsquared
  if current_r2 > best_r2:
    best_r2 = current_r2
    best_sample_idx = sample_idx
    best_model = model
  
  # 打印循环进度（可选）
  if (i+1) % 100 == 0:
    print(f"完成第{i+1}次循环，当前最优R²：{round(best_r2, 4)}")

# 输出结果
print("最优样本索引：", best_sample_idx)
print(best_model.summary())

注意事项 #

• 若循环后未找到完全满足所有假设的样本，可适当放宽假设阈值（如VIF放宽至10，Shapiro检验p值放宽至0.01）
• 最优指标可替换为AIC、BIC或交叉验证误差，根据业务需求调整
• 算力充足时可增大loop_times，提升找到优质样本组合的概率

内容的提问来源于stack exchange，提问作者Berryl Bagas Wirawan