R語言怎么使用缺失數(shù)據(jù)的Bootstrap與Jackknife方法

本篇內(nèi)容介紹了“R語言怎么使用缺失數(shù)據(jù)的Bootstrap與Jackknife方法”的有關(guān)知識，在實際案例的操作過程中，不少人都會遇到這樣的困境，接下來就讓小編帶領(lǐng)大家學(xué)習(xí)一下如何處理這些情況吧！希望大家仔細(xì)閱讀，能夠?qū)W有所成！

一、題目

下面再加入缺失的情況來繼續(xù)深入探討，同樣還是如習(xí)題1.6的構(gòu)造方式來加入缺失值，其中a=2， b = 0

我們將進(jìn)行如下幾種操作：

二、解答

a）Bootstrap與Jackknife進(jìn)行估計

首先構(gòu)建生成數(shù)據(jù)函數(shù)。

# 生成數(shù)據(jù)
# 生成數(shù)據(jù)
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的函數(shù)：

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))
}

進(jìn)行重復(fù)試驗所需的函數(shù)：

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))
}

下面重復(fù)100次實驗進(jìn)行 Y2的均值與變異系數(shù)標(biāo)準(zhǔn)差的估計：

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)))

## 變異系數(shù)
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)))

從上面可以發(fā)現(xiàn)，Bootstrap與Jackknife兩者估計結(jié)果較為相近，其中對均值標(biāo)準(zhǔn)差的估計，Jackknife的方差更小。這其實較為符合常識：Jackknife估計每次只取出一個樣本，用剩下的樣本來作為樣本整體；而Bootstrap每次都會比較隨機(jī)地重抽樣，隨機(jī)性相對較高，所以重復(fù)100次模擬實驗，導(dǎo)致其方差相對較大。

下面我們用計算公式來進(jìn)行推導(dǎo)。

b）均值與變異系數(shù)（大樣本）的標(biāo)準(zhǔn)差解析式推導(dǎo)與計算

均值

變異系數(shù)（大樣本近似）

## 變異系數(shù)
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)
}))

變異系數(shù)大樣本近似值為：0.03717648，說明前面的Bootstrap與Jackknife兩種方法估計的都較為準(zhǔn)確。

c）缺失插補(bǔ)后的Bootstrap與Jackknife

構(gòu)造線性填補(bǔ)的函數(shù)，并進(jìn)行線性填補(bǔ)。

DatImputation <- function(dat_incomp) {
  dat_imp <- dat_incomp
  lm_model = lm(y2 ~ y1, data = na.omit(dat_incomp))
  
  # 找出y2缺失對應(yīng)的那部分data
  na_ind = is.na(dat_incomp$y2)
  na_dat = dat_incomp[na_ind, ]
  
  # 將缺失數(shù)據(jù)進(jìn)行填補(bǔ)
  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的填補(bǔ)結(jié)果，很大一部分是由于數(shù)據(jù)的缺失會造成距離真實值較遠(yuǎn)。但單從兩種方法估計出來的值比較接近。

c）缺失插補(bǔ)前的Bootstrap與Jackknife

先構(gòu)建相關(guān)的函數(shù)：

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))
}

然后看看兩種方式估計出來的結(jié)果：

Bootstrap_imp(dat$dat_incomp)$sd

Jackknife_imp(dat$dat_incomp)$sd

缺失插補(bǔ)前進(jìn)行Bootstrap與Jackknife也還是有一定的誤差，標(biāo)準(zhǔn)差都相對更大，表示波動會比較大。具體表現(xiàn)情況下面我們多次重復(fù)模擬實驗，通過90%置信區(qū)間來看各個方法的優(yōu)劣。

d）比較各種方式的90%置信區(qū)間情況（重復(fù)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))
}

重復(fù)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

變異系數(shù)情況：

result_cv <- apply(sapply(1:nrep, RepSimulationCI, 'cv'), 1, sum)
names(result_cv) <- c('boots_after', 'boots_before', 'jack_after', 'jack_before')
result_cv

上面的數(shù)字越表示90%置信區(qū)間覆蓋真實值的個數(shù)，數(shù)字越大表示覆蓋的次數(shù)越多，也就說明該方法會相對更好。

填補(bǔ)之前進(jìn)行Bootstrap或Jackknife

無論是均值還是變異系數(shù)，通過模擬實驗都能看出，在填補(bǔ)之前進(jìn)行Bootstrap或Jackknife，其估計均會遠(yuǎn)優(yōu)于在填補(bǔ)之后進(jìn)行Bootstrap或Jackknife。而具體到Bootstrap或Jackknife，這兩種方法相差無幾。

填補(bǔ)之后進(jìn)行Bootstrap或Jackknife

在填補(bǔ)之后進(jìn)行Bootstrap或Jackknife，效果都會很差，其實仔細(xì)思考后也能夠理解，本身缺失了近一半的數(shù)據(jù)，然后填補(bǔ)會帶來很大的偏差，此時我們再從中抽樣，有很大可能抽出來的絕大多數(shù)都是原本填補(bǔ)的有很大偏差的樣本，這樣估計就會更為不準(zhǔn)了。

當(dāng)然，從理論上說，填補(bǔ)之前進(jìn)行Bootstrap或Jackknife是較為合理的，這樣對每個Bootstrap或Jackknife樣本，都可以用當(dāng)前的觀測值去填補(bǔ)當(dāng)前的缺失值，這樣每次填補(bǔ)可能花費(fèi)的時間將對較長，但實際卻更有效。

“R語言怎么使用缺失數(shù)據(jù)的Bootstrap與Jackknife方法”的內(nèi)容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業(yè)相關(guān)的知識可以關(guān)注億速云網(wǎng)站，小編將為大家輸出更多高質(zhì)量的實用文章！