r - 如何在 R 中模拟 MAR 缺失数据？

Question

我想在 R 中模拟一些丢失的数据，但遇到了麻烦。我创建了两个变量('pre' 和 'post')，它们代表相同个体治疗前后的测量值(即配对数据)。我已经能够对完全随机丢失 (MCAR) 的数据执行此操作 - 见下文，但我无法弄清楚如何为随机丢失 (MAR) 对其进行编码。对于 MAR 缺失数据，我想根据治疗前观察创建 3 个类别，这些类别将确定有多少治疗后观察缺失。 IE。

对于前 > 25，40% 后丢失
对于前 > 21 和 ≤ 25，30% 后遗漏
对于前 ≤ 21, 20% 后遗漏

任何人都可以帮忙吗？ (我真的很感激!)
谢谢

set.seed(80122)
n <- 1000

# Simulate 1000 people with high pre-treatment (mean 28, sd 3) and normal (mean 18, sd 3) post-treatment. Correlation between paired data = 0.7.
data <- rmvnorm(n,mean=c(28,18),sigma=matrix(c(9,0.7*sqrt(81),0.7*sqrt(81),9),2,2)) # Covariance matrix

# Split into pre and post treatment and check correlation is what was specified
pre <- data[, 1]
post <- data[, 2]
cor.test(pre,post)

# Simulate MCAR
mcar <- 1 - rbinom(n, 1, 0.2) # Will create ~ 20% zero's which we'll convert to NA's
post_mcar <- post
post_mcar[mcar == 0] <- mcar[mcar==0] # Replace post data with random zero's from mcar vector
post_mcar[mcar == 0] <- NA # Change zero's to NAs

Answer 1

这是一个老问题，但我想我会尝试一下。

模拟 OP 中的假数据:

library(tidyverse)
library(mvtnorm)

# Number of data values
n <- 1000

# Simulate 1000 people with high pre-treatment (mean 28, sd 3) and normal (mean 18, sd 3) post-treatment. Correlation between paired data = 0.7.
set.seed(80122)
data <- rmvnorm(n, mean=c(28,18),
                sigma=matrix(c(9,0.7*sqrt(81),0.7*sqrt(81),9),2,2)) # Covariance matrix

转换为数据框:

data = as.data.frame(data)
names(data) = c("pre", "post")

模拟完全随机丢失 (MCAR) 数据:

data$post_mcar <- data$post

set.seed(2)
data$post_mcar[sample(1:nrow(data), 0.2*nrow(data))] = NA

模拟随机缺失 (MAR) 数据:首先，我们将创建一个分组变量， frac ，其值是我们要设置为缺失的组的分数。我们将使用 cut函数来创建这些组并设置标签值，然后我们将标签转换为数字以供以后使用:

data = data %>% 
  mutate(post_mar = post,
         frac = as.numeric(as.character(cut(pre, breaks=c(-Inf, 21, 25, Inf),
                                            labels=c(0.2,0.3,0.4)))))

现在，分组 frac并将随机选择的一部分值设置为 NA , 使用 frac确定设置为 NA 的值的分数.

set.seed(3)
data = data %>% 
  group_by(frac) %>% 
  mutate(post_mar=replace(post_mar, row_number(post_mar) %in% sample(1:n(), round(unique(frac)*n())), NA)) %>% 
  ungroup

以下是结果数据框的最后 6 行:

          pre     post post_mcar post_mar frac
995  28.63476 19.35081  19.35081 19.35081  0.4
996  32.86278 24.16119        NA       NA  0.4
997  28.25965 16.64538  16.64538 16.64538  0.4
998  24.35255 17.80365  17.80365 17.80365  0.3
999  28.12426 18.25222  18.25222       NA  0.4
1000 27.55075 14.47757  14.47757 14.47757  0.4

这是对每组中缺失值比例的检查。请注意，设置为缺失值的实际百分比可能与 frac 不同。如果请求的百分比不会导致整数行。例如，这里无法选择 8 个值中的 20%。它可以是 12.5%(1 个值)或 25%(2 个值)。

data %>% group_by(frac) %>% 
  summarise(N=n(), 
            N_missing=sum(is.na(post_mar)), 
            Frac_missing=N_missing/N)

  frac   N N_missing Frac_missing
1  0.2   8         2    0.2500000
2  0.3 138        41    0.2971014
3  0.4 854       342    0.4004684