Apakah prediksi model Random Forest memiliki interval prediksi?

Jika saya menjalankan randomForestmodel, saya kemudian dapat membuat prediksi berdasarkan model. Apakah ada cara untuk mendapatkan interval prediksi dari masing-masing prediksi sehingga saya tahu seberapa "yakin" model dari jawabannya. Jika ini mungkin, apakah hanya didasarkan pada variabilitas variabel dependen untuk seluruh model atau akankah ada interval yang lebih luas dan sempit tergantung pada pohon keputusan tertentu yang diikuti untuk prediksi tertentu?

Ini sebagian merupakan respons terhadap @Sashikanth Dareddy (karena tidak akan cocok dengan komentar) dan sebagian lagi merupakan respons terhadap pos asli.

Ingat apa itu interval prediksi, itu adalah interval atau serangkaian nilai di mana kami memperkirakan bahwa pengamatan di masa depan akan terletak. Umumnya interval prediksi memiliki 2 bagian utama yang menentukan lebarnya, bagian yang mewakili ketidakpastian tentang prediksi rata-rata (atau parameter lain) ini adalah bagian interval kepercayaan, dan bagian yang mewakili variabilitas pengamatan individu di sekitar rata-rata itu. Interval kepercayaan cukup kuat karena Teorema Limit Sentral dan dalam kasus hutan acak, bootstrap juga membantu. Tetapi interval prediksi sepenuhnya tergantung pada asumsi tentang bagaimana data didistribusikan mengingat variabel prediktor, CLT dan bootstrap tidak berpengaruh pada bagian itu.

Interval prediksi harus lebih luas di mana interval kepercayaan yang sesuai juga akan lebih luas. Hal lain yang akan mempengaruhi lebar interval prediksi adalah asumsi tentang varians yang sama atau tidak, ini harus berasal dari pengetahuan peneliti, bukan model hutan acak.

Interval prediksi tidak masuk akal untuk hasil kategoris (Anda bisa melakukan set prediksi daripada interval, tetapi sebagian besar waktu itu mungkin tidak akan sangat informatif).

Kita bisa melihat beberapa masalah di sekitar interval prediksi dengan mensimulasikan data di mana kita tahu kebenaran sebenarnya. Pertimbangkan data berikut:

set.seed(1)

x1 <- rep(0:1, each=500)
x2 <- rep(0:1, each=250, length=1000)

y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

Data khusus ini mengikuti asumsi untuk regresi linier dan cukup lurus ke depan untuk kesesuaian hutan acak. Kita tahu dari model "benar" bahwa ketika kedua prediktor adalah 0 yang berarti 10, kita juga tahu bahwa masing-masing titik mengikuti distribusi normal dengan standar deviasi 1. Ini berarti bahwa interval prediksi 95% berdasarkan pada pengetahuan sempurna untuk titik-titik ini adalah dari 8 hingga 12 (sebenarnya 8,04 hingga 11,96, tetapi pembulatan membuatnya lebih sederhana). Setiap perkiraan interval prediksi harus lebih lebar dari ini (tidak memiliki informasi yang sempurna menambah lebar untuk mengkompensasi) dan termasuk kisaran ini.

Mari kita lihat interval dari regresi:

fit1 <- lm(y ~ x1 * x2)

newdat <- expand.grid(x1=0:1, x2=0:1)

(pred.lm.ci <- predict(fit1, newdat, interval='confidence'))
#        fit       lwr      upr
# 1 10.02217  9.893664 10.15067
# 2 14.90927 14.780765 15.03778
# 3 20.02312 19.894613 20.15162
# 4 21.99885 21.870343 22.12735

(pred.lm.pi <- predict(fit1, newdat, interval='prediction'))
#        fit      lwr      upr
# 1 10.02217  7.98626 12.05808
# 2 14.90927 12.87336 16.94518
# 3 20.02312 17.98721 22.05903
# 4 21.99885 19.96294 24.03476

Kita dapat melihat ada beberapa ketidakpastian dalam estimasi rata-rata (interval kepercayaan) dan itu memberi kita interval prediksi yang lebih luas (tetapi termasuk) kisaran 8 hingga 12.

Sekarang mari kita lihat interval berdasarkan prediksi individu dari masing-masing pohon (kita harus berharap ini menjadi lebih luas karena hutan acak tidak mendapat manfaat dari asumsi (yang kita tahu benar untuk data ini) yang dilakukan oleh regresi linier):

library(randomForest)
fit2 <- randomForest(y ~ x1 + x2, ntree=1001)

pred.rf <- predict(fit2, newdat, predict.all=TRUE)

pred.rf.int <- apply(pred.rf$individual, 1, function(x) {
  c(mean(x) + c(-1, 1) * sd(x), 
  quantile(x, c(0.025, 0.975)))
})

t(pred.rf.int)
#                           2.5%    97.5%
# 1  9.785533 13.88629  9.920507 15.28662
# 2 13.017484 17.22297 12.330821 18.65796
# 3 16.764298 21.40525 14.749296 21.09071
# 4 19.494116 22.33632 18.245580 22.09904

Interval lebih lebar daripada interval prediksi regresi, tetapi mereka tidak mencakup seluruh rentang. Mereka memang memasukkan nilai-nilai sebenarnya dan karena itu mungkin sah sebagai interval kepercayaan, tetapi mereka hanya memprediksi di mana rata-rata (nilai prediksi) berada, tidak ada bagian tambahan untuk distribusi di sekitar rata-rata itu. Untuk kasus pertama di mana x1 dan x2 keduanya 0 interval tidak pergi di bawah 9,7, ini sangat berbeda dari interval prediksi yang sebenarnya turun ke 8. Jika kita menghasilkan titik data baru maka akan ada beberapa titik (lebih banyak lagi dari 5%) yang berada dalam interval benar dan regresi, tetapi jangan jatuh dalam interval hutan acak.

Untuk menghasilkan interval prediksi, Anda perlu membuat beberapa asumsi kuat tentang distribusi titik-titik individual di sekitar rata-rata yang diprediksi, maka Anda bisa mengambil prediksi dari masing-masing pohon (potongan interval kepercayaan bootstrap) kemudian menghasilkan nilai acak dari asumsi distribusi dengan pusat itu. Kuantil untuk potongan-potongan yang dihasilkan dapat membentuk interval prediksi (tapi saya masih akan mengujinya, Anda mungkin perlu mengulangi proses beberapa kali lebih banyak dan menggabungkannya).

Berikut adalah contoh melakukan ini dengan menambahkan penyimpangan normal (karena kami tahu data asli menggunakan normal) ke prediksi dengan standar deviasi berdasarkan pada perkiraan UMK dari pohon itu:

pred.rf.int2 <- sapply(1:4, function(i) {
  tmp <- pred.rf$individual[i, ] + rnorm(1001, 0, sqrt(fit2$mse))
  quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
#           2.5%    97.5%
# [1,]  7.351609 17.31065
# [2,] 10.386273 20.23700
# [3,] 13.004428 23.55154
# [4,] 16.344504 24.35970

Interval ini berisi interval berdasarkan pengetahuan yang sempurna, jadi terlihat masuk akal. Tetapi, mereka akan sangat bergantung pada asumsi yang dibuat (asumsi tersebut valid di sini karena kami menggunakan pengetahuan tentang bagaimana data disimulasikan, mereka mungkin tidak valid dalam kasus data nyata). Saya masih akan mengulangi simulasi beberapa kali untuk data yang lebih mirip data asli Anda (tetapi disimulasikan sehingga Anda tahu kebenarannya) beberapa kali sebelum sepenuhnya mempercayai metode ini.

Saya menyadari ini adalah posting lama tetapi saya telah menjalankan beberapa simulasi tentang ini dan berpikir saya akan membagikan temuan saya.

@GregSnow membuat posting yang sangat rinci tentang ini, tetapi saya percaya ketika menghitung interval menggunakan prediksi dari pohon individu dia sedang melihat yang hanya interval prediksi 70%. Kita perlu melihat untuk mendapatkan interval prediksi 95%.[ μ + 1,96 ∗ σ , μ - 1,96 ∗ σ $[ \mu + \sigma, \mu - \sigma]$$[\mu+1.96*\sigma,\mu-1.96*\sigma]$

Membuat perubahan ini ke kode @GregSnow, kami mendapatkan hasil berikut

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

library(randomForest)
fit2 <- randomForest(y~x1+x2)
pred.rf <- predict(fit2, newdat, predict.all=TRUE)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))

pred.rf.int
                          2.5%    97.5%
1  7.826896 16.05521  9.915482 15.31431
2 11.010662 19.35793 12.298995 18.64296
3 14.296697 23.61657 14.749248 21.11239
4 18.000229 23.73539 18.237448 22.10331

Sekarang, membandingkan ini dengan interval yang dihasilkan dengan menambahkan deviasi normal ke prediksi dengan deviasi standar seperti yang disarankan oleh MSE seperti @GregSnow.

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
          2.5%    97.5%
[1,]  7.486895 17.21144
[2,] 10.551811 20.50633
[3,] 12.959318 23.46027
[4,] 16.444967 24.57601

Interval dari kedua pendekatan ini sekarang terlihat sangat dekat. Merencanakan interval prediksi untuk tiga pendekatan terhadap distribusi kesalahan dalam kasus ini terlihat seperti di bawah ini

• Garis hitam = interval prediksi dari regresi linier,
• Garis merah = Interval hutan acak yang dihitung berdasarkan Prediksi Individual,
• Garis biru = Interval hutan acak dihitung dengan menambahkan deviasi normal ke prediksi

Sekarang, mari kita jalankan kembali simulasi tetapi kali ini meningkatkan varian istilah kesalahan. Jika perhitungan interval prediksi kami baik, kami harus berakhir dengan interval yang lebih luas daripada yang kami dapatkan di atas.

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000,mean=0,sd=5)

fit1 <- lm(y~x1+x2)
newdat <- expand.grid(x1=0:1,x2=0:1)
predict(fit1,newdata=newdat,interval = "prediction")
      fit       lwr      upr
1 10.75006  0.503170 20.99695
2 13.90714  3.660248 24.15403
3 19.47638  9.229490 29.72327
4 22.63346 12.386568 32.88035

set.seed(1)
fit2 <- randomForest(y~x1+x2,localImp=T)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))
pred.rf.int
                          2.5%    97.5%
1  7.889934 15.53642  9.564565 15.47893
2 10.616744 18.78837 11.965325 18.51922
3 15.024598 23.67563 14.724964 21.43195
4 17.967246 23.88760 17.858866 22.54337

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
         2.5%    97.5%
[1,] 1.291450 22.89231
[2,] 4.193414 25.93963
[3,] 7.428309 30.07291
[4,] 9.938158 31.63777

Sekarang, ini memperjelas bahwa penghitungan interval prediksi dengan pendekatan kedua jauh lebih akurat dan menghasilkan hasil yang cukup dekat dengan interval prediksi dari regresi linier.

Mengambil asumsi normalitas, ada cara lain yang lebih mudah untuk menghitung interval prediksi dari hutan acak. Dari masing-masing pohon kami memiliki nilai prediksi ( ) dan juga kesalahan kuadrat rata-rata ( ). Jadi prediksi dari masing-masing pohon dapat dianggap sebagai . Dengan menggunakan properti distribusi normal prediksi kami dari hutan acak akan memiliki distribusi . Menerapkan ini pada contoh yang kita diskusikan di atas, kita mendapatkan hasil di bawah ini M S E i N ( μ i , R M S E i ) N ( ∑ μ i / n , ∑ R M S E i / n $\mu_i$$MSE_i$$N(\mu_i,RMSE_i)$$N(\sum \mu_i/n, \sum RMSE_i/n)$

mean.rf <- pred.rf$aggregate
sd.rf <- mean(sqrt(fit2$mse))
pred.rf.int3 <- cbind(mean.rf - 1.96*sd.rf, mean.rf + 1.96*sd.rf)
pred.rf.int3
1  1.332711 22.09364
2  4.322090 25.08302
3  8.969650 29.73058
4 10.546957 31.30789

Penghitungan ini sangat baik dengan interval model linier dan juga pendekatan yang disarankan @GregSnow. Tetapi perhatikan bahwa asumsi mendasar dalam semua metode yang kami diskusikan adalah bahwa kesalahan mengikuti distribusi Normal.

Jika Anda menggunakan R, Anda dapat dengan mudah menghasilkan interval prediksi untuk prediksi regresi hutan acak: Cukup gunakan paket quantregForest(tersedia di CRAN ) dan baca makalah oleh N. Meinshausen tentang bagaimana kuantil kondisional dapat disimpulkan dengan hutan regresi kuantil dan bagaimana mereka dapat digunakan untuk membangun interval prediksi. Sangat informatif bahkan jika Anda tidak bekerja dengan R!

Ini mudah diselesaikan dengan randomForest.

Pertama biarkan saya berurusan dengan tugas regresi (dengan asumsi hutan Anda memiliki 1.000 pohon). Dalam predictfungsi ini, Anda memiliki opsi untuk mengembalikan hasil dari setiap pohon. Ini berarti bahwa Anda akan menerima output 1000 kolom. Kita dapat mengambil rata-rata 1000 kolom untuk setiap baris - ini adalah output RF reguler yang akan dihasilkan dengan cara apa pun. Sekarang untuk mendapatkan interval prediksi katakanlah +/- 2 std. penyimpangan yang perlu Anda lakukan adalah, untuk setiap baris, dari 1000 nilai, hitung +/- 2 std. penyimpangan dan membuat ini batas atas dan bawah Anda pada prediksi Anda.

Kedua, dalam kasus klasifikasi, ingatlah bahwa setiap pohon menghasilkan 1 atau 0 (secara default) dan jumlah dari semua 1000 pohon yang dibagi dengan 1000 memberikan kemungkinan kelas (dalam kasus klasifikasi biner). Untuk menempatkan interval prediksi pada probabilitas, Anda perlu memodifikasi min. opsi nodesize (lihat dokumentasi RandomForest untuk nama persis opsi itu) setelah Anda menetapkan nilai >> 1 maka pohon individu akan menampilkan angka antara 1 dan 0. Sekarang, dari sini Anda dapat mengulangi proses yang sama seperti yang dijelaskan di atas untuk tugas regresi.

Saya harap itu masuk akal.

Saya sudah mencoba beberapa opsi (ini semua WIP):

1. Saya benar-benar membuat variabel dependen masalah klasifikasi dengan hasil sebagai rentang, bukan nilai tunggal. Hasil yang saya dapatkan buruk, dibandingkan dengan menggunakan nilai polos. Saya menyerah dengan pendekatan ini.

2. Saya kemudian mengonversinya menjadi beberapa masalah klasifikasi, yang masing-masing merupakan batas bawah untuk kisaran (hasil model adalah apakah akan melintasi batas bawah atau tidak) dan kemudian menjalankan semua model (~ 20), dan kemudian menggabungkan hasilnya untuk mendapatkan jawaban akhir sebagai rentang. Ini berfungsi lebih baik dari 1 di atas tetapi tidak sebagus yang saya butuhkan. Saya masih bekerja untuk meningkatkan pendekatan ini.

Saya menggunakan perkiraan OOB dan tinggalkan satu untuk memutuskan seberapa baik / buruk model saya.

Masalah membangun interval prediksi untuk prediksi hutan acak telah dibahas dalam makalah berikut:

Zhang, Haozhe, Joshua Zimmerman, Dan Nettleton, dan Daniel J. Nordman. "Interval Prediksi Hutan Acak." The American Statistician, 2019.

Paket R "rfinterval" adalah implementasinya yang tersedia di CRAN.

Instalasi

Untuk menginstal rfinterval paket R :

#install.packages("devtools")
#devtools::install_github(repo="haozhestat/rfinterval")
install.packages("rfinterval")
library(rfinterval)
?rfinterval

Pemakaian

Mulai cepat:

train_data <- sim_data(n = 1000, p = 10)
test_data <- sim_data(n = 1000, p = 10)

output <- rfinterval(y~., train_data = train_data, test_data = test_data,
                     method = c("oob", "split-conformal", "quantreg"),
                     symmetry = TRUE,alpha = 0.1)

### print the marginal coverage of OOB prediction interval
mean(output$oob_interval$lo < test_data$y & output$oob_interval$up > test_data$y)

### print the marginal coverage of Split-conformal prediction interval
mean(output$sc_interval$lo < test_data$y & output$sc_interval$up > test_data$y)

### print the marginal coverage of Quantile regression forest prediction interval
mean(output$quantreg_interval$lo < test_data$y & output$quantreg_interval$up > test_data$y)

Contoh data:

oob_interval <- rfinterval(pm2.5 ~ .,
                            train_data = BeijingPM25[1:1000, ],
                            test_data = BeijingPM25[1001:2000, ],
                            method = "oob",
                            symmetry = TRUE,
                            alpha = 0.1)
str(oob_interval)