模拟 (Y1, Y2) 上的数据，其中 Y2 有缺失值答案

【问题标题】：Simulating Data on (Y1, Y2) where Y2 has missing values模拟 (Y1, Y2) 上的数据，其中 Y2 有缺失值
【发布时间】：2020-10-22 19:58:17
【问题描述】：

考虑一个二变量（Y₁，Y₂）问题，每个变量定义如下：

Y₁ = 1 + Z₁，完全观察到Y₁
Y₂ = 5 + 2*(Z₁) + Z₂，Y₂为如果 2*(Y₁ − 1) + Z₃ 则缺失
Z₁、Z₂ 和 Z₃ 遵循独立的标准正态分布。

我们将如何在 (Y₁, Y₂) 上模拟大小为 500 的（完整）数据集？这是我在下面写的：

    n <- 500
    y <- rnorm(n)

我们将如何模拟相应的观察数据集（通过施加缺失在 Y₂)？我不知道该去哪里回答这个问题。

    n <- 500
    z1 <- rnorm(n)
    z2 <- rnorm(n)
    z3 <- rnorm(n)

    y1 <- 1 + z1
    y2 <- 5 + 2*z1 + z2

显示完整（最初模拟的）和观察到的（施加缺失后）数据的 Y₂ 的边际分布。

【问题讨论】：

你几乎拥有它。您只需要将满足条件的 Y2 的值替换为 NA。

标签： r simulation distribution missing-data

【解决方案1】：

除了对 @jay.sf 的精彩解释之外，另一种显示分布的方法是在新变量中构建缺失数据机制并比较 y2 和 y2_missing：

library(ggplot2)
library(dplyr)
library(tidyr)
set.seed(123)
#Data
n <- 500
#Random vars
z1 <- rnorm(n)
z2 <- rnorm(n)
z3 <- rnorm(n)
#Design Y1 and Y2
y1 <- 1+z1
y2 = 5 + 2*(z1) + z2
#For missing
y2_missing <- y2
#Set missing
index <- which(((2*(y1-1))+z3)<0)
y2_missing[index]<-NA
#Complete dataset
df <- data.frame(y1,y2,y2_missing)
#Plot distributions
df %>% select(-y1) %>%
  pivot_longer(everything()) %>%
  ggplot(aes(x=value,fill=name))+
  geom_density(alpha=0.5)+
  ggtitle('Distribution for y2 and y2_missing')+
  labs(fill='Variable')+
  theme_bw()

输出：

【讨论】：

【解决方案2】：

您可能希望在数据模拟中包含一个误差项，因此应在等式中包含另一个均值为 0 的向量，再次使用 rnorm(n)。

seed <- sample(1:1e3, 1)

set.seed(635)  ## for sake of reproducibility

n <- 500
z1 <- rnorm(n)
z2 <- rnorm(n)

要获取缺失值，您可以对向量进行采样并将其设置为NA。

y2 <- 5 + 2*z1 + z2 + rnorm(n)  ## add error term independent of the `z`s

pct.mis <- .1  ## percentage missings
y2[sample(length(y2), length(y2)*pct.mis)] <- NA

## check 1: resulting missings
prop.table(table(is.na(y2)))
# FALSE  TRUE 
#   0.9   0.1 

summary(y2)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
# -2.627   3.372   5.123   4.995   6.643  13.653      50

## check 2: rounded coefficients resemble equation
fit <- lm(y2 ~ z1 + z2)
round(fit$coe)
# (Intercept)          z1          z2 
#           5           2           1 

## check 3: number of fitted values equals number of non-missing obs.
length(fit$fitted.values) / length(y2)
# [1] 0.9

【讨论】：