R语言｜使用 mice包进行多重插补与敏感性分析

`mice` 包简介

mice (Multivariate Imputation by Chained Equations) 是 R 语言中一个强大的多重插补包。它通过链式方程 (chained equations) 的方式，对缺失数据进行插补，生成多个完整的数据集，从而减少因缺失数据带来的偏差。mice 包支持多种插补方法，包括预测均值匹配 (Predictive Mean Matching, PMM)、线性回归、Logistic 回归等，可以灵活应用于各种类型的数据。

主要功能

多重插补: 生成多个完整的数据集，每个数据集都对缺失值进行了不同的插补。
链式方程: 使用链式方程的方式，对每个包含缺失值的变量，使用其他变量作为预测变量进行插补。
多种插补方法: 支持多种插补方法，包括 PMM, 线性回归, Logistic 回归等。
诊断工具: 提供多种诊断工具，用于评估插补效果，例如缺失值模式图、密度图等。
结果合并: 提供函数用于合并多个插补数据集的分析结果，例如回归系数、标准差等。

`mice` 函数详解

mice() 函数是 mice 包的核心函数，用于生成多重插补数据集。

参数说明

data: 包含不完整数据的数据框或矩阵。缺失值用 NA 表示。
m: 多重插补的数量，默认为 5。
method: 指定每个变量的插补方法。可以是单个字符串（所有变量使用相同方法），也可以是长度等于 length(blocks) 的字符串向量。如果为 NULL，则根据变量类型使用默认方法（由 defaultMethod 参数控制）。
predictorMatrix: 一个数值矩阵，指定每个目标变量的预测变量集。行数等于 length(blocks)，列数等于 ncol(data)。值为 1 表示该列变量用作行变量（目标变量）的预测变量。默认情况下，predictorMatrix 是一个对角线为 0，其余为 1 的方阵。
ignore: 一个逻辑向量，指示在创建插补模型时忽略哪些行。默认为 NULL，即包括所有具有观测值的行。
where: 一个与 data 具有相同维度的逻辑矩阵，指示应在何处创建插补。默认为 is.na(data)，即对缺失值进行插补。
blocks: 一个向量列表，包含每个块的变量名。列表元素可以命名以标识块。默认情况下，每个变量都放在自己的块中，这实际上是使用单变量模型进行完全条件规范 (FCS)。
visitSequence: 一个块名称向量，指定在 Gibbs 抽样的一次迭代期间插补块的顺序。默认值 visitSequence = "roman" 按照块在 blocks 中出现的顺序（从左到右）访问块。也可以使用以下关键字：“arabic”（从右到左）、“monotone”（按缺失数据比例从低到高排序）和 “revmonotone”（单调的逆序）。
formulas: 一个命名公式列表，或可以通过 as.formula 转换为公式的表达式。列表元素对应于块。应用列表元素的块由其名称标识，因此列表名称必须与块名称对应。
blots: 一个命名 alist 列表，可用于将参数传递给较低级别的插补函数。
post: 一个字符串向量，其长度为 ncol(data)，指定表达式作为字符串。每个字符串都在 sampler() 函数中解析和执行，以在迭代期间对插补值进行后处理。
defaultMethod: 一个长度为 4 的向量，包含以下默认插补方法：1) 数值数据，2) 具有 2 个水平的因子数据，3) 具有 > 2 个无序水平的因子数据，以及 4) 具有 > 2 个有序水平的因子数据。
maxit: 一个标量，给出迭代次数。默认为 5。
printFlag: 如果为 TRUE，则 mice 将在控制台上打印历史记录。使用 print=FALSE 进行静默计算。
seed: 一个整数，用作 set.seed() 的参数，用于偏移随机数生成器。默认值是不管随机数生成器。
data.init: 一个与 data 具有相同大小和类型的数据框，没有缺失数据，用于在迭代过程开始之前初始化插补。
…: 传递给单变量插补函数的命名参数。

`mice`包其他函数

with 函数用于对数据中的每个插补数据集执行计算。

with(data, expr, …)

本示例演示如何使用 mice 包对包含缺失值的 iris 数据集进行多重插补，并进行敏感性分析，评估插补结果的稳健性。

0. 加载必要的包

# 0. 加载必要的包
# install.packages(c("mice", "miceadds", "dplyr", "tidyr", "ggplot2", "cowplot", "naniar", "plotly"))
library(mice)
library(miceadds)
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(naniar)
library(plotly)

1. 数据准备与缺失值引入

# 1. 数据准备与缺失值引入
data <- iris
set.seed(42)
missing_indices <- sample(1:nrow(data), size = 30)
data$Sepal.Length[missing_indices] <- NA
data$Sepal.Width[sample(1:nrow(data), size = 20)] <- NA
summary(data)

2. 缺失机制检验

# 2. 缺失机制检验
mcar_test_result <- mcar_test(data)
print(mcar_test_result) # p < 0.05 为非随机缺失

3. MICE 插补

# 3. MICE 插补
md.pattern(data) # 缺失值模式

imputed_data <- mice(data, m = 5, maxit = 20, printFlag = FALSE) # 多重插补
summary(imputed_data)

plot(imputed_data) # 诊断图


densityplot(imputed_data) # 密度图

4. 插补过程追踪

# 4. 插补过程追踪
imputed_data <- mice(data, m=5, maxit=20, printFlag=TRUE, seed=42)


plot(imputed_data, c("Sepal.Length", "Sepal.Width"))

5. 回归分析与结果合并 (基准模型)

# 5. 回归分析与结果合并 (基准模型)
fit <- with(imputed_data, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_results <- pool(fit)
summary(pooled_results)

6. 模型评价 (基准模型)

# 6. 模型评价 (基准模型)
rsquared <- pool.r.squared(fit)
print(rsquared)

7. 模型评估指标 (基准模型)

# 7. 模型评估指标 (基准模型)
model_metrics <- lapply(1:imputed_data$m, function(i) {
  model <- lm(Sepal.Length ~ ., data=complete(imputed_data, i))
  data.frame(
    AIC = AIC(model),
    BIC = BIC(model),
    R2_adj = summary(model)$adj.r.squared
  )
})

bind_rows(model_metrics) %>%
  summarise_all(mean) %>%
  print()

8. 残差分析 (仅展示第一个插补数据集的残差图)

# 8. 残差分析 (仅展示第一个插补数据集的残差图)
imp_data1 <- complete(imputed_data, 1)
model1 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = imp_data1)
par(mfrow = c(2, 2))
plot(model1, which = 1:4)

9. 插补方法敏感性分析

# 9. 插补方法敏感性分析
imp_pmm <- mice(data, method = "pmm", m = 5, print = FALSE)
fit_pmm <- with(imp_pmm, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_pmm <- pool(fit_pmm)

imp_rf <- mice(data, method = "rf", m = 5, print = FALSE)
fit_rf <- with(imp_rf, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_rf <- pool(fit_rf)

results <- bind_rows(
  summary(pooled_pmm, conf.int = TRUE) %>% mutate(Method = "PMM"),
  summary(pooled_rf, conf.int = TRUE) %>% mutate(Method = "RF")
) %>%
  filter(term != "(Intercept)")

print(results)

coef_plot <- ggplot(results, aes(x = term, y = estimate, color = Method)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "变量", y = "回归系数 (95% CI)", title = "不同插补方法的系数对比") +
  theme_minimal() +
  coord_flip()

print(coef_plot)


ggplotly(coef_plot) # 动态图

10. 插补次数敏感性分析

# 10. 插补次数敏感性分析
imp_m10 <- mice(data, method = "pmm", m = 10, print = FALSE)
fit_m10 <- with(imp_m10, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_m10 <- pool(fit_m10)

results_m <- bind_rows(
  summary(pooled_results, conf.int = TRUE) %>% mutate(m = 5),
  summary(pooled_m10, conf.int = TRUE) %>% mutate(m = 10)
) %>%
  filter(term != "(Intercept)")

m_plot <- ggplot(results_m, aes(x = term, y = estimate, color = as.factor(m))) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "变量", y = "回归系数 (95% CI)", color = "插补次数 (m)",
       title = "不同插补次数的系数对比") +
  theme_minimal() +
  coord_flip()

print(m_plot)

11. 缺失比例敏感性分析

# 11. 缺失比例敏感性分析
set.seed(123)
data_30miss <- data
data_30miss$Sepal.Length[sample(1:nrow(data), 0.3 * nrow(data))] <- NA

imp_30miss <- mice(data_30miss, method = "pmm", m = 5, print = FALSE)
fit_30miss <- with(imp_30miss, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_30miss <- pool(fit_30miss)

results_miss <- bind_rows(
  summary(pooled_results, conf.int = TRUE) %>% mutate(Missing = "10%"),
  summary(pooled_30miss, conf.int = TRUE) %>% mutate(Missing = "30%")
) %>%
  filter(term != "(Intercept)")

miss_plot <- ggplot(results_miss, aes(x = term, y = estimate, color = Missing)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", color= "red") +
  labs(x = "变量", y = "回归系数 (95% CI)", title = "不同缺失比例的系数对比") +
  theme_minimal() +
  coord_flip()

print(miss_plot)

12. 模型设定敏感性分析

# 12. 模型设定敏感性分析
fit_main <- with(imputed_data, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Species))
fit_interaction <- with(imputed_data, lm(Sepal.Length ~ Sepal.Width * Petal.Length + Species))

results_model <- bind_rows(
  summary(pool(fit_main), conf.int = TRUE) %>% mutate(Model = "主效应"),
  summary(pool(fit_interaction), conf.int = TRUE) %>% mutate(Model = "交互效应")
) %>%
  filter(term %in% c("Sepal.Width", "Petal.Length", "Speciesversicolor", "Speciesvirginica"))

model_plot <- ggplot(results_model, aes(x = term, y = estimate, color = Model)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "变量", y = "回归系数 (95% CI)", title = "不同模型设定的系数对比") +
  theme_minimal() +
  coord_flip()

print(model_plot)

13. 结果可视化

# 13. 结果可视化
combined_plot <- plot_grid(
  coef_plot, m_plot, miss_plot, model_plot,
  ncol = 2, labels = "AUTO", align = "hv"
)
print(combined_plot)

14. 其他插补方法

# 14. 其他插补方法
imp_norm <- mice(data, method="norm", m=5, print=FALSE)
fit_norm <- with(imp_norm, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_norm <- pool(fit_norm)

imp_lasso <- mice(data, method="lasso.norm", m=5, print=FALSE)
fit_lasso <- with(imp_lasso, lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species))
pooled_lasso <- pool(fit_lasso)

mice 包提供了多种单变量插补方法，可以根据变量类型和数据特点选择合适的方法。以下是 mice 包中常用的插补方法及其简要说明：

数值型变量 (numeric):

pmm (Predictive Mean Matching): 预测均值匹配。从观测值中找到预测值与缺失值预测值最接近的样本，用该样本的观测值进行插补。适用于连续变量，能保持数据的分布形态。
midastouch (Weighted Predictive Mean Matching): 加权预测均值匹配。类似于 pmm，但使用加权平均的方式选择匹配样本。
mean (Unconditional Mean Imputation): 无条件平均插补。用观测值的平均值填充缺失值。简单但可能引入偏差，不推荐作为主要插补方法。
norm (Bayesian Linear Regression): 贝叶斯线性回归。使用贝叶斯线性回归模型进行插补，考虑了模型的不确定性。
norm.nob (Linear Regression Ignoring Model Error): 忽略模型误差的线性回归。使用线性回归模型进行插补，但不考虑模型误差。
norm.boot (Linear Regression Using Bootstrap): 使用 Bootstrap 的线性回归。使用 Bootstrap 方法估计线性回归模型的参数，并进行插补。
norm.predict (Linear Regression, Predicted Values): 线性回归，预测值。使用线性回归模型进行插补，直接使用预测值。
lasso.norm (Lasso Linear Regression): Lasso 线性回归。使用 Lasso 正则化的线性回归模型进行插补，可以进行变量选择。
lasso.select.norm (Lasso Select + Linear Regression): Lasso 选择 + 线性回归。先使用 Lasso 进行变量选择，然后使用线性回归模型进行插补。
quadratic (Imputation of Quadratic Terms): 二次项插补。用于插补二次项，例如 x^2。
ri (Random Indicator for Nonignorable Data): 不可忽略数据的随机指示符。用于处理 MNAR 数据，但需要谨慎使用。
2l.norm (Level-1 Normal Heteroscedastic): 水平 1 正态异方差。用于多层数据，假设水平 1 的误差项服从正态分布且具有异方差性。
2l.lmer (Level-1 Normal Homoscedastic, lmer): 水平 1 正态同方差，lmer。用于多层数据，使用 lmer 函数进行建模，假设水平 1 的误差项服从正态分布且具有同方差性。
2l.pan (Level-1 Normal Homoscedastic, pan): 水平 1 正态同方差，pan。用于多层数据，使用 pan 函数进行建模，假设水平 1 的误差项服从正态分布且具有同方差性。
2lonly.mean (Level-2 Class Mean): 二级平均值。用于多层数据，使用二级类别均值进行插补。
2lonly.norm (Level-2 Class Normal): 二级正态。用于多层数据，假设二级类别服从正态分布。
2lonly.pmm (Level-2 Class Predictive Mean Matching): 二级预测均值匹配。用于多层数据，使用二级类别预测均值匹配进行插补。

二元变量 (binary):

logreg (Logistic Regression): Logistic 回归。使用 Logistic 回归模型进行插补。
logreg.boot (Logistic Regression with Bootstrap): Logistic 回归 with Bootstrap。使用 Bootstrap 方法估计 Logistic 回归模型的参数，并进行插补。
lasso.logreg (Lasso Logistic Regression): Lasso Logistic 回归。使用 Lasso 正则化的 Logistic 回归模型进行插补，可以进行变量选择。
lasso.select.logreg (Lasso Select + Logistic Regression): Lasso 选择 + Logistic 回归。先使用 Lasso 进行变量选择，然后使用 Logistic 回归模型进行插补。
2l.bin (Level-1 Logistic, glmer): 一级 Logistic，glmer。用于多层数据，使用 glmer 函数进行建模，假设水平 1 的误差项服从 Logistic 分布。

有序分类变量 (ordered):

polr (Proportional Odds Model): 比例优势模型。使用比例优势模型进行插补。

无序分类变量 (unordered):

polyreg (Polytomous Logistic Regression): 多分类 Logistic 回归。使用多分类 Logistic 回归模型进行插补。
lda (Linear Discriminant Analysis): 线性判别分析。使用线性判别分析模型进行插补。

其他

sample (Random Sample from Observed Values): 从观测值中随机抽样。简单但可能破坏变量间的关系，不推荐作为主要插补方法。
cart (Classification and Regression Trees): 分类和回归树。使用 CART 模型进行插补，可以处理各种类型的变量。
rf (Random Forest Imputations): 随机森林插补。使用随机森林模型进行插补，可以处理各种类型的变量，且通常效果较好。

原创文章（本站视频密码：66668888），作者：xujunzju，如若转载，请注明出处：https://zyicu.cn/?p=20089

R语言｜使用 mice包进行多重插补与敏感性分析

`mice` 包简介

主要功能

`mice` 函数详解

参数说明

`mice`包其他函数

0. 加载必要的包

1. 数据准备与缺失值引入

2. 缺失机制检验

3. MICE 插补

4. 插补过程追踪

5. 回归分析与结果合并 (基准模型)

6. 模型评价 (基准模型)

7. 模型评估指标 (基准模型)

8. 残差分析 (仅展示第一个插补数据集的残差图)

9. 插补方法敏感性分析

11. 缺失比例敏感性分析

12. 模型设定敏感性分析

13. 结果可视化

14. 其他插补方法

数值型变量 (numeric):

二元变量 (binary):

有序分类变量 (ordered):

无序分类变量 (unordered):

其他

发表回复

邮箱：

xujunzju@gmail.com

公众号：

xujunzju6174

R语言｜使用 mice包进行多重插补与敏感性分析

mice 包简介

主要功能

mice 函数详解

参数说明

mice包其他函数

0. 加载必要的包

1. 数据准备与缺失值引入

2. 缺失机制检验

3. MICE 插补

4. 插补过程追踪

5. 回归分析与结果合并 (基准模型)

6. 模型评价 (基准模型)

7. 模型评估指标 (基准模型)

8. 残差分析 (仅展示第一个插补数据集的残差图)

9. 插补方法敏感性分析

11. 缺失比例敏感性分析

12. 模型设定敏感性分析

13. 结果可视化

14. 其他插补方法

数值型变量 (numeric):

二元变量 (binary):

有序分类变量 (ordered):

无序分类变量 (unordered):

其他

相关推荐

如何调整视频中的音量大小

R语言:逆概率处理加权IPTW

Pandas数据挖掘分析神器Bamboolib

发表回复

邮箱：

xujunzju@gmail.com

公众号：

xujunzju6174

`mice` 包简介

`mice` 函数详解

`mice`包其他函数