R语言mice多重插补后的数据处理与分析问题

背景 #

我们有含V1、V2、V3（存在缺失值）、sex、Outcome的数据集original_data，已通过mice包完成多重插补得到imp_data，以下是针对两个问题的具体解决方法：

参考代码（数据构建与插补） #

原始数据集构建 #

set.seed(123456)
original_data <- data.table(
    V1 = rnorm(20),
    V2 = rnorm(20),
    V3 = rnorm(20),
    sex = sample(c("M", "W"), 20, replace = TRUE),
    Outcome = rbinom(20, 1, 0.5)
)
original_data$V1[sample(20, 2)] <- NA
original_data$V2[sample(20, 5)] <- NA
original_data$V3[sample(20, 7)] <- NA
original_data

多重插补执行 #

imp_data <- mice::mice(original_data, m = 5) 

问题1：在插补数据集imp_data中生成衍生变量Vsum #

未插补数据中直接赋值的方式不适用于mice的插补对象，需针对每个插补后的数据集计算Vsum，推荐两种方法：

方法1：直接变换插补对象（推荐） #

使用transform()函数对imp_data进行操作，自动在所有插补数据集中添加Vsum：

imp_data_with_Vsum <- transform(imp_data, Vsum = V1 + V2 + V3)

处理后的imp_data_with_Vsum可直接用于后续的多重插补分析。

方法2：提取插补数据集后修改 #

若需要更灵活的操作，可先提取所有插补后的完整数据集，添加变量后再整合为mice对象：

# 提取全部5个插补数据集
completed_datasets <- complete(imp_data, "all")

# 遍历每个数据集添加Vsum
completed_datasets <- lapply(completed_datasets, function(df) {
  df$Vsum <- df$V1 + df$V2 + df$V3
  df
})

# 重新转换为mice对象
imp_data_with_Vsum <- mice::as.mids(completed_datasets)

问题2：计算性别特异性的比值比（OR） #

可通过分层拟合模型或添加交互项两种方式实现：

方法1：按性别分层拟合模型 #

分别对男性、女性子集拟合逻辑回归模型，再用pool()合并多重插补的结果：

# 女性组模型
model_female <- with(imp_data, glm(Outcome ~ V1, 
                                   family = binomial(link = 'logit'),
                                   data = subset(., sex == "W")))
pooled_female <- pool(model_female)
# 输出OR值（exponentiate=TRUE将对数OR转换为OR）
summary(pooled_female, exponentiate = TRUE)

# 男性组模型
model_male <- with(imp_data, glm(Outcome ~ V1, 
                                 family = binomial(link = 'logit'),
                                 data = subset(., sex == "M")))
pooled_male <- pool(model_male)
summary(pooled_male, exponentiate = TRUE)

方法2：添加性别与V1的交互项 #

拟合包含交互项的模型，一次性得到不同性别下V1的效应及交互效应：

# 拟合含交互项的模型
model_interaction <- with(imp_data, glm(Outcome ~ V1 * sex, 
                                        family = binomial(link = 'logit')))
pooled_interaction <- pool(model_interaction)
# 输出OR值
summary(pooled_interaction, exponentiate = TRUE)

结果中，V1的系数对应参考性别（默认是因子水平靠前的M）的OR，V1:sexW的系数则是女性相对于男性的效应差异，两者结合可得到女性组的OR。

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