R語言學習筆記缺失數據的Bootstrap與Jackknife方法
一、題目
下面再加入缺失的情況來繼續深入探討,同樣還是如習題1.6的構造方式來加入缺失值,其中a=2, b = 0
我們將進行如下幾種操作:
二、解答
a)Bootstrap與Jackknife進行估計
首先構建生成數據函數。
# 生成數據 # 生成數據 GenerateData <- function(a = 0, b = 0) { y <- matrix(nrow = 3, ncol = 100) z <- matrix(rnorm(300), nrow = 3) y[1, ] <- 1 + z[1, ] y[2, ] <- 5 + 2 * z[1, ] + z[2, ] u <- a * (y[1, ] - 1) + b * (y[2, ] - 5) + z[3, ] # m2 <- 1 * (u < 0) y[3, ] <- y[2, ] y[3, u < 0] <- NA dat_comp <- data.frame(y1 = y[1, ], y2 = y[2, ]) dat_incomp <- data.frame(y1 = y[1, ], y2 = y[3, ]) # dat_incomp <- na.omit(dat_incomp) return(list(dat_comp = dat_comp, dat_incomp = dat_incomp)) }
Bootstrap與Jackknife的函數:
Bootstrap1 <- function(Y, B = 200, fun) { Y_len <- length(Y) mat_boots <- matrix(sample(Y, Y_len * B, replace = T), nrow = B, ncol = Y_len) statis_boots <- apply(mat_boots, 1, fun) boots_mean <- mean(statis_boots) boots_sd <- sd(statis_boots) return(list(mean = boots_mean, sd = boots_sd)) } Jackknife1 <- function(Y, fun) { Y_len <- length(Y) mat_jack <- sapply(1:Y_len, function(i) Y[-i]) redu_samp <- apply(mat_jack, 2, fun) jack_mean <- mean(redu_samp) jack_sd <- sqrt(((Y_len - 1) ^ 2 / Y_len) * var(redu_samp)) return(list(mean = jack_mean, sd = jack_sd)) }
進行重復試驗所需的函數:
RepSimulation <- function(seed = 2018, fun) { set.seed(seed) dat <- GenerateData() dat_comp_y2 <- dat$dat_comp$y2 boots_sd <- Bootstrap1(dat_comp_y2, B = 200, fun)$sd jack_sd <- Jackknife1(dat_comp_y2, fun)$sd return(c(boots_sd = boots_sd, jack_sd = jack_sd)) }
下面重復100次實驗進行 Y2的均值與變異系數標準差的估計:
nrep <- 100 ## 均值 fun = mean mat_boots_jack <- sapply(1:nrep, RepSimulation, fun) apply(mat_boots_jack, 1, function(x) paste(round(mean(x), 3), '±', round(sd(x), 3)))
## 變異系數 fun = function(x) sd(x) / mean(x) mat_boots_jack <- sapply(1:nrep, RepSimulation, fun) apply(mat_boots_jack, 1, function(x) paste(round(mean(x), 3), '±', round(sd(x), 3)))
從上面可以發現,Bootstrap與Jackknife兩者估計結果較為相近,其中對均值標準差的估計,Jackknife的方差更小。這其實較為符合常識:Jackknife估計每次隻取出一個樣本,用剩下的樣本來作為樣本整體;而Bootstrap每次都會比較隨機地重抽樣,隨機性相對較高,所以重復100次模擬實驗,導致其方差相對較大。
下面我們用計算公式來進行推導。
b)均值與變異系數(大樣本)的標準差解析式推導與計算
均值
變異系數(大樣本近似)
## 變異系數 sd(sapply(1:10000, function(x) { set.seed(x) dat <- GenerateData(a = 0, b = 0) sd(dat$dat_comp$y2) / mean(dat$dat_comp$y2) }))
變異系數大樣本近似值為:0.03717648,說明前面的Bootstrap與Jackknife兩種方法估計的都較為準確。
c)缺失插補後的Bootstrap與Jackknife
構造線性填補的函數,並進行線性填補。
DatImputation <- function(dat_incomp) { dat_imp <- dat_incomp lm_model = lm(y2 ~ y1, data = na.omit(dat_incomp)) # 找出y2缺失對應的那部分data na_ind = is.na(dat_incomp$y2) na_dat = dat_incomp[na_ind, ] # 將缺失數據進行填補 dat_imp[na_ind, 'y2'] = predict(lm_model, na_dat) return(dat_imp) } dat <- GenerateData(a = 2, b = 0) dat_imp <- DatImputation(dat$dat_incomp)
fun = mean Bootstrap1(dat_imp$y2, B = 200, fun)$sd
Jackknife1(dat_imp$y2, fun)$sd
fun = function(x) sd(x) / mean(x) Bootstrap1(dat_imp$y2, B = 200, fun)$sd
Jackknife1(dat_imp$y2, fun)$sd
Bootstrap與Jackknife的填補結果,很大一部分是由於數據的缺失會造成距離真實值較遠。但單從兩種方法估計出來的值比較接近。
c)缺失插補前的Bootstrap與Jackknife
先構建相關的函數:
Array2meancv <- function(j, myarray) { dat_incomp <- as.data.frame(myarray[, j, ]) names(dat_incomp) <- c('y1', 'y2') dat_imp <- DatImputation(dat_incomp) y2_mean <- mean(dat_imp$y2) y2_cv <- sd(dat_imp$y2) / y2_mean return(c(mean = y2_mean, cv = y2_cv)) } Bootstrap_imp <- function(dat_incomp, B = 200) { n <- nrow(dat_incomp) array_boots <- array(dim = c(n, B, 2)) mat_boots_ind <- matrix(sample(1:n, n * B, replace = T), nrow = B, ncol = n) array_boots[, , 1] <- sapply(1:B, function(i) dat_incomp$y1[mat_boots_ind[i, ]]) array_boots[, , 2] <- sapply(1:B, function(i) dat_incomp$y2[mat_boots_ind[i, ]]) mean_cv_imp <- sapply(1:B, Array2meancv, array_boots) boots_imp_mean <- apply(mean_cv_imp, 1, mean) boots_imp_sd <- apply(mean_cv_imp, 1, sd) return(list(mean = boots_imp_mean, sd = boots_imp_sd)) } Jackknife_imp <- function(dat_incomp) { n <- nrow(dat_incomp) array_jack <- array(dim = c(n - 1, n, 2)) array_jack[, , 1] <- sapply(1:n, function(i) dat_incomp[-i, 'y1']) array_jack[, , 2] <- sapply(1:n, function(i) dat_incomp[-i, 'y2']) mean_cv_imp <- sapply(1:n, Array2meancv, array_jack) jack_imp_mean <- apply(mean_cv_imp, 1, mean) jack_imp_sd <- apply(mean_cv_imp, 1, function(x) sqrt(((n - 1) ^ 2 / n) * var(x))) return(list(mean = jack_imp_mean, sd = jack_imp_sd)) }
然後看看兩種方式估計出來的結果:
Bootstrap_imp(dat$dat_incomp)$sd
Jackknife_imp(dat$dat_incomp)$sd
缺失插補前進行Bootstrap與Jackknife也還是有一定的誤差,標準差都相對更大,表示波動會比較大。具體表現情況下面我們多次重復模擬實驗,通過90%置信區間來看各個方法的優劣。
d)比較各種方式的90%置信區間情況(重復100次實驗)
RepSimulationCI <- function(seed = 2018, stats = 'mean') { mean_true <- 5 cv_true <- sqrt(5) / 5 myjudge <- function(x, value) { return(ifelse((x$mean - qnorm(0.95) * x$sd < value) & (x$mean + qnorm(0.95) * x$sd > value), 1, 0)) } if(stats == 'mean') { fun = mean value = mean_true } else if(stats == 'cv') { fun = function(x) sd(x) / mean(x) value = cv_true } set.seed(seed) boots_after_ind <- boots_before_ind <- jack_after_ind <- jack_before_ind <- 0 dat <- GenerateData(a = 2, b = 0) dat_incomp <- dat$dat_incomp # after imputation dat_imp <- DatImputation(dat_incomp) boots_after <- Bootstrap1(dat_imp$y2, B = 200, fun) boots_after_ind <- myjudge(boots_after, value) jack_after <- Jackknife1(dat_imp$y2, fun) jack_after_ind <- myjudge(jack_after, value) # before imputation boots_before <- Bootstrap_imp(dat_incomp) jack_before <- Jackknife_imp(dat_incomp) if(stats == 'mean') { boots_before$mean <- boots_before$mean[1] boots_before$sd <- boots_before$sd[1] jack_before$mean <- jack_before$mean[1] jack_before$sd <- jack_before$sd[1] } else if(stats == 'cv') { boots_before$mean <- boots_before$mean[2] boots_before$sd <- boots_before$sd[2] jack_before$mean <- jack_before$mean[2] jack_before$sd <- jack_before$sd[2] } boots_before_ind <- myjudge(boots_before, value) jack_before_ind <- myjudge(jack_before, value) return(c(boots_after = boots_after_ind, boots_before = boots_before_ind, jack_after = jack_after_ind, jack_before = jack_before_ind)) }
重復100次實驗,均值情況:
nrep <- 100 result_mean <- apply(sapply(1:nrep, RepSimulationCI, 'mean'), 1, sum) names(result_mean) <- c('boots_after', 'boots_before', 'jack_after', 'jack_before') result_mean
變異系數情況:
result_cv <- apply(sapply(1:nrep, RepSimulationCI, 'cv'), 1, sum) names(result_cv) <- c('boots_after', 'boots_before', 'jack_after', 'jack_before') result_cv
上面的數字越表示90%置信區間覆蓋真實值的個數,數字越大表示覆蓋的次數越多,也就說明該方法會相對更好。
填補之前進行Bootstrap或Jackknife
無論是均值還是變異系數,通過模擬實驗都能看出,在填補之前進行Bootstrap或Jackknife,其估計均會遠優於在填補之後進行Bootstrap或Jackknife。而具體到Bootstrap或Jackknife,這兩種方法相差無幾。
填補之後進行Bootstrap或Jackknife
在填補之後進行Bootstrap或Jackknife,效果都會很差,其實仔細思考後也能夠理解,本身缺失瞭近一半的數據,然後填補會帶來很大的偏差,此時我們再從中抽樣,有很大可能抽出來的絕大多數都是原本填補的有很大偏差的樣本,這樣估計就會更為不準瞭。
當然,從理論上說,填補之前進行Bootstrap或Jackknife是較為合理的,這樣對每個Bootstrap或Jackknife樣本,都可以用當前的觀測值去填補當前的缺失值,這樣每次填補可能花費的時間將對較長,但實際卻更有效。
以上就是R語言學習筆記缺失數據的Bootstrap與Jackknife方法的詳細內容,更多關於R語言學習筆記的資料請關註WalkonNet其它相關文章!
推薦閱讀:
- python 實現以相同規律打亂多組數據
- python中pandas操作apply返回多列的實現
- R語言中矩陣matrix和數據框data.frame的使用詳解
- R語言學習Rcpp基礎知識全面整理
- R語言 實現將1對多數據與1對1數據互換