缺失資料構造置信區間:《Statistical Analysis with Missing Data》習題7.9
一、題目
7.9
依題意,我們用下述方法生成模擬資料:
其中
、
均服從標準正態分佈,
。
缺失資料我們按照下述方法進行生成:
若
,則
缺失的概率為
;若
,則
缺失的概率為
我們將進行如下幾種操作:
a)構造一個檢驗來檢驗刪失是否為MCAR,並用於構造的資料集。
b)通過後面三種方法來計算
的
%置信區間,並比較效果。(i)完整的資料;(ii)有缺失值的資料;(iii)用
分佈近似來構造置信區間。
二、解答
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,所以進一步思考,懷疑出現這種結果是由於樣本量 太小,故我們設定 ,再看看檢驗情況。
# 重複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%置信區間,並比較效果
後面我們將用三種方法來計算 的 置信區間,並比較效果。(i)完整的資料;(ii)有缺失值的資料;(iii)用 分佈近似來構造置信區間(同樣我們也分別針對有缺失與無缺失的資料來進行比較)。
首先我們進行相關函式的構造
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)