Referensi Utama

Judul Asli
Prediction, Bootstrapping and Monte Carlo Analyses Based on Linear Mixed Models with QAPE 2.0 Package

Terjemahan Judul
Prediksi, Bootstrap, dan Analisis Monte Carlo Berdasarkan Model Campuran Linier (LMM) dengan Paket QAPE 2.0

Penulis
- Alicja Wolny-Dominiak
Department of Statistical and Mathematical Methods in Economics
- Tomasz Żądło
Department of Statistics, Econometrics and Mathematics

Publikasi
- Diterbitkan: 10 Januari 2025
- Diterima: 12 Agustus 2022
- DOI: 10.32614/RJ-2024-004
- Jurnal: The R Journal, Volume 16/1, Halaman 67–82

1. Pendahuluan

Dalam banyak bidang terapan seperti survei, pertanian, asuransi, dan kesehatan, sering kali kita ingin melakukan prediksi karakteristik populasi atau subpopulasi, seperti rata-rata pendapatan, jumlah hasil panen, atau premi asuransi.

Namun, prediksi ini tidak cukup hanya menghasilkan angka. Kita juga perlu: - Memperkirakan seberapa akurat prediksi tersebut - Menyesuaikan model dengan struktur data yang berjenjang atau berkelompok

Paper ini memperkenalkan paket qape versi 2.0 di R, yang menyediakan alat lengkap untuk:

🔁 Inti Penjelasan Konseptual:

  • LMM → menjadi fondasi statistik utama yang menangani struktur data yang memiliki efek tetap dan efek acak.

  • EBLUP / EBP / PLUG-IN → adalah berbagai metode prediksi berbasis LMM, masing-masing dengan kelebihan dan keterbatasan dalam menangani bentuk fungsi karakteristik populasi.

  • Bootstrap → digunakan sebagai cara berbasis data (data-driven) untuk mengukur seberapa besar ketidakpastian dari prediksi, terutama saat formula analitik sulit atau bias.

  • Monte Carlo → adalah metode simulasi acak berulang-ulang untuk menguji apakah prediksi dan estimasi akurasinya tetap konsisten dan andal dalam jangka panjang.

Secara keseluruhan, qape memungkinkan pengguna R untuk: - Membuat prediksi yang fleksibel - Mengukur akurasi secara menyeluruh - Menguji kekuatan metode mereka melalui simulasi berulang (Monte Carlo)

2. Ukuran Akurasi Prediksi: RMSE dan QAPE

Setelah diperoleh nilai prediksi \(\hat{\theta}\) dari metode seperti EBLUP, EBP, atau PLUG-IN, langkah berikutnya adalah mengevaluasi akurasi prediksi terhadap nilai sebenarnya \(\theta\).

2.1 Root Mean Squared Error (RMSE)

RMSE merupakan ukuran umum untuk menilai akurasi prediksi. RMSE menghitung akar dari rata-rata kuadrat selisih antara prediksi dan nilai sebenarnya.

\[ \text{RMSE}(\hat{\theta}) = \sqrt{\mathbb{E}[(\hat{\theta} - \theta)^2]} = \sqrt{\mathbb{E}(U^2)} \]

dengan: - \(\hat{\theta}\) = nilai prediksi
- \(\theta\) = nilai sebenarnya
- \(U = \hat{\theta} - \theta\) = error prediksi

Catatan:
RMSE sangat sensitif terhadap pencilan. Jika distribusi error sangat miring atau mengandung nilai ekstrem, maka RMSE bisa memberikan gambaran yang menyesatkan karena terlalu dipengaruhi oleh nilai-nilai besar.

2.2 Quantile of Absolute Prediction Errors (QAPE)

Sebagai alternatif dari RMSE, digunakan QAPE yang mengukur akurasi berdasarkan kuantil dari galat absolut.

\[ QAPE_p(\hat{\theta}) = \inf \left\{ x : P\left( \left| \hat{\theta} - \theta \right| \le x \right) \ge p \right\} \]

Artinya: - Untuk \(p = 0.75\), 75% galat absolut berada di bawah QAPE\(_{0.75}\). - QAPE lebih tahan terhadap pencilan dibanding RMSE. - QAPE memberikan informasi tentang penyebaran error pada level kuantil tertentu.

2.3 Implementasi dalam Paket qape

Paket qape menyediakan fungsi-fungsi untuk menghitung estimasi dari RMSE dan QAPE, baik: - Dengan bootstrap (parametrik, residual, double) - Maupun melalui simulasi Monte Carlo

Fungsi utama untuk ini meliputi: - bootRes() → residual bootstrap - bootPar() → parametric bootstrap - bootResMis(), bootParMis() → untuk model salah spesifikasi - mcLMMmis() → simulasi Monte Carlo

Dengan demikian, pengguna memiliki dua sudut pandang dalam mengukur kualitas prediksi: 1. Rata-rata kesalahan kuadrat (RMSE) 2. Distribusi error absolut (QAPE)

Penggunaan keduanya membuat evaluasi lebih menyeluruh, terutama pada distribusi error yang tidak normal.

3. Prediksi Berdasarkan Linear Mixed Model (LMM)

Bagian ini menunjukkan bagaimana melakukan prediksi karakteristik populasi dengan menggunakan dua metode yang berbasis LMM, yaitu: - EBLUP (Empirical Best Linear Unbiased Predictor) - PLUG-IN Predictor

Data yang digunakan adalah dataset radon, yang berisi pengukuran kadar radon (dalam log skala) di rumah-rumah pada 85 county di negara bagian Minnesota, AS.

# Set CRAN mirror supaya gak error saat install package
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Muat library dan data
library(qape)
## Warning: package 'qape' was built under R version 4.4.3
library(lme4)
## Loading required package: Matrix
library(psych)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(HLMdiag)
## Warning: package 'HLMdiag' was built under R version 4.4.3
## 
## Attaching package: 'HLMdiag'
## The following object is masked from 'package:stats':
## 
##     covratio
data("radon")
head(radon)
##   log.radon basement    uranium county county.name
## 1 0.7884574        1 -0.6890476      1      AITKIN
## 2 0.7884574        0 -0.6890476      1      AITKIN
## 3 1.0647107        0 -0.6890476      1      AITKIN
## 4 0.0000000        0 -0.6890476      1      AITKIN
## 5 1.1314021        0 -0.8473129      2       ANOKA
## 6 0.9162907        0 -0.8473129      2       ANOKA
str(radon)
## 'data.frame':    919 obs. of  5 variables:
##  $ log.radon  : num  0.788 0.788 1.065 0 1.131 ...
##  $ basement   : int  1 0 0 0 0 0 0 0 0 0 ...
##  $ uranium    : num  -0.689 -0.689 -0.689 -0.689 -0.847 ...
##  $ county     : int  1 1 1 1 2 2 2 2 2 2 ...
##  $ county.name: Factor w/ 85 levels "AITKIN","ANOKA",..: 1 1 1 1 2 2 2 2 2 2 ...

Dataset radon terdiri dari 919 observasi dan 5 variabel sebagai berikut:

3.1 Struktur Model LMM

Model ini bertujuan memprediksi kadar log-radon berdasarkan dua variabel prediktor: - uranium: kandungan uranium tanah - basement: apakah pengukuran dilakukan di basement atau lantai 1

Karena data berasal dari 85 county yang berbeda, kita gunakan model Linear Mixed Model dengan intercept dan slope acak untuk tiap county.

library(lme4)

# Membuat model LMM
radon.model <- lmer(log.radon ~ basement + uranium + (basement | county), data = radon)

# Melihat ringkasan hasil model
summary(radon.model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log.radon ~ basement + uranium + (basement | county)
##    Data: radon
## 
## REML criterion at convergence: 2128.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1161 -0.6127  0.0288  0.6379  3.5101 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  county   (Intercept) 0.01668  0.1291       
##           basement    0.12753  0.3571   0.21
##  Residual             0.56006  0.7484       
## Number of obs: 919, groups:  county, 85
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  1.46265    0.03565  41.026
## basement    -0.64239    0.08737  -7.353
## uranium      0.76801    0.08888   8.641
## 
## Correlation of Fixed Effects:
##          (Intr) basmnt
## basement -0.260       
## uranium   0.181 -0.045

1. Kriteria REML
- Nilai REML: 2128.6
- Cocok untuk membandingkan model campuran dengan fixed effects yang sama

2. Scaled Residuals
- Outlier kuat: residual -5.1 dan +3.5
- Sebagian besar residual kecil dan seimbang (1Q dan 3Q dekat nol)
- Model fit cukup baik, ada observasi ekstrem yang perlu diperhatikan

3. Random Effects (Efek Acak)
- Intercept antar county: SD ≈ 0.1291
- Slope basement antar county: SD ≈ 0.3571
- Korelasi intercept & slope basement: 0.21
- Residual (galat tidak dijelaskan): SD ≈ 0.7484

4. Fixed Effects (Efek Tetap)
- Intercept: 1.46265 (baseline saat basement = 0, uranium = 0)
- Basement: -0.64239 (rumah dengan basement → kadar radon lebih rendah)
- Uranium: 0.76801 (setiap +1 unit uranium → log.radon naik 0.768)
- Semua efek tetap signifikan (|t| > 2)

5. Korelasi Antar Efek Tetap
- Korelasi antar prediktor kecil (|r| < 0.3)
- Tidak ada masalah multikolinearitas

6. Kesimpulan Umum
- Uranium berpengaruh positif terhadap kadar radon
- Model secara umum fit dengan baik, meskipun ada outlier yang perlu perhatian

3.2 Prediksi Menggunakan EBLUP

Tujuan kita: memprediksi rata-rata log-radon di county ke-26, dengan asumsi bahwa data dari lantai 1 di county tersebut tidak tersedia.

library(qape)
library(dplyr)

# 1. Vektor populasi (log-radon)
Ypop <- radon$log.radon

# 2. Vektor penanda: data dari lantai 1 di county 26 dianggap tidak tersedia
con <- rep(1, length(Ypop))
con[radon$county == 26 & radon$basement == 1] <- 0

# 3. Ambil sampel berdasarkan penanda
YS <- Ypop[con == 1]

# 4. Matriks prediktor dari semua data
reg <- select(radon, -log.radon)

# 5. Komponen model
fixed.part <- 'basement + uranium'
random.part <- '(basement|county)'

# 6. Bobot gamma untuk menghitung rata-rata di county 26
gamma <- (1 / sum(radon$county == 26)) * ifelse(radon$county == 26, 1, 0)

# 7. Jalankan EBLUP
myeblup <- EBLUP(YS, fixed.part, random.part, reg, con, gamma,
                 weights = NULL, estMSE = TRUE)

# 8. Tampilkan hasil
myeblup$thetaP      # prediksi mean log-radon
## [1] 1.306916
sqrt(myeblup$neMSE) # RMSE naive
## [1] 0.04788248
Prediksi Rata-Rata Kadar Log-Radon (County ke-26)

Prediksi rata-rata kadar log-radon pada county ke-26 dilakukan menggunakan metode EBLUP (Empirical Best Linear Unbiased Prediction), dengan asumsi bahwa data dari rumah ber-basement (basement = 1) di county tersebut tidak tersedia.

  • Hasil prediksi rata-rata log-radon: 1.307
  • Jika dikonversi ke skala asli: sekitar 3.7 unit kadar radon
  • Tingkat ketelitian prediksi tinggi, ditunjukkan oleh nilai galat standar (RMSE) sebesar 0.048
Interpretasi

Metode EBLUP terbukti mampu menghasilkan estimasi yang stabil dan akurat, bahkan dalam kondisi data tidak lengkap. Hal ini karena: - EBLUP memperhitungkan variasi antar-wilayah melalui efek acak, dan - Mengakomodasi pengaruh variabel prediktor secara keseluruhan melalui efek tetap.

Dengan demikian, EBLUP menjadi alat prediksi yang efektif dalam analisis data hierarkis atau data dengan struktur kelompok seperti pada kasus ini.

3.3 Prediksi Menggunakan PLUG-IN

Metode PLUG-IN memberikan fleksibilitas yang lebih tinggi dibanding EBLUP karena memungkinkan prediksi untuk berbagai fungsi karakteristik populasi, seperti: - Rata-rata (mean) - Rata-rata geometrik (geometric mean) - Median

Selain itu, prediksi dilakukan dalam skala asli (picoCurie per liter) melalui proses back-transformation dari skala log.

Berikut adalah kode implementasinya:

library(psych)

# 1. Fungsi untuk menghitung tiga karakteristik di county 26
thetaFun <- function(x) {
  c(
    mean(x[radon$county == 26]),
    geometric.mean(x[radon$county == 26]),
    median(x[radon$county == 26])
  )
}

# 2. Fungsi untuk back-transform dari log ke skala asli
backTransExp <- function(x) exp(x)

# 3. Jalankan PLUG-IN predictor
myplugin <- plugInLMM(
  YS, fixed.part, random.part, reg, con,
  weights = NULL,
  backTrans = backTransExp,
  thetaFun = thetaFun
)

# 4. Tampilkan hasil prediksi
myplugin$thetaP
## [1] 4.553745 3.694761 3.900000
Hasil Prediksi dengan Metode Plug-in (County ke-26)

Dari hasil analisis prediksi, diperoleh:

  • Rata-rata kadar radon: 4.554
  • Prediksi rata-rata geometrik: 3.695
  • Prediksi median kadar radon: 3.900
Interpretasi

Metode Plug-in menunjukkan keunggulannya karena mampu memprediksi berbagai karakteristik populasi, tidak terbatas hanya pada rata-rata.

Prediksi dilakukan dalam skala asli melalui proses back-transformation dari skala logaritmik, sehingga hasilnya langsung dapat digunakan untuk interpretasi nyata di lapangan (misalnya, untuk kebijakan atau tindakan kesehatan lingkungan).

Hal ini menjadikan metode Plug-in sebagai pendekatan yang fleksibel dan informatif, terutama dalam konteks distribusi yang tidak simetris seperti kadar radon.

4. Estimasi Akurasi Prediksi Menggunakan Bootstrap

Setelah menghasilkan prediksi menggunakan EBLUP dan PLUG-IN, kita perlu mengetahui seberapa akurat prediksi tersebut. Kita akan menggunakan metode bootstrap untuk mengestimasi dua ukuran akurasi:

4.1 Residual Bootstrap untuk PLUG-IN Predictor

Kita akan menggunakan metode residual bootstrap dengan koreksi, sesuai dengan rekomendasi dari penulis qape.

# Tentukan jumlah iterasi bootstrap dan kuantil QAPE yang ingin dihitung
B <- 500                      # Jumlah iterasi bootstrap
p <- c(0.75, 0.9)             # Kuantil QAPE: 75% dan 90%

set.seed(1056)                # Reproducible
residBoot <- bootRes(
  predictor = myplugin,
  B = B,
  p = p,
  correction = TRUE
)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00204383 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0071418 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00349342 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
# Tampilkan RMSE estimasi
residBoot$estRMSE
## [1] 0.1848028 0.2003681 0.2824359
# Tampilkan QAPE estimasi
residBoot$estQAPE
##          [,1]      [,2]      [,3]
## 75% 0.1533405 0.2135476 0.2908988
## 90% 0.2813886 0.3397411 0.4374534
Evaluasi Akurasi Prediksi (RMSE & QAPE)
RMSE Estimasi

Nilai RMSE per domain: - Domain 1: 0.1848 → prediksi sangat akurat
- Domain 2: 0.2004 → prediksi cukup akurat
- Domain 3: 0.2824 → galat prediksi relatif besar

Semakin kecil RMSE, semakin baik akurasi model.

QAPE Estimasi (Galat Absolut Kuantil)

Kuantil ke-75 (75% galat < nilai berikut): - Domain 1: 0.1533
- Domain 2: 0.2135
- Domain 3: 0.2909

Kuantil ke-90 (90% galat < nilai berikut): - Domain 1: 0.2814
- Domain 2: 0.3397
- Domain 3: 0.4375

Kesimpulan

Metode Plug-in predictor cukup akurat, terutama pada domain 1 (nilai RMSE & QAPE paling rendah).
Sebaliknya, domain 3 menunjukkan akurasi lebih rendah.
Bootstrap residual dengan koreksi terbukti memberikan estimasi akurasi yang realistis dan layak dijadikan dasar evaluasi model.

5. Perbandingan Model yang Salah Spesifikasi

Tujuan dari bagian ini adalah untuk membandingkan akurasi prediksi antara: - Model awal (myplugin) → model lengkap (benar) - Model yang lebih sederhana (myplugin.mis) → model salah spesifikasi

Kita tetap menggunakan metode residual bootstrap dengan koreksi untuk mengestimasi RMSE dan QAPE dari kedua model tersebut.

5.1 Definisikan Model Sederhana

# Model lebih sederhana: hanya intercept acak per county
fixed.part.mis <- '1'
random.part.mis <- '(1|county)'

# Estimasi PLUG-IN dari model salah spesifikasi
myplugin.mis <- plugInLMM(
  YS, fixed.part.mis, random.part.mis, reg, con,
  weights = NULL, backTrans = backTransExp, thetaFun = thetaFun
)

5.2 Bandingkan Akurasi Kedua Model

Setelah mendefinisikan dua model — model benar (myplugin) dan model salah spesifikasi (myplugin.mis) — kita gunakan residual bootstrap dengan koreksi untuk membandingkan akurasi prediksi keduanya.

# Menentukan jumlah iterasi dan kuantil QAPE (jika belum ditentukan sebelumnya)
B <- 500
p <- c(0.75, 0.9)

# Jalankan residual bootstrap perbandingan
set.seed(1056)
residBootMis <- bootResMis(
  predictorLMM = myplugin,
  predictorLMMmis = myplugin.mis,
  B = B,
  p = p,
  correction = TRUE
)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00204383 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0071418 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00349342 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00204383 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0071418 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00349342 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
# Tampilkan hasil RMSE
residBootMis$estRMSElmm      # Model benar
## [1] 0.1848028 0.2003681 0.2824359
residBootMis$estRMSElmmMis   # Model salah spesifikasi
## [1] 0.1919184 0.3192304 0.2762137
# Tampilkan hasil QAPE
residBootMis$estQAPElmm      # Model benar
##          [,1]      [,2]      [,3]
## 75% 0.1533405 0.2135476 0.2908988
## 90% 0.2813886 0.3397411 0.4374534
residBootMis$estQAPElmmMis   # Model salah spesifikasi
##          [,1]      [,2]      [,3]
## 75% 0.2267062 0.3802836 0.3255197
## 90% 0.2813787 0.4970726 0.4489399
Perbandingan Akurasi Dua Model (RMSE dan QAPE)
RMSE
  • Model benar (myplugin) menunjukkan galat prediksi yang lebih kecil:
    • Domain 1: 0.1848
    • Domain 2: 0.2004
    • Domain 3: 0.2824
  • Model salah (myplugin.mis) menghasilkan galat lebih besar, terutama pada domain 2:
    • Domain 1: 0.1919
    • Domain 2: 0.3192
    • Domain 3: 0.2762
QAPE (Kuantil Galat Absolut)

Kuantil 75% (75% galat di bawah nilai berikut): - Model benar: 0.1533 (D1), 0.2135 (D2), 0.2909 (D3)
- Model salah: 0.2267 (D1), 0.3803 (D2), 0.3255 (D3)

Kuantil 90% (90% galat di bawah nilai berikut): - Model benar: 0.2814 (D1), 0.3397 (D2), 0.4375 (D3)
- Model salah: 0.2814 (D1), 0.4971 (D2), 0.4489 (D3)

Kesimpulan

Model myplugin lebih akurat dibandingkan myplugin.mis, baik dari sisi RMSE maupun QAPE.
Perbedaan paling mencolok terlihat pada domain 2, di mana model salah menghasilkan galat yang jauh lebih besar.
Pemilihan model yang tepat sangat penting untuk menghasilkan prediksi yang akurat dan stabil.

6. Simulasi Monte Carlo: Pengujian Konsistensi Prediksi

Tujuan dari simulasi Monte Carlo ini adalah: - Menilai akurasi sejati dari metode prediksi - Melihat bias relatif dan RMSE sejati - Mengestimasi kuantil dari error (QAPE) secara stabil

Langkah-langkah ini dilakukan dengan menggunakan seluruh data populasi (radon) untuk: 1. Mengestimasi parameter model 2. Mensimulasikan ratusan populasi buatan 3. Melakukan prediksi ulang dan mengukur error tiap iterasi

6.1 Menjalankan Simulasi Monte Carlo

# Menentukan objek-objek prediktor untuk digunakan dalam simulasi
predictorLMMmis <- myplugin     # Digunakan untuk definisi model (populasi buatan)
predictorLMM <- myplugin        # Digunakan untuk menilai akurasi
predictorLMM2 <- myplugin       # (bisa sama untuk uji tunggal)

# Menentukan parameter simulasi
K <- 500           # Jumlah iterasi Monte Carlo
p <- c(0.75, 0.9)  # Kuantil untuk QAPE
ratioR <- 1        # Tanpa modifikasi residual variance
ratioG <- 1        # Tanpa modifikasi random effects variance

# Jalankan simulasi
set.seed(1086)
MC <- mcLMMmis(
  Ypop = radon$log.radon,
  predictorLMMmis = predictorLMMmis,
  predictorLMM = predictorLMM,
  predictorLMM2 = predictorLMM2,
  K = K,
  p = p,
  ratioR = ratioR,
  ratioG = ratioG
)
## Warning: package 'future' was built under R version 4.4.3
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00217376 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00419829 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0469687 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0023342 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00263891 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0023342 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00263891 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
# Lihat hasil simulasi
MC$rBlmm        # Bias relatif (%)
## [1] -1.73208393 -0.04053178 -5.22355236
MC$rRMSElmm     # RMSE relatif (%)
## [1] 3.429465 4.665810 7.146678
MC$QAPElmm      # QAPE 75% dan 90%
##          [,1]      [,2]      [,3]
## 75% 0.1491262 0.1989504 0.2919221
## 90% 0.2895684 0.2959457 0.4728064
Kesimpulan Simulasi Monte Carlo

Hasil simulasi menunjukkan bahwa:

  • Bias relatif kecil: -1.73% (D1), -0.04% (D2), -5.22% (D3) → prediksi tidak bias secara signifikan.
  • RMSE relatif rendah: 3.43% (D1), 4.67% (D2), 7.15% (D3) → prediksi cukup akurat.
  • QAPE:
    • Kuantil 75%: 0.15 (D1), 0.20 (D2), 0.29 (D3)
    • Kuantil 90%: 0.29 (D1), 0.30 (D2), 0.47 (D3)

Kesimpulan:
Metode prediksi bekerja dengan baik di semua domain, terutama di domain 1 dan 2.
Domain 3 sedikit kurang akurat, tetapi hasil masih dapat diterima.

7. Kesimpulan

Dalam studi ini, telah dilakukan penerapan dan evaluasi metode prediksi berbasis Linear Mixed Model (LMM) pada data kadar radon rumah di berbagai county. Analisis berfokus pada prediksi karakteristik di county ke-26 dengan menggunakan pendekatan EBLUP dan PLUG-IN predictor. Berikut poin-poin kesimpulan utama:

7.1 Prediksi

  • Metode EBLUP berhasil memprediksi rata-rata kadar log-radon di county ke-26 sebesar 1.307 (setara dengan 3.7 dalam skala asli).
  • Metode PLUG-IN mampu memberikan prediksi lebih lengkap, termasuk:
    • Rata-rata: 4.554
    • Geometric mean: 3.695
    • Median: 3.900
  • PLUG-IN unggul karena bisa digunakan untuk prediksi karakteristik non-linier dan bekerja langsung dalam skala asli.

7.2 Evaluasi Akurasi (RMSE & QAPE)

  • Residual bootstrap dengan koreksi digunakan untuk mengukur galat prediksi.
  • Hasil menunjukkan bahwa model benar (myplugin) memiliki RMSE dan QAPE lebih kecil dibanding model salah (myplugin.mis), terutama di domain 2.
  • PLUG-IN menunjukkan akurasi prediksi cukup baik, khususnya:
    • RMSE model benar: 0.18–0.28
    • QAPE 75%: < 0.29, dan 90%: < 0.44 di semua domain.

7.3 Model Salah Spesifikasi

  • Model dengan spesifikasi salah (tanpa efek slope acak) menunjukkan akurasi lebih rendah.
  • Namun, model tersebut tetap dapat digunakan jika komputasi perlu disederhanakan dan ketelitian masih dapat ditoleransi.

7.4 Simulasi Monte Carlo

  • Hasil simulasi Monte Carlo menunjukkan bahwa metode PLUG-IN:
    • Memiliki bias relatif kecil (sekitar -1.7% hingga -5.2%)
    • RMSE relatif rendah (sekitar 3.4% hingga 7.1%)
    • QAPE stabil pada kuantil 75% dan 90%
  • Hal ini menunjukkan bahwa metode cukup andal dan akurat dalam banyak pengulangan.

7.5 Rekomendasi

  • Gunakan PLUG-IN predictor untuk prediksi karakteristik non-linier (seperti median dan geometric mean) dalam skala asli.
  • Gunakan EBLUP untuk prediksi linier dan interpretasi sederhana.
  • Penting untuk melakukan evaluasi akurasi menggunakan bootstrap dan simulasi Monte Carlo agar prediksi lebih meyakinkan.