含缺失值的效应量数据集线性回归分析可行性咨询

用户问题 #

我拥有14个变量（Var1-Var14）的效应量数据，每个数值代表特定处理对对应变量的效应，正值表示促进作用，负值表示抑制作用，缺失值因部分文献未涉及对应变量所致。我需要完成两项分析：

1. 开展两两线性回归，分析所有变量间的关联；
2. 以Var1为因变量、Var2-Var14为自变量，使用glmulti包筛选最优拟合模型，确定对Var1变化影响最大的变量。

以下是样本数据代码：

set.seed(123)

# Create the dataset with effect sizes and missing values

mydata <- data.frame(
  Var1 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var2 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var3 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var4 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var5 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var6 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var7 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var8 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var9 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var10 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var11 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var12 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var13 = sample(c(-20:14, NA), 64, replace = TRUE),
  Var14 = sample(c(-20:14, NA), 64, replace = TRUE)
)

# Set more than 50% missing values in each column
for (col in 1:14) {
  missing_indices <- sample(1:64, size = 32)
  mydata[missing_indices, col] <- NA
}

请问基于这类含缺失值的数据集，能否完成上述分析？谢谢！

解答 #

一、两两线性回归分析：可以完成，需先处理缺失值 #

你的数据每列缺失率超过50%，直接用线性回归默认的列表删除法（删除所有含缺失值的观测）会导致有效样本量极低，结果可靠性差。推荐两种更合理的处理方式：

1. 多重插补（首选，适合高缺失率）

通过生成多个完整数据集进行分析，再合并结果，能最大程度保留数据信息：

library(mice)
library(broom)

# 生成5个插补后的完整数据集
imputed_data <- mice(mydata, m = 5, method = "pmm", seed = 123)

# 定义两两回归函数
pairwise_reg <- function(data) {
  vars <- colnames(data)
  results <- list()
  for (i in 1:(length(vars)-1)) {
    for (j in (i+1):length(vars)) {
      fit <- lm(as.formula(paste(vars[i], "~", vars[j])), data = data)
      results[[paste(vars[i], "_vs_", vars[j])]] <- tidy(fit)
    }
  }
  return(results)
}

# 在插补数据集上运行回归并合并结果
pooled_res <- with(imputed_data, pairwise_reg(data))
pooled_summary <- lapply(pooled_res, pool)
lapply(pooled_summary, summary)

2. 成对删除

仅删除当前分析的两个变量中存在缺失值的观测，相比列表删除能保留更多样本，但不同配对的样本量可能不一致：

vars <- colnames(mydata)
pairwise_results <- list()

for (i in 1:(length(vars)-1)) {
  for (j in (i+1):length(vars)) {
    subset_data <- na.omit(mydata[, c(vars[i], vars[j])])
    if (nrow(subset_data) >= 3) { # 保证样本量满足回归要求
      fit <- lm(as.formula(paste(vars[i], "~", vars[j])), data = subset_data)
      pairwise_results[[paste(vars[i], "_vs_", vars[j])]] <- tidy(fit)
    } else {
      pairwise_results[[paste(vars[i], "_vs_", vars[j])]] <- "样本量不足，无法拟合回归"
    }
  }
}

# 查看结果
pairwise_results

二、glmulti最优模型筛选：可以完成，需先生成完整数据集 #

glmulti本身不支持带缺失值的数据，必须先处理缺失值，推荐两种方案：

1. 多重插补+合并模型结果（首选）

library(glmulti)
library(mice)

# 多重插补生成完整数据集
imputed_data <- mice(mydata, m = 5, method = "pmm", seed = 123)

# 定义glmulti分析函数
glmulti_analyze <- function(data) {
  fit <- glmulti(Var1 ~ ., data = data[, -1], 
                 level = 1, # 仅考虑主效应
                 method = "h", # 穷举法（13个自变量，计算量可接受）
                 crit = "aic", # 用AIC筛选最优模型
                 confsetsize = 5) # 保留前5个最优模型
  return(fit)
}

# 在插补数据集上运行并合并结果
imputed_glmulti <- with(imputed_data, glmulti_analyze(data))
summary(imputed_glmulti)

2. 列表删除（谨慎使用）

仅当删除后剩余样本量足够（至少大于自变量数13）时适用：

library(glmulti)

clean_data <- na.omit(mydata)

if (nrow(clean_data) >= 15) {
  fit <- glmulti(Var1 ~ ., data = clean_data[, -1],
                 level = 1,
                 method = "h",
                 crit = "aic",
                 confsetsize = 5)
  summary(fit)
} else {
  cat("列表删除后样本量不足，无法进行模型筛选")
}

关键注意事项 #

• 高缺失率下多重插补是首选方法，列表删除会导致样本量骤降，结果偏差大；
• 若缺失不是随机的（比如部分文献系统性不报告特定变量），需先检验缺失机制，必要时调整插补方法或做敏感性分析；
• 模型筛选后，可通过变量的回归系数绝对值或标准化系数判断对Var1影响最大的变量。

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