一、題目

本題基於之前習題1.6產生關於 $(Y_1,Y_2,$

我們將分別考慮三種情況：(1) $a=0,b=0$；(2) $a=2,b=0$；(3) $a=0,b=2$。

在按照上述方法生成模擬資料後，我們將通過下述不同方法，計算與比較 $Y_2$均值與方差的估計：

a) 完整樣本情形的分析；
b) Buck’s方法。基於完整樣本，通過線性迴歸給定 $Y_1$插補 $Y_2$的條件均值（c)基於此問，這裡忽略此問）；
c) 基於正態模型的隨機迴歸插補，其中，為b)中的每個條件均值新增隨機正態偏差 $N(0,s_{22.1}^2)$；
d) Hot-deck插補，通過基於 $Y_1$分佈的四分位數對完整樣本進行歸類，進而按照這個歸類，對有缺失的變數中缺失樣本進行插補。

二、解答

首先我們還是按照第一章作業，生成資料。

# 生成資料
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)

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

1) a = 0, b = 0

為了更好的比較各種插補方法的優劣，這裡通過改變隨機種子的方式，產生10個不同的資料，並以10個不同的資料來分析比較各種方法的優劣。

a) 完整樣本情形

RepGenDat = function(a, b, seed = seed) {
  dat_all = GenerateData(a = a, b = a, seed = seed)
  
  # 均值
  y2_mean = mean(dat_all$dat_comp$y2)
  # 方差
  y2_var = var(dat_all$dat_comp$y2)
  
  return(c(y2_mean = round(y2_mean, 3), y2_var = round(y2_var, 3)))
}

(rep_dat_10 = sapply(1:10, function(x) RepGenDat(a = 0, b = 0, seed = x)))

上面通過修改隨機種子，產生了10份不同的資料（列），第一行表示 $Y_2$的均值，第二行表示其方差。

c) 基於正態模型的隨機迴歸插補

首先構造一個函式，我們對缺失資料進行填補，並對填補後 $y_2$的均值和方差進行估計。

StoRegImp = function(dat_all, seed = 2018) {
  set.seed(seed)
  
  dat_comp = dat_all$dat_comp
  dat_incomp = dat_all$dat_incomp
  lm_model = lm(y2 ~ y1, data = na.omit(dat_incomp))
  
  n = nrow(dat_incomp)
  r = nrow(na.omit(dat_incomp))
  s_22.1 = sum(lm_model$residuals ^ 2) / r
  
  # 找出y2缺失對應的那部分data
  na_ind = is.na(dat_incomp$y2)
  na_dat = dat_incomp[na_ind, ]
  
  new_value = predict(lm_model, na_dat)
  
  # 將缺失資料進行填補
  dat_incomp[na_ind, 'y2'] = new_value + rnorm(n - r, sd = sqrt(s_22.1))
  
  # y2的均值與方差
  y2_mean = mean(dat_incomp$y2)
  y2_var = var(dat_incomp$y2)
  
  return(c(y2_mean = y2_mean, y2_var = y2_var))
}

由於填補具有隨機性。通過一次填補，很難得到一個讓人信服的結果，所以我們重複100次，並計算出兩個估計的均值與標準差，以便後續多種方法的比較與分析。

RepGenDat_Imp = function(a = 0, b = 0, seed_dat = seed_dat, fun = StoRegImp, B = 100) {
  dat_all = GenerateData(a = a, b = b, seed = seed_dat)
  result_B = sapply(1:B, function(i) fun(dat_all, seed = i))
  result_mean = round(apply(result_B, 1, mean), 3)
  result_std = round(apply(result_B, 1, sd), 3)
  result = paste(result_mean, result_std, sep = ' ± ')
  return(result)
}

rep_dat_sto_reg_10 = sapply(1:10, function(x) RepGenDat_Imp(a = 0, b = 0, seed_dat = x, fun = StoRegImp))
row.names(rep_dat_sto_reg_10) = c('隨機迴歸均值估計', '隨機迴歸方差估計')
rbind(rep_dat_10, rep_dat_sto_reg_10)

比較10次資料的生成後，進行隨機迴歸填補的均值與方差估計。發現 $a = 0, b = 0$的情況，線性迴歸新增噪聲項的方法估計均值與方差均較為準確。

下面我們再來關注Hot-deck插補：

d) Hot-deck插補

這裡具體的做法是先利用 $Y_1$的四分位數將樣本歸為4類，然後我們再按照4類中完整的 $Y_2$樣本隨機插補其缺失的樣本。

HotDeckImpY2 = function(dat_incomp, group) {
  y2_group = dat_incomp[dat_incomp$quantile_group == group, ]
  ind_na = is.na(y2_group$y2)
  len_na = sum(ind_na)
  
  # 用有值的樣本插補缺失的樣本
  y2_group$y2[ind_na] = sample(y2_group$y2[!ind_na], len_na, replace = TRUE)
  return(y2_group$y2)
}

HotDeckImp = function(dat_all, seed = 2018) {
  dat_incomp = dat_all$dat_incomp
  quantile_y1 = quantile(dat_incomp$y1)
  dat_incomp$quantile_group = as.numeric(cut(dat_incomp$y1, breaks = c(quantile_y1[1] - 1, quantile_y1[2:5])))
  y2_imp = unlist(lapply(1:4, function(group) HotDeckImpY2(dat_incomp, group)))
  y2_mean = mean(y2_imp)
  y2_var = var(y2_imp)
  return(c(y2_mean = y2_mean, y2_var = y2_var))
}

# rep_dat_hotdeck_10 = sapply(1:10, function(x) RepGenDat_Imp(a = 0, b = 0, seed_dat = x, fun = HotDeckImp))
# row.names(rep_dat_hotdeck_10) = c('Hot-Deck均值估計', 'Hot-Deck方差估計')
# rbind(rep_dat_10, rep_dat_sto_reg_10, rep_dat_hotdeck_10)

其實執行上述註釋的部分可以彙總檢視結果，但是由於後面我們還需要關注 $a = 2, b = 0$以及 $a = 0, b = 2$這兩種情況，所以我們先整合成一個函式，然後再通過這個函式來看最終10次資料生成展示的結果。

SummaryAll = function(a = 0, b = 0, rep = 10, B = 100) {
  rep_dat = sapply(1:rep, function(x) RepGenDat(a = a, b = b, seed = x))
  rep_dat_sto_reg = sapply(1:rep, function(x) RepGenDat_Imp(a = a, b = b, seed_dat = x, fun = StoRegImp, B = 100))
  rep_dat_hotdeck = sapply(1:rep, function(x) RepGenDat_Imp(a = a, b = b, seed_dat = x, fun = HotDeckImp, B = 100))
  
  result_all = rbind(rep_dat, rep_dat_sto_reg, rep_dat_hotdeck)
  row.names(result_all) = c('完整資料y2均值', '完整資料y2方差', '隨機迴歸均值估計', '隨機迴歸方差估計', 'Hot-Deck均值估計', 'Hot-Deck方差估計')
  colnames(result_all) = paste0('第', 1:rep, '次模擬')
  return(result_all)
}

SummaryAll()

通過上表中兩種方法與原始完整資料之間的比較我們可以發現：當 $a = 0, b = 0$（MCAR），隨機迴歸插補均值與Hot-Deck插補均值估計兩者相差不多。但對方差的估計，Hot-Deck插補之後的波動較大，表現的沒有隨機迴歸插補方差估計的好。

2) a = 2, b = 0

SummaryAll(a = 2, b = 0, rep = 1, B = 100)

接著我們考慮 $a = 2, b = 0$（NAR）。注意這裡取rep = 1，而不是為10，是因為此種情況是NAR，如果取10，有可能對某種資料在進行Hot-deck進行插補時，某一個group全部是NA，這樣就沒有辦法進行填補，從而就會出錯。同樣的 $a = 0, b = 2$也是一種NAR的缺失，所以也只用一個數據來看。

從上面的結果可以看出，Hot-Deck方法估計的均值與方差會比迴歸後隨機化插補的均值與方差的估計要差。並且方法估計會偏小，這也是比較容易理解的。因為我們就從現有資料中選擇進行插補，未新增隨機擾動，方差會比真實情況要小一些。

3) a = 0, b = 2

SummaryAll(a = 0, b = 2, rep = 1, B = 100)

接著我們再來考慮 $a = 0, b = 2$（NAR）的情況。同上面一種情況，Hot-Deck方法估計結果，無論是均值還是方差，都要略遜色於迴歸後隨機化插補方法。注意到，這時迴歸後隨機化插補方法估計的方差也遠遠小於真實的方差。因為這種刪失完全依賴於 $Y_2$具體的值，原本的波動較大，但缺失這部分資料後，波動會變小很多，此時再用這種情況插補，及時加了擾動項，也還是遠小於真實的方差。

最後，d) Hot-Deck方法雖然在這種模擬中表現的效果差於c) 隨機化線性迴歸插補的方法，但這是由於我們的模擬資料是線性構造的原因。如果我們用非線性的方法（二次、三次、指數函式等）產生資料，則可能Hot-Deck方法會更有優勢。