題目

解答

由於題目要求需要重複三次類似的操作，故首先載入所需要的包，構造生成資料的函式以及繪圖的函式：

library(tidyr)    # 繪圖所需
library(ggplot2)  # 繪圖所需

# 生成資料
GenerateData <- function(a = 0, b = 0, seed = 2018) {
  set.seed(seed)
  z1 <- rnorm(100)
  z2 <- rnorm(100)
  z3 <- rnorm(100)

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

  u <- a * (y1 - 1) + b * (y2 - 5) + z3
  m2 <- 1 * (u < 0)

  y2_na <- y2
  y2_na[u < 0] <- NA
  # y2_na[as.logical(m2)] <- NA

  dat_comp <- data.frame(y1 = y1, y2 = y2)
  dat_incomp <- data.frame(y1 = y1, y2 = y2_na)
  dat_incomp <- na.omit(dat_incomp)

  return(list(dat_comp = dat_comp, dat_incomp = dat_incomp))
}

# 展現缺失出具與未缺失資料的分佈情況
PlotTwoDistribution <- function(dat) {
  p1 <- dat_comp %>%
    gather(y1, y2, key = "var", value = "value") %>%
    ggplot(aes(x = value)) +
    geom_histogram(aes(fill = factor(var), y = ..density..),
                   alpha = 0.3, colour = 'black') +
    stat_density(geom = 'line', position = 'identity', size = 1.5,
                 aes(colour = factor(var))) +
    facet_wrap(~ var, ncol = 2) +
    labs(y = '直方圖與密度曲線', x = '值',
         title = '完整無缺失資料', fill = '變數') +
    theme(plot.title = element_text(hjust = 0.5)) +
    guides(color = FALSE)

  p2 <- dat_incomp %>%
    gather(y1, y2, key = "var", value = "value") %>%
    ggplot(aes(x = value)) +
    geom_histogram(aes(fill = factor(var), y = ..density..),
                   alpha = 0.3, colour = 'black') +
    stat_density(geom = 'line', position = 'identity', size = 1.5,
                 aes(colour = factor(var))) +
    facet_wrap(~ var, ncol = 2) +
    labs(y = '直方圖與密度曲線', x = '值',
         title = '有缺失資料', fill = '變數') +
    theme(plot.title = element_text(hjust = 0.5)) +
    guides(color = FALSE)

  return(list(p_comp = p1, p_incomp = p2))
}
 

下面考慮三種情況：

1. a = 0, b = 0

a) 生成資料並繪圖展示

# 生成資料並檢視資料樣式
dat <- GenerateData(a = 0, b = 0)
dat_comp <- dat$dat_comp
dat_incomp <- dat$dat_incomp

head(dat_comp)
head(dat_incomp)

# 繪圖展示
p <- PlotTwoDistribution(dat)
p$p_comp
p$p_incomp

缺失資料與未缺失資料的分佈如上圖所示。可以發現，對於完整資料與缺失資料之間的$Y_1$的分佈與$Y_2$

b) 進行t檢驗

題設條件中說的是$Y_1$的均值,所以考慮完整資料與缺失資料（這裡的缺失指的是若$Y_2$有缺失，$Y_1$也會進行相應地缺失處理）

t.test(dat_comp$y1, dat_incomp$y1)

這裡進行t檢驗（其實不是非常嚴謹，因為不一定滿足正態假設），比較缺失與否$Y_{}$

2. a = 2, b = 0

a) 生成資料並繪圖展示

# 生成資料並檢視資料樣式
dat <- GenerateData(a = 2, b = 0)
dat_comp <- dat$dat_comp
dat_incomp <- dat$dat_incomp

head(dat_comp)
head(dat_incomp)

# 繪圖展示
p <- PlotTwoDistribution(dat)
p$p_comp
p$p_incomp

缺失資料與未缺失資料的分佈如上圖所示。可以發現，兩個資料的期望以及分佈（無論$Y_1$還是$Y_2$），整體都有一定差異。在採用$a = 2, b = 0$這種構造時，從構造的公式可以看出，$Y_2$中樣本的缺失情況與$Y_1$有關，所以這種缺失機制是：MAR。

b) 進行t檢驗

t.test(dat_comp$y1, dat_incomp$y1)

這裡進行t檢驗（其實不是非常嚴謹，因為不一定滿足正態假設），比較缺失與否$Y_1$的均值，這裡p-value = $2.398 \times 10^{-5}$，p-value非常小，說明不是MCAR，但有可能是NMAR, MAR這兩種情況。NMAR自不必提，有可能為MAR是因為，雖然是$Y_2$缺失，但其如果為MAR是有可能與$Y_1$有關的，這樣就會出現對$Y_1$進行t檢驗為顯著的情況。

3. a = 0, b = 2

a) 生成資料並繪圖展示

# 生成資料並檢視資料樣式
dat <- GenerateData(a = 0, b = 2)
dat_comp <- dat$dat_comp
dat_incomp <- dat$dat_incomp

head(dat_comp)
head(dat_incomp)

# 繪圖展示
p <- PlotTwoDistribution(dat)
p$p_comp
p$p_incomp

缺失資料與未缺失資料的分佈如上圖所示。可以發現與上一種情況一樣，兩個資料的期望以及分佈（無論$Y_1$還是$Y_2$），整體都有一定差異。在採用$a = 0, b = 2$這種構造時，從構造的公式可以看出，$Y_2$中樣本的缺失情況與$Y_2$本身有關，所以這種缺失機制是：NMAR。

b) 進行t檢驗

t.test(dat_comp$y1, dat_incomp$y1)

這裡進行t檢驗（其實不是非常嚴謹，因為不一定滿足正態假設），比較缺失與否$Y_1$的均值，這裡p-value = $3.012 \times 10^{-5}$，p-value同樣非常小，說明不是MCAR，但有可能是NMAR, MAR這兩種情況。