Cara: Interval prediksi untuk regresi linier melalui bootstrap

Saya mengalami kesulitan untuk memahami cara menggunakan bootstrap untuk menghitung interval prediksi untuk model regresi linier. Adakah yang bisa menjelaskan prosedur langkah demi langkah? Saya mencari melalui Google tetapi tidak ada yang benar-benar masuk akal bagi saya.

Saya mengerti bagaimana menggunakan bootstrap untuk menghitung interval kepercayaan untuk parameter model.

Interval kepercayaan memperhitungkan ketidakpastian estimasi. Interval prediksi menambah ketidakpastian mendasar ini. R's predict.lmakan memberi Anda interval prediksi untuk model linier. Dari sana, yang harus Anda lakukan adalah menjalankannya berulang kali pada sampel bootstrap.

n <- 100
n.bs <- 30

dat <- data.frame( x<-runif(n), y=x+runif(n) )
plot(y~x,data=dat)


regressAndPredict <- function( dat ) {
  model <- lm( y~x, data=dat )
  predict( model, interval="prediction" )
}

regressAndPredict(dat)

replicate( n.bs, regressAndPredict(dat[ sample(seq(n),replace=TRUE) ,]) )

Hasilnya replicateadalah array 3 dimensi ( nx 3x n.bs). Dimensi panjang 3 terdiri dari nilai pas untuk setiap elemen data, dan batas bawah / atas dari interval prediksi 95%.

Metode Gary King

Tergantung pada apa yang Anda inginkan, ada metode keren oleh King, Tomz, dan Wittenberg . Ini relatif mudah diterapkan, dan menghindari masalah bootstrap untuk perkiraan tertentu (misalnya max(Y)).

Saya akan mengutip dari definisinya tentang ketidakpastian mendasar di sini, karena ini cukup bagus:

Bentuk kedua variabilitas, ketidakpastian mendasar yang diwakili oleh komponen stokastik (distribusi f) dalam Persamaan 1, hasil dari peristiwa kebetulan yang tak terhitung jumlahnya seperti cuaca atau penyakit yang dapat mempengaruhi Y tetapi tidak termasuk dalam X. Bahkan jika kami tahu nilai eksak parameter (dengan demikian menghilangkan ketidakpastian estimasi), ketidakpastian mendasar akan mencegah kami memprediksi Y tanpa kesalahan.

Bootstrap tidak mengasumsikan pengetahuan tentang bentuk distribusi induk yang mendasari dari mana sampel muncul. Estimasi parameter statistik klasik tradisional didasarkan pada asumsi normalitas. Bootstrap berurusan dengan non-normalitas dan lebih akurat dalam praktik daripada metode klasik.

Bootstrapping menggantikan daya komputasi mentah komputer untuk analisis teoretis yang ketat. Ini adalah perkiraan untuk distribusi sampling dari istilah kesalahan kumpulan data. Bootstrap meliputi: pengambilan sampel kembali data yang ditetapkan beberapa kali, menghitung rata-rata dari setiap sampel dan menemukan kesalahan standar rata-rata.

Kode "R" berikut menunjukkan konsep:

Contoh praktis ini menunjukkan manfaat bootstrap dan memperkirakan kesalahan standar. Kesalahan standar diperlukan untuk menghitung interval kepercayaan.

Biarkan kami menganggap Anda memiliki set data miring "a":

a<-rexp(395, rate=0.1)          # Create skewed data

visualisasi dari kumpulan data yang miring

plot(a,type="l")                # Scatter plot of the skewed data
boxplot(a,type="l")             # Box plot of the skewed data
hist(a)                         # Histogram plot of the skewed data

Lakukan prosedur bootstrap:

n <- length(a)                  # the number of bootstrap samples should equal the original data set
    xbarstar <- c()                 # Declare the empty set “xbarstar” variable which will be holding the mean of every bootstrap iteration
    for (i in 1:1000) {             # Perform 1000 bootstrap iteration
        boot.samp <- sample(a, n, replace=TRUE) #”Sample” generates the same number of elements as the original data set
    xbarstar[i] <- mean(boot.samp)} # “xbarstar” variable  collects 1000 averages of the original data set
    ## 
    plot(xbarstar)                  # Scatter plot of the bootstrapped data
    boxplot(xbarstar)               # Box plot of the bootstrapped data
    hist(xbarstar)                  # Histogram plot of the bootstrapped data

    meanOfMeans <- mean(xbarstar)
    standardError <- sd(xbarstar)    # the standard error is the standard deviation of the mean of means
    confidenceIntervalAboveTheMean <- meanOfMeans + 1.96 * standardError # for 2 standard deviation above the mean 
    confidenceIntervalBelowTheMean <- meanOfMeans - 1.96 * standardError # for 2 standard deviation above the mean 
    confidenceInterval <- confidenceIntervalAboveTheMean + confidenceIntervalBelowTheMean
    confidenceInterval