一、題目

7.9

依題意，我們用下述方法生成模擬資料：
$y_{i1}=z_i$

缺失資料我們按照下述方法進行生成：
若 $y_{i1}<0$，則 $y_{i2}$缺失的概率為 $0.2$；若 $y_{i1}\geq 0$，則 $y_{i2}$缺失的概率為 $0.2$

我們將進行如下幾種操作：
a）構造一個檢驗來檢驗刪失是否為MCAR，並用於構造的資料集。
b）通過後面三種方法來計算 $\mu_2$的 $95$%置信區間，並比較效果。（i）完整的資料；（ii）有缺失值的資料；（iii）用 $t$分佈近似來構造置信區間。

二、解答

a）構造一個檢驗來檢驗刪失是否為MCAR，並用於構造的資料集

首先載入依賴包：

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

構建生成資料函式。

# 生成資料
GenerateData <- function(n = 20) {
  z <- matrix(rnorm(n * 2), nrow = 2)
  
  y1 <- z[1, ]
  y2 <- z[1, ] + z[2, ]
  
  ind_na1 <- sample(which(y1 < 0), round(length(which(y1 < 0)) * 0.2))
  ind_na2 <- sample(which(y1 >= 0), round(length(which(y1 >= 0)) * 0.8))
  
  y2_na <- y2
  y2_na[c(ind_na1, ind_na2)] <- 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))
}

接著，再展現缺失出具與未缺失資料的分佈情況，來決定使用什麼檢驗。若近似正態則使用t-檢驗，若分佈非正態，則使用非引數wilcoxon-檢驗。

PlotTwoDistribution <- function(dat_incomp) {
  
  y1 <- dat_incomp$y1
  y1_na <- y1
  y1_na[is.na(dat_incomp$y2)] <- NA
  dat_plot <- data.frame(`y2未缺失` = y1, `y2缺失` = y1_na)
  
  p <- dat_plot %>%
    gather(`y2未缺失`, `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 = '兩種資料y1的分佈', fill = '變數') +
    theme(plot.title = element_text(hjust = 0.5)) +
    guides(color = FALSE)
  
  return(p)
}

# 看看兩種不同資料y1的分佈情況
dat_incomp <- GenerateData()$dat_incomp
PlotTwoDistribution(dat_incomp)

由上圖可以看出，由於樣本量較少，不服從正態分佈，所以實際的檢驗我們用非引數的檢驗會更好。不過下面我們兩種檢驗都會進行。

這裡構建進行重複test所需的函式，並同時進行t-檢驗與wilcoxon-檢驗，這裡我們只看兩種檢驗的p值。

MyTest <- function(dat = dat_incomp, method = 't.test') {
  my_test <- get(method)(na.omit(dat)$y1, dat$y1)
  p <- my_test$p.value
  return(p)
}

RepSimTest <- function(seed = 2018, n = 20) {
  set.seed(seed)
  dat_incomp <- GenerateData(n = n)$dat_incomp
  t_result <- MyTest(dat_incomp, method = 't.test')
  wilcox_result <- MyTest(dat_incomp, method = 'wilcox.test')
  return(c(t = t_result, wilcox = wilcox_result))
}

下面重複100次實驗：

# 重複100次
mat_t_wilcox <- sapply(1:100, RepSimTest)
cat('重複100次兩種檢驗，p值小於0.05的次數：\n')
rowSums(mat_t_wilcox < 0.05)
cat('重複100次兩種檢驗，p值的平均：\n')
rowMeans(mat_t_wilcox)

由兩個檢驗結果可以發現，即使我們重複了100次，也並沒有充足的理由來拒絕假設，即通過結果發現貌似是MCAR，只是p值平均都偏小。

但由於我們實際生成資料時知道，我們並不是MCAR，所以進一步思考，懷疑出現這種結果是由於樣本量 $n = 20$太小，故我們設定 $n = 100$，再看看檢驗情況。

# 重複100次
mat_t_wilcox <- sapply(1:100, RepSimTest, n = 100)
cat('重複100次兩種檢驗，p值小於0.05的次數：\n')
rowSums(mat_t_wilcox < 0.05)
cat('重複100次兩種檢驗，p值的平均：\n')
rowMeans(mat_t_wilcox)

發現前面檢驗不顯著確實是樣本量過小的問題。實際上我們的缺失是並非是MCAR，樣本量增加後我們還是可以檢驗而出的。所以今後需要注意，在樣本量過小時，很多檢驗並不靠譜，需要觀察資料本身來尋找差異。

b）通過三種方法來計算95%置信區間，並比較效果

後面我們將用三種方法來計算 $\mu_2$的 $95$置信區間，並比較效果。（i）完整的資料；（ii）有缺失值的資料；（iii）用 $t$分佈近似來構造置信區間（同樣我們也分別針對有缺失與無缺失的資料來進行比較）。

首先我們進行相關函式的構造

CalculateMissVar <- function(dat_incomp = dat_incomp) {
  y1 <- dat_incomp$y1
  y2 <- dat_incomp$y2
  dat_incomp_ <- na.omit(dat_incomp)
  y1_ <- dat_incomp_$y1
  y2_ <- dat_incomp_$y2
  
  n <- length(y2)
  r <- length(y2_)
  
  mu1 <- mean(y1)
  y1_bar <- mean(y1_)
  y2_bar <- mean(y2_)
  s11 <- var(y1_) * (r - 1) / r
  s22 <- var(y2_) * (r - 1) / r
  s12 <- cov(y1_, y2_) * (r - 1) / r
  
  sigma22.1 <- s22 - s12 ^ 2 / s11
  
  beta21.1 <- s12 / s11
  beta20.1 <- y2_bar - beta21.1 * y1_bar
  
  sigma11_hat <- var(y1) * (n - 1) / n
  sigma22_hat <- s22 + beta21.1 ^ 2 * (sigma11_hat - s11)
  
  rho <- (s12 * (s11 * s22) ^ (-1/2)) * ((sigma11_hat / s11) * (s22 / sigma22_hat)) ^ (1/2)
  
  return(sigma22.1 * (1 / r + rho ^ 2 / n / (1 - rho ^ 2) + (y1_bar - mu1) ^ 2 / r / s11))
}

CalculateCI <- function(n = 20, deleted = FALSE, t.dis = FALSE) {
  
  dat_all <- GenerateData(n = n)
  dat_comp <- dat_all$dat_comp
  dat_incomp <- dat_all$dat_incomp
  r <- length(na.omit(dat_all$dat_incomp$y2))
  
  if (t.dis == FALSE) {
    quan <- qnorm(0.975)
  } else {
    quan <- qt(0.975, df = r - 1)
  }
  
  if (deleted == FALSE) {
    y2 <- dat_all$dat_comp$y2
    ci <- c(left_value = mean(y2) - quan * sd(y2) / sqrt(n), right_value = mean(y2) + quan * sd(y2) / sqrt(n))
  } else {
    y2 <- na.omit(dat_all$dat_incomp$y2)
    dat_incomp <- dat_all$dat_incomp
    y2_sd <- sqrt(CalculateMissVar(dat_incomp))
    ci <- c(left_value = mean(y2) - quan * y2_sd, right_value = mean(y2) + quan * y2_sd)
  }
  return(ci)
}

PlotCI <- function(rep = 100, n = 20, deleted = FALSE, t.dis = FALSE) {
  
  mat_ci <- sapply(1:rep, function(x) CalculateCI(n = n, deleted = deleted, t.dis = t.dis))
  cover_times <- sum(mat_ci[1, ] * mat_ci[2, ] < 0)
  
  ind <- rep(seq(0, 1, length.out = rep), each = 2)
  value <- matrix(mat_ci, ncol = 1)
  group <- rep(1:rep, each = 2)
  in_ci <- rep(mat_ci[1, ] * mat_ci[2, ] < 0, each = 2)
  
  dat_ci <- data.frame(ind, value, group, in_ci)
  
  p <- ggplot(dat = dat_ci) +
    geom_line(aes(x = ind, y = value, group = group, color = in_ci), size = 1) +
    geom_hline(yintercept = 0, size = 1) + 
    labs(y = '置信區間', x = paste0(rep, '次'),
         title = paste0('重複', rep, '次，覆蓋次數為：', cover_times)) +
    theme_bw() + theme(panel.grid = element_blank(), 
                       axis.text = element_blank(),
                       axis.ticks = element_blank(),
                       axis.title = element_blank(),
                       panel.border = element_blank(),
                       plot.title = element_text(hjust = 0.5),
                       legend.position = "none")
  return(p)
}

###（i）無刪失 + 正態

PlotCI(200, deleted = F, t.dis = F)

（ii）刪失 + 正態

PlotCI(200, deleted = T, t.dis = F)

（iii）刪失 + t分佈

PlotCI(200, deleted = T, t.dis = T)

（iv）無刪失 + t分佈

PlotCI(200, deleted = F, t.dis = T)

從上面四種情況發現（題目要求應該是前三種，後面自己補了第四種），第一種“無刪失+正態”的情況，發現置信區間相對較準；但第二種“有刪失+正態”的情況，置信區間的覆蓋頻率在70%左右，相對來說差了一些；第三種情況“有刪失+t分佈”，相比第二種，效果又有了一定的提升，因為樣本量本身就小，並且還刪失了接近一半的值，所以用t分佈的分位數來構造置信區間效果更好；至於最後一種“無刪失+t分佈”的情況，只是“畫蛇添足”加上去的，發現效果更差，因為我們本身的例子就是用正態分佈構造的。

錯誤的方法

最後補一下之前做錯的方法：沒有用極大似然估計，針對缺失資料直接使用sd函式來計算標準差。這樣做會發現第二種與第三種效果都很差，覆蓋頻率在60%左右，達不到用極大似然估計出來的結果。

下面放上之前錯誤的第二第三種情況：

CalculateCI <- function(n = 20, deleted = FALSE, t.dis = FALSE) {
  
  dat_all <- GenerateData(n = n)
  dat_comp <- dat_all$dat_comp
  dat_incomp <- dat_all$dat_incomp
  r <- length(na.omit(dat_all$dat_incomp$y2))
  
  if (t.dis == FALSE) {
    quan <- qnorm(0.975)
  } else {
    quan <- qt(0.975, df = r - 1)
  }
  
  if (deleted == FALSE) {
    y2 <- dat_all$dat_comp$y2
    ci <- c(left_value = mean(y2) - quan * sd(y2) / sqrt(n), right_value = mean(y2) + quan * sd(y2) / sqrt(n))
  } else {
    y2 <- na.omit(dat_all$dat_incomp$y2)
    ci <- c(left_value = mean(y2) - quan * sd(y2) / sqrt(n), right_value = mean(y2) + quan * sd(y2) / sqrt(n))
  }
  return(ci)
}

（ii）刪失 + 正態

PlotCI(200, deleted = T, t.dis = F)

（iii）刪失 + t分佈

PlotCI(200, deleted = T, t.dis = T)