从0开始学习R语言--Day60--EM插补法-优快云博客

本文链接：https://blog.youkuaiyun.com/Chef_Chen/article/details/149727771

虽然我们常常在数据处理中做数据分布以及异常值的处理，但实际上对于缺失值，很多时候我们都不能简单地去删除或赋予0值，毕竟很多都是有意义的数据，只是可能因为各种原因没有在数据面板中显示，直接删除或赋予0这种忽略数据本身意义的做法，会破坏数据的属性，扭曲数据的性质。

一般来说，对于各种缺失的数据，我们都会用EM插补法来填补数据，原理是根据缺失所属的数据列，粗略估计一个数据后，计算等于这个值的概率，然后重复这个过程指到数值不再发生变化，简单来说就是根据已有的数据列参考来回归一个数据。

以下是一个例子：

# 加载必要的包
install.packages('mice',type = 'binary')
library(mice)    # 提供EM插补功能
library(mvtnorm) # 用于生成多元正态分布数据
library(norm)    # 提供EM算法实现

# 1. 生成模拟数据集
set.seed(123)
n <- 200  # 样本量
p <- 5    # 变量数

# 生成完整的多元正态分布数据
mu <- c(5, 10, 15, 20, 25)  # 均值向量
sigma <- matrix(c(    # 协方差矩阵
  4, 2, 1, 0.5, 0.1,
  2, 9, 3, 1, 0.5,
  1, 3, 16, 4, 1,
  0.5, 1, 4, 25, 5,
  0.1, 0.5, 1, 5, 36
), ncol = p)

complete_data <- rmvnorm(n, mean = mu, sigma = sigma)
colnames(complete_data) <- paste0("X", 1:p)

# 2. 人为制造缺失值 (MCAR机制)
missing_data <- complete_data
for (j in 1:p) {
  # 每个变量随机缺失20%
  missing_indices <- sample(1:n, size = n * 0.2)
  missing_data[missing_indices, j] <- NA
}

# 查看缺失模式
summary(missing_data)
md.pattern(missing_data)

# 3. 使用norm包进行EM插补
# 首先需要对数据进行预处理
s <- prelim.norm(missing_data)  # 预处理
thetahat <- em.norm(s)          # EM算法估计参数

# 获取插补后的数据集
em_imputed <- imp.norm(s, thetahat, missing_data)

# 查看插补后的完整数据集
head(em_imputed)

# 4. 使用mice包进行EM插补 (更简单的方法)
em_mice <- mice(missing_data, method = "norm", m = 1, maxit = 50)
complete_mice <- complete(em_mice)

# 比较原始数据和插补数据
par(mfrow = c(2, 3))
for (i in 1:p) {
  plot(density(complete_data[, i], na.rm = TRUE), 
       main = paste("X", i), col = "blue")
  lines(density(complete_mice[, i], na.rm = TRUE), col = "red")
  legend("topright", legend = c("Original", "Imputed"), 
         col = c("blue", "red"), lty = 1)
}

# 5. 评估插补质量
# 计算均方误差 (仅对缺失部分)
mse <- numeric(p)
for (i in 1:p) {
  missing_ind <- is.na(missing_data[, i])
  mse[i] <- mean((complete_data[missing_ind, i] - complete_mice[missing_ind, i])^2)
}
print(paste("MSE for each variable:", paste(round(mse, 3), collapse = ", ")))

输出：