1 Motivasi EBLUP

Dapat diterapkan pada taksiran pada Model Linier Campuran (Linear mixed model)
Pengembangan dari Best Linear Unbiased Prediction (BLUP)
Dikembangkan untuk masalah pendugaan area kecil (SAE)
Pendugaan ragam melalui data contoh (sampel) seperti pada model Fay-Herriot (FH)
Mensubstitusikan komponen ragam yang tidak diketahui dengan nilai penduganya pada model dugaan BLUP

2 Contoh 1

Contoh pertama menggunakan data yang berasal dari data survey jagung dan kacang kedelai, dan data satelit di 12 Kabupaten di Iowa.
Dataset ini berisi informasi mengenai produksi jagung dan kedelai di berbagai area kecil.
Data ini tersedia di R pada package sae. Akan dilakukan pendugaan area kecil terhadap data ini.

Dataset tersebut merupakan data yang terdiri atas 37 observasi dengan 5 variabel.

County: numeric county code.
CornHec: reported hectares of corn from the survey.
SoyBeansHec: reported hectares of soy beans from the survey.
CornPix: number of pixels of corn in sample segment within county, from satellite data.
SoyBeansPix: number of pixels of soy beans in sample segment within county, from satellite data.

Penggunaan dalam Model

Model Regresi Linear Biasa

\[\text{Cornhec} = \beta_0 + \beta_1 \text{CornPix} + \epsilon \]

Model EBLUP (SAE)

\[ \text{Cornhec} = \beta_0 + \beta_1 \text{CornPix} + u_{\text{County}} + \epsilon \]

Efek Tetap: CornPix, SoyBeansPix
Efek Acak: County (untuk mengakomodasi variasi antar area kecil)

rm(list=ls())   # Menghapus semua dataset dan variabel
graphics.off()  # Menutup semua gambar

library(sae)

## Warning: package 'MASS' was built under R version 4.4.3

## Warning: package 'lme4' was built under R version 4.4.3

## Warning: package 'Matrix' was built under R version 4.4.3

data("cornsoybean")
cornsoybean

# Melihat struktur dan beberapa baris pertama dari dataset
str(cornsoybean)

## 'data.frame':    37 obs. of  5 variables:
##  $ County     : int  1 2 3 4 4 5 5 5 6 6 ...
##  $ CornHec    : num  165.8 96.3 76.1 185.3 116.4 ...
##  $ SoyBeansHec: num  8.09 106.03 103.6 6.47 63.82 ...
##  $ CornPix    : int  374 209 253 432 367 361 288 369 206 316 ...
##  $ SoyBeansPix: int  55 218 250 96 178 137 206 165 218 221 ...

data("cornsoybeanmeans")
cornsoybeanmeans

# Melihat struktur dan beberapa baris pertama dari dataset
str(cornsoybeanmeans)

## 'data.frame':    12 obs. of  6 variables:
##  $ CountyIndex          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ CountyName           : Factor w/ 12 levels "CerroGordo","Franklin",..: 1 3 11 6 2 8 10 12 9 4 ...
##  $ SampSegments         : int  1 1 1 2 3 3 3 3 4 5 ...
##  $ PopnSegments         : num  545 566 394 424 564 570 402 567 687 569 ...
##  $ MeanCornPixPerSeg    : num  295 300 290 291 318 ...
##  $ MeanSoyBeansPixPerSeg: num  190 197 205 220 188 ...

Disisi lain, didata tersebut juga terdapat data auxialary (data pendukung) pada file cornsoybeanmeans.

cornsoybeanmeans: Data frame dengan 12 observasi dan 6 variabel.

CountyIndex:numeric county code.
CountyName: name of the county.
SampSegments: number of sample segments in the county (sample size).
PopnSegments: number of population segments in the county (population size).
MeanCornPixPerSeg: mean number of corn pixels per segment in the county.
MeanSoyBeansPixPerSeg: mean number of soy beans pixels per segment in the county.

# Mengambil beberapa variabel
(Xmean <- data.frame(cornsoybeanmeans[, c("CountyIndex","MeanCornPixPerSeg","MeanSoyBeansPixPerSeg")]))

# Variabel Xmean sebenarnya digunakan sebagai data auxiliary dalam model EBLUP (Empirical Best Linear Unbiased Prediction)
# terutama saat menggunakan model BHF (Battese, Harter, and Fuller) dalam fungsi pbmseBHF dari paket sae.

(Popn <- data.frame(cornsoybeanmeans[, c("CountyIndex","PopnSegments")]))

# popnsize = Popn: Memastikan bahwa ukuran populasi di setiap county diperhitungkan dalam model.
# Tanpa Popn, model mungkin akan menganggap setiap area kecil memiliki bobot yang sama,
# padahal kenyataannya populasi setiap area bisa sangat berbeda.

# Memanggil fungsi untuk menghitung estimasi EBLUP (Empirical Best Linear Unbiased Prediction) dari rata-rata luas area tanaman jagung di semua kabupaten, serta menghitung estimasi Mean Squared Error (MSE) menggunakan metode bootstrap parametrik dengan 100 pengulangan (replicates), tanpa menghapus data apapun. 
# Proses fitting model menggunakan metode default yaitu REML (Restricted Maximum Likelihood).

set.seed(2025)
BHF <- pbmseBHF(CornHec ~ CornPix + SoyBeansPix, 
                dom = County,
                meanxpop = Xmean, 
                popnsize = Popn,
                B = 100, data = cornsoybean) #y = CorHec , x1=CornPix, x2=SoyBeansPix

## 
## Bootstrap procedure with B = 100 iterations starts.
## b = 1 
## b = 2 
## b = 3 
## b = 4 
## b = 5

## boundary (singular) fit: see help('isSingular')

## b = 6 
## b = 7 
## b = 8 
## b = 9 
## b = 10 
## b = 11 
## b = 12 
## b = 13 
## b = 14

## boundary (singular) fit: see help('isSingular')

## b = 15 
## b = 16 
## b = 17 
## b = 18 
## b = 19 
## b = 20 
## b = 21

## boundary (singular) fit: see help('isSingular')

## b = 22 
## b = 23

## boundary (singular) fit: see help('isSingular')

## b = 24

## boundary (singular) fit: see help('isSingular')

## b = 25 
## b = 26 
## b = 27 
## b = 28

## boundary (singular) fit: see help('isSingular')

## b = 29

## boundary (singular) fit: see help('isSingular')

## b = 30 
## b = 31

## boundary (singular) fit: see help('isSingular')

## b = 32 
## b = 33 
## b = 34 
## b = 35

## boundary (singular) fit: see help('isSingular')

## b = 36 
## b = 37 
## b = 38 
## b = 39 
## b = 40

## boundary (singular) fit: see help('isSingular')

## b = 41 
## b = 42 
## b = 43 
## b = 44 
## b = 45

## boundary (singular) fit: see help('isSingular')

## b = 46 
## b = 47 
## b = 48 
## b = 49 
## b = 50 
## b = 51

## boundary (singular) fit: see help('isSingular')

## b = 52 
## b = 53 
## b = 54

## boundary (singular) fit: see help('isSingular')

## b = 55 
## b = 56 
## b = 57

## boundary (singular) fit: see help('isSingular')

## b = 58 
## b = 59

## boundary (singular) fit: see help('isSingular')

## b = 60 
## b = 61 
## b = 62 
## b = 63 
## b = 64 
## b = 65 
## b = 66 
## b = 67

## boundary (singular) fit: see help('isSingular')

## b = 68 
## b = 69 
## b = 70 
## b = 71 
## b = 72 
## b = 73 
## b = 74 
## b = 75 
## b = 76

## boundary (singular) fit: see help('isSingular')

## b = 77 
## b = 78

## boundary (singular) fit: see help('isSingular')

## b = 79 
## b = 80 
## b = 81

## boundary (singular) fit: see help('isSingular')

## b = 82 
## b = 83 
## b = 84 
## b = 85 
## b = 86 
## b = 87 
## b = 88

## boundary (singular) fit: see help('isSingular')

## b = 89 
## b = 90 
## b = 91

## boundary (singular) fit: see help('isSingular')

## b = 92 
## b = 93 
## b = 94 
## b = 95 
## b = 96

## boundary (singular) fit: see help('isSingular')

## b = 97 
## b = 98 
## b = 99 
## b = 100

# Menampilkan hasil EBLUP
print(BHF)

## $est
## $est$eblup
##    domain    eblup sampsize
## 1       1 122.5825        1
## 2       2 123.5274        1
## 3       3 113.0343        1
## 4       4 114.9901        2
## 5       5 137.2660        3
## 6       6 108.9807        3
## 7       7 116.4839        3
## 8       8 122.7711        3
## 9       9 111.5648        4
## 10     10 124.1565        5
## 11     11 112.4626        5
## 12     12 131.2515        6
## 
## $est$fit
## $est$fit$summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ys ~ -1 + Xs + (1 | dom)
## 
## REML criterion at convergence: 322
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9288 -0.5711  0.1041  0.5953  1.5333 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  dom      (Intercept)  63.31    7.957  
##  Residual             297.71   17.254  
## Number of obs: 37, groups:  dom, 12
## 
## Fixed effects:
##                 Estimate Std. Error t value
## XsXs(Intercept) 17.96398   30.97450   0.580
## XsXsCornPix      0.36634    0.06496   5.640
## XsXsSoyBeansPix -0.03036    0.06758  -0.449
## 
## Correlation of Fixed Effects:
##             XsX(I) XsXsCP
## XsXsCornPix -0.947       
## XsXsSyBnsPx -0.897  0.734
## 
## $est$fit$fixed
## XsXs(Intercept)     XsXsCornPix XsXsSoyBeansPix 
##      17.9639791       0.3663352      -0.0303638 
## 
## $est$fit$random
##    (Intercept)
## 1     2.184574
## 2     1.475118
## 3    -4.730863
## 4    -2.764825
## 5     8.370915
## 6     4.274827
## 7    -2.705540
## 8     1.156682
## 9     5.026852
## 10   -2.883398
## 11   -8.652532
## 12   -0.751808
## 
## $est$fit$errorvar
## [1] 297.7128
## 
## $est$fit$refvar
## [1] 63.3149
## 
## $est$fit$loglike
## [1] -161.0058
## 
## $est$fit$residuals
##  [1]  10.2720798   6.9361476 -22.2449806  14.8089506 -27.8094281   7.6579280
##  [7]  26.4555017   5.2474324   1.7954444  18.6496604  -0.3444511 -14.3508799
## [13]  14.6861461 -13.0569810  21.4594809 -10.4389338  -5.5817170  -6.1470789
## [19]  16.8133454   7.2348743   5.7356095  -9.8543869   3.5991772 -11.2943444
## [25]   9.4655170  -5.4739832   8.2792213 -17.6896617 -17.4938455  -6.4118846
## [31]  -7.3688850  -1.2606073 -50.5344992  22.9470294  -2.1152674  13.3886503
## [37]  14.0396194
## 
## 
## 
## $mse
##    domain      mse
## 1       1 69.28657
## 2       2 71.60253
## 3       3 90.37036
## 4       4 62.73799
## 5       5 64.72651
## 6       6 54.81490
## 7       7 43.70106
## 8       8 56.32491
## 9       9 39.13320
## 10     10 41.72201
## 11     11 32.41771
## 12     12 30.46472

# Menampilkan MSE hasil EBLUP
BHF$mse

Penjelasan sintaks diatas adalah sebagai berikut.

1.Formula Model: CornHec ~ CornPix + SoyBeansPix

Respon (y): CornHec (luas area jagung dalam hektar).
Prediktor (x1, x2): CornPix (piksel jagung) dan SoyBeansPix (piksel kedelai).

2.dom = County

County sebagai variabel domain atau area kecil.
EBLUP akan dihitung untuk setiap County dalam dataset.

3.meanxpop = Xmean

Xmean berisi rata-rata populasi untuk setiap variabel prediktor dalam setiap domain (County).
Membantu dalam menghitung estimasi untuk area kecil dengan memanfaatkan data tambahan.

4.popnsize = Popn

Popn: Ukuran populasi dalam setiap County.
Penting untuk menyesuaikan estimasi EBLUP pada skala populasi.

5.B = 100

Menggunakan 100 replikasi bootstrap untuk memperkirakan MSE EBLUP.
Bootstrap memberikan estimasi variabilitas yang lebih baik, terutama pada data area kecil.

#Menghitung Standard Error dari Estimasi EBLUP menggunakan akar kkuadrat dari MSE
sqrtmse.BHF <- sqrt(BHF$mse$mse)
rescorn <- data.frame(CountyIndex = BHF$est$eblup[, "domain"],
                       CountyName = cornsoybeanmeans[, "CountyName"],
                       BHF$est$eblup[, c("sampsize", "eblup")],sqrtmse=sqrtmse.BHF)


#Menampilkan hasil
print(rescorn, row.names = FALSE)

##  CountyIndex CountyName sampsize    eblup  sqrtmse
##            1 CerroGordo        1 122.5825 8.323855
##            2   Hamilton        1 123.5274 8.461828
##            3      Worth        1 113.0343 9.506333
##            4   Humboldt        2 114.9901 7.920732
##            5   Franklin        3 137.2660 8.045279
##            6 Pocahontas        3 108.9807 7.403708
##            7  Winnebago        3 116.4839 6.610678
##            8     Wright        3 122.7711 7.504992
##            9    Webster        4 111.5648 6.255653
##           10    Hancock        5 124.1565 6.459257
##           11    Kossuth        5 112.4626 5.693656
##           12     Hardin        6 131.2515 5.519486

MSE disimpan dalam BHF$mse$mse

Standar Error (SE) dihitung dengan akar kuadrat MSE, karena

\[SE=\sqrt{MSE}\]

Menghasilkan SE EBLUP setiap County.

# Membuat plot perbandingan Standar Error EBLUP dan Direct Estimation
plot(rescorn$CountyIndex, rescorn$sqrtmse, type = "o", col = "blue", pch = 16,
     ylab = "Standar Error", xlab = "County", main = "Standar Error EBLUP (sqrt MSE)")

Dapat dilihat nilai eblup dari masing-masing Kabupaten. Berdasarkan hasil tersebut terlihat bahwasanya untuk Kabupaten dengan sample size yang lebih besar menghasilkan nilai sqrtmse yang lebih kecil.

# Menggabungkan data rescorn dengan data asli cornsoybean berdasarkan CountyIndex
rescorn <- merge(rescorn, cornsoybean[, c("County", "CornHec")], 
                 by.x = "CountyIndex", by.y = "County", all.x = TRUE)

# Mengubah nama kolom untuk memudahkan interpretasi
colnames(rescorn)[ncol(rescorn)] <- "y_aktual"

# Menambahkan kolom y_prediksi menggunakan kolom eblup
rescorn$y_prediksi <- rescorn$eblup

# Menampilkan data frame lengkap dengan y prediksi dan y aktual
print(rescorn, row.names = FALSE)

##  CountyIndex CountyName sampsize    eblup  sqrtmse y_aktual y_prediksi
##            1 CerroGordo        1 122.5825 8.323855   165.76   122.5825
##            2   Hamilton        1 123.5274 8.461828    96.32   123.5274
##            3      Worth        1 113.0343 9.506333    76.08   113.0343
##            4   Humboldt        2 114.9901 7.920732   185.35   114.9901
##            4   Humboldt        2 114.9901 7.920732   116.43   114.9901
##            5   Franklin        3 137.2660 8.045279   162.08   137.2660
##            5   Franklin        3 137.2660 8.045279   152.04   137.2660
##            5   Franklin        3 137.2660 8.045279   161.75   137.2660
##            6 Pocahontas        3 108.9807 7.403708    92.88   108.9807
##            6 Pocahontas        3 108.9807 7.403708   149.94   108.9807
##            6 Pocahontas        3 108.9807 7.403708    64.75   108.9807
##            7  Winnebago        3 116.4839 6.610678   127.07   116.4839
##            7  Winnebago        3 116.4839 6.610678   133.55   116.4839
##            7  Winnebago        3 116.4839 6.610678    77.70   116.4839
##            8     Wright        3 122.7711 7.504992   206.39   122.7711
##            8     Wright        3 122.7711 7.504992   108.33   122.7711
##            8     Wright        3 122.7711 7.504992   118.17   122.7711
##            9    Webster        4 111.5648 6.255653    99.96   111.5648
##            9    Webster        4 111.5648 6.255653   140.43   111.5648
##            9    Webster        4 111.5648 6.255653    98.95   111.5648
##            9    Webster        4 111.5648 6.255653   131.04   111.5648
##           10    Hancock        5 124.1565 6.459257   114.12   124.1565
##           10    Hancock        5 124.1565 6.459257   100.60   124.1565
##           10    Hancock        5 124.1565 6.459257   127.88   124.1565
##           10    Hancock        5 124.1565 6.459257   116.90   124.1565
##           10    Hancock        5 124.1565 6.459257    87.41   124.1565
##           11    Kossuth        5 112.4626 5.693656    93.48   112.4626
##           11    Kossuth        5 112.4626 5.693656   121.00   112.4626
##           11    Kossuth        5 112.4626 5.693656   109.91   112.4626
##           11    Kossuth        5 112.4626 5.693656   122.66   112.4626
##           11    Kossuth        5 112.4626 5.693656   104.21   112.4626
##           12     Hardin        6 131.2515 5.519486    88.59   131.2515
##           12     Hardin        6 131.2515 5.519486    88.59   131.2515
##           12     Hardin        6 131.2515 5.519486   165.35   131.2515
##           12     Hardin        6 131.2515 5.519486   104.00   131.2515
##           12     Hardin        6 131.2515 5.519486    88.63   131.2515
##           12     Hardin        6 131.2515 5.519486   153.70   131.2515

# Membuat plot untuk membandingkan y prediksi dan y aktual
plot(rescorn$y_aktual, type = "o", col = "red", pch = 16, 
     ylab = "Nilai", xlab = "County", main = "Perbandingan y Prediksi vs y Aktual")
lines(rescorn$y_prediksi, type = "o", col = "blue", pch = 16)

# Menambahkan legenda
legend("topright", legend = c("y Aktual", "y Prediksi (EBLUP)"), 
       col = c("red", "blue"), pch = 16, lty = 1)

3 Contoh 2

Dipunyai data Susenas tahun dan data Potensi Desa (PODES) tahun 2008 yang diperoleh dari data Badan Pusat Statistik (BPS) dengan 35 sampel desa di kabupaten Jember. Jumlah rumah tangga sampel sebanyak 549 RT.

Tabel Variabel Penelitian
No	Variabel	Keterangan
1	Y	Rata-rata pengeluaran rumah tangga per kapita per bulan di desa se kabupaten Jember
2	$X_1$	Persentase keluarga yang menggantungkan hidupnya pada pertanian
3	$X_2$	Jumlah keluarga yang menerima askeskin dalam satu tahun terakhir
4	$X_3$	Jumlah keluarga pengguna listrik PLN
5	$X_4$	Jumlah keluarga yang mengenyam pendidikan SD, SMP, SMA, dan PT
6	$X_5$	Jumlah keluarga yang tinggal di pemukiman kumuh
7	$X_6$	Jumlah keluarga yang memiliki Surat Keterangan Tidak Mampu (SKTM) dan satu tahun terakhir
8	$X_7$	Jumlah keluarga yang pernah belajar di lembaga pendidikan dan keterampilan
9	$X_8$	Jumlah keluarga yang memiliki anggota keluarga sebagai Tenaga Kerja Indonesia

Permasalahan

Pendugaan langsung seringkali menimbulkan error yang besar karena menggunakan sampel desa kecil sehingga dibutuhkan model SAE sebagai pendugaan tidak langsung dengan memanfaatkan informasi peubah penyerta (auxiliary).

Carilah:

Taksiran pengeluaran rumah tangga per kapita per bulan di desa di Kabupaten Jember ($\tilde {\theta}_i^{EBLUP}$)
Taksiran pengaruh tetap ($\hat{\beta}$)
g1, g2, dan g3
Taksiran Mean Squared Error (MSE)
Taksiran komponen acak** ($\hat{v}$)
Perbandingan SE pendugaan langsung dengan pendugaan EBLUP

Penyelesaian

# Install Package
rm(list=ls())   # Menghapus semua dataset dan variabel
graphics.off()  # Menutup semua gambar
library(readxl) #membaca data excel

## Warning: package 'readxl' was built under R version 4.4.3

library(dplyr) #mmebuat plot
library(ggplot2) #membuat plot

# Input Data
pengeluaran <- read_excel("pengeluaran.xlsx")
pengeluaran=as.matrix(pengeluaran)
prediktor<-read_excel("prediktor.xlsx")
prediktor=as.matrix(prediktor)

#Pelabelan peubah
ydirect<-pengeluaran[,1] #y dari pendugaan langsung
Xpop<-prediktor #prediktor atau peubah penyerta
vardirect=pengeluaran[,2]^2 #MSE
ydirect=as.matrix(ydirect) #Matriks y
Xpop=as.matrix(Xpop) #matriks Z
vardirect=as.matrix(vardirect) #Matriks MSE
m=35 #Jumlah Desa

#EBLUB dengan metode ML, scoring algorithm
EBLUP.area.ML <- function(ydirect, Xpop, vardirect, m, 
                          tol=10e-5, maxiter=50)
{ 
  sigma2.v.estimate <- 0 #komponen variansi
  sigma2.v.estimate[1] <- 5000000 #inisial value
  k <- 0 
  diff <- tol+1 
  while ((diff>tol) & (k<maxiter)) 
  { 
    k <- k+1
    V <- diag(1,m) #matriks varians-kovarians
    
    for (i in 1:m) 
    { 
      V[i,i] <- 
        (sigma2.v.estimate[k]+vardirect[i,1]) 
    } 
   
 #vardirect merupakan vektor variansi sampelberukuran mx1 
    Vinvers <- solve(V)
    Beta <- 
      solve(t(Xpop)%*%Vinvers%*%Xpop)%*%t(Xpop)%*%Vinvers%*%ydirect[,1] #ydirect merupakan vektor dari taksiran langsung,berukuran mx1
    sdev <-(-0.5)*sum(diag(Vinvers))-(0.5)*t(ydirect[,1]-(Xpop)%*%Beta)%*%((-1)*Vinvers%*%Vinvers)%*%(ydirect[,1]-(Xpop)%*%Beta) #Xpop adalah matriks dari auxiliary variables, berukuran mxp
    Idev <-((0.5)*(sum(diag(Vinvers%*%Vinvers))))^(-1) #Idev is invers dari matriks informasi
    #scoring algorithm
    sigma2.v.estimate[k+1] <- 
      sigma2.v.estimate[k]+Idev*sdev 
    diff <- abs(sigma2.v.estimate[k+1]-
                  sigma2.v.estimate[k]) #batas toleransi 
  } 
  
 
 if (sigma2.v.estimate[k+1]<0) 
    sigma2v <- 0
  else
    sigma2v <- sigma2.v.estimate[k+1]
  
  V <- diag(1,m)
  G <- diag(1,m)*sigma2v
  
  for (i in 1:m)
  { 
    V[i,i] <-(sigma2v+vardirect[i,1]) 
  } 
  Vinvers <- solve(V)
  Bestimate <- 
    solve(t(Xpop)%*%Vinvers%*%Xpop)%*%t(Xpop)%*%Vinvers%*%
    ydirect[,1]
  
  m1 <- diag(1,m)
  thetaEBLUP <- 
    Xpop%*%Bestimate+m1%*%G%*%Vinvers%*%(ydirect[,1]-(Xpop%*%Bestimate)) #EBLUP estimator
  efekrandom <- m1%*%G%*%Vinvers%*%(ydirect[,1]-(Xpop%*%Bestimate))#Taksiran efek random
  #MSE EBLUP small area
  #g1
  g1 <- diag((G-G%*%Vinvers%*%G))
  #g2
  g2 <- matrix(0,m,1)
  for (i in 1:m) 
  { 
    m1 <- matrix(0,m,1)
    m1[i] <- 1 
    g2[i] <-(Xpop[i,]-t(m1)%*%G%*%Vinvers%*%Xpop)%*%solve(t(Xpop)%*%Vinvers%*%Xpop)%*%t(Xpop[i,]-t(m1)%*%G%*%Vinvers%*%Xpop) 
  }   
 Idev <-((0.5)*(sum(diag(Vinvers%*%Vinvers))))^(-1)
  #g3
  g3 <- matrix(0,m,1)
  for (i in 1:m) 
  { m1 <- matrix(0,m,1)
  m1[i] <- 1 
  g3[i] <-(t(m1)%*%(Vinvers+G%*%((-1)*Vinvers%*%Vinvers))%*%V%*%t((t(m1)%*%(Vinvers+G%*%((-1)*Vinvers%*%Vinvers)))))*Idev 
  } 
  #bias
  bdist <- 0 
  btr <- 0 
  gradg1 <- 0 
  I <- diag(1,m)
  distorsione <- matrix(0,m,1)
  for (i in 1:m) 
  { 
    
m1 <- matrix(0,m,1)
    m1[i] <- 1 
    gradg1 <- t(m1)%*%(I-((I%*%Vinvers%*%G)+(G%*%((-1)*Vinvers%*%I%*%Vinvers)%*%G)+(G%*%Vinvers%*%I)))%*%m1
    btr <- 
      sum(diag(solve(t(Xpop)%*%Vinvers%*%Xpop)%*%t(Xpop)%*%((-1)*Vinvers%*%I%*%Vinvers)%*%Xpop))
    bdist <- (1/(m*2))*Idev*btr
    distorsione[i,1] <- bdist*gradg1 
  } 
#Taksiran MSE EBLUP
  mseestimate <- g1-distorsione+g2+2*g3 
  #list(beta=Bestimate, EBLUP=thetaEBLUP, g1=g1, g2=g2,g3=g3, 
  #MSE=mseestimate)
  list(EBLUP=thetaEBLUP, beta=Bestimate, 
       sigma2v=sigma2.v.estimate[k+1],
       g1=g1, g2=g2, g3=g3, mse=mseestimate, 
       efekrandom=efekrandom)
} 
eblup=EBLUP.area.ML(ydirect, Xpop, vardirect, m, tol=10e-5, 
              maxiter=50)
eblup #print estimasi eblup

## $EBLUP
##           [,1]
##  [1,] 668349.8
##  [2,] 655743.6
##  [3,] 577792.5
##  [4,] 580859.3
##  [5,] 596793.0
##  [6,] 686255.2
##  [7,] 678949.9
##  [8,] 656077.9
##  [9,] 627436.6
## [10,] 620438.7
## [11,] 624152.3
## [12,] 608647.0
## [13,] 668187.1
## [14,] 623690.7
## [15,] 657171.8
## [16,] 615005.6
## [17,] 631231.9
## [18,] 775181.2
## [19,] 633064.7
## [20,] 588383.7
## [21,] 609636.9
## [22,] 577281.0
## [23,] 608394.0
## [24,] 566551.5
## [25,] 605534.2
## [26,] 660425.9
## [27,] 730834.3
## [28,] 586503.1
## [29,] 596075.3
## [30,] 586283.5
## [31,] 568858.3
## [32,] 707470.8
## [33,] 513974.3
## [34,] 516977.3
## [35,] 481776.3
## 
## $beta
##          [,1]
## x1 -63874.548
## x2 -47021.432
## x3  59036.835
## x4  -6977.181
## x5   3048.766
## x6 -51071.893
## x7 -17283.963
## x8 -59079.144
## 
## $sigma2v
## [1] 2.80108e+11
## 
## $g1
##  [1]  2070782720  1424658278   397716341   246555468   149468640  1968718789
##  [7]  1249824985   738107503  1025868786   611028701  1286965549   289211306
## [13]  1035894860   446904486   867085630   460546019   462303093  7878984109
## [19]   744724992   178269208   376319486   118498393   334636159    63220032
## [25]   605473081   927837512   814896114  1013229875  1179200550  1416055982
## [31]   375777604   650521938  2442056207 20460702275  4500855847
## 
## $g2
##               [,1]
##  [1,] 5.182709e+06
##  [2,] 3.231516e+06
##  [3,] 6.093301e+04
##  [4,] 5.541615e+04
##  [5,] 7.872998e+03
##  [6,] 1.270364e+06
##  [7,] 4.443142e+05
##  [8,] 4.480948e+05
##  [9,] 6.798062e+05
## [10,] 4.823632e+05
## [11,] 9.695203e+05
## [12,] 1.696555e+05
## [13,] 3.094147e+05
## [14,] 1.517314e+05
## [15,] 3.127642e+05
## [16,] 2.827769e+05
## [17,] 2.055500e+05
## [18,] 1.791288e+07
## [19,] 1.461332e+05
## [20,] 9.192777e+03
## [21,] 1.659199e+05
## [22,] 1.171995e+03
## [23,] 1.803806e+05
## [24,] 8.589199e+02
## [25,] 3.357939e+05
## [26,] 8.619534e+05
## [27,] 4.644401e+05
## [28,] 6.176638e+05
## [29,] 2.919732e+05
## [30,] 2.750442e+06
## [31,] 1.134682e+05
## [32,] 7.046970e+05
## [33,] 8.494305e+06
## [34,] 6.594652e+08
## [35,] 7.796955e+06
## 
## $g3
##               [,1]
##  [1,] 8.786985e+05
##  [2,] 4.168697e+05
##  [3,] 3.260793e+04
##  [4,] 1.253835e+04
##  [5,] 4.609584e+03
##  [6,] 7.945067e+05
##  [7,] 3.210330e+05
##  [8,] 1.121726e+05
##  [9,] 2.164629e+05
## [10,] 7.690742e+04
## [11,] 3.403512e+05
## [12,] 1.724945e+04
## [13,] 2.207068e+05
## [14,] 4.116512e+04
## [15,] 1.547287e+05
## [16,] 4.371443e+04
## [17,] 4.404835e+04
## [18,] 1.245497e+07
## [19,] 1.141903e+05
## [20,] 6.556463e+03
## [21,] 2.919597e+04
## [22,] 2.897574e+03
## [23,] 2.308979e+04
## [24,] 8.249066e+02
## [25,] 7.551675e+04
## [26,] 1.771317e+05
## [27,] 1.366887e+05
## [28,] 2.111716e+05
## [29,] 2.858489e+05
## [30,] 4.118634e+05
## [31,] 2.911201e+04
## [32,] 8.715805e+04
## [33,] 1.220399e+06
## [34,] 8.011111e+07
## [35,] 4.114798e+06
## 
## $mse
##              [,1]
##  [1,]  2077823284
##  [2,]  1428771083
##  [3,]   397846196
##  [4,]   246637385
##  [5,]   149486256
##  [6,]  1971668966
##  [7,]  1250947960
##  [8,]   738792706
##  [9,]  1027006173
## [10,]   611673626
## [11,]  1288654574
## [12,]   289417420
## [13,]  1036670827
## [14,]   447143226
## [15,]   867725465
## [16,]   460921194
## [17,]   462601747
## [18,]  7923261250
## [19,]   745112499
## [20,]   178292258
## [21,]   376547116
## [22,]   118505689
## [23,]   334865342
## [24,]    63222634
## [25,]   605968497
## [26,]   929073896
## [27,]   815649488
## [28,]  1014293933
## [29,]  1180096797
## [30,]  1419677126
## [31,]   375952604
## [32,]   651410865
## [33,]  2453131020
## [34,] 21290197170
## [35,]  4517356977
## 
## $efekrandom
##             [,1]
##  [1,]  166389.08
##  [2,] -109468.74
##  [3,]  367437.45
##  [4,]  763553.25
##  [5,]  346565.43
##  [6,]  426915.75
##  [7,]  504496.02
##  [8,]   96543.53
##  [9,]  790675.94
## [10,]  854161.07
## [11,]  746069.51
## [12,]  343442.90
## [13,]  843733.35
## [14,]  288032.77
## [15,]  579096.22
## [16,]  560106.92
## [17,] 1091283.46
## [18,]  809113.32
## [19,]  674222.09
## [20,]  581511.43
## [21,]  472121.51
## [22,]  420615.13
## [23,]  464819.41
## [24,]  547012.73
## [25,]   33597.31
## [26,]  344855.36
## [27,]  256938.59
## [28,]  -66698.65
## [29,]  489578.41
## [30,] -154403.72
## [31,]  314671.81
## [32,]  492838.02
## [33,]  622253.17
## [34,]  283873.02
## [35,]  361509.56

SEeblup=sqrt(eblup$mse) #standar error EBLUP
mean.eblup=mean(SEeblup) #Rataan SE EBLUP
SEdiret=pengeluaran[,2] # Standar error penduga langsung
mean.SEdirect=mean(SEdiret) #rataan SE penduga langsg

#Plot direct estimasi dan EBLUP
x=c(1:35)
#line SEeblup
plot(x, SEeblup, type="o", col="blue", pch="o", ylab="Standar Error", lty=1)
#line SEdiret
points(x, SEdiret, col="red", pch="*")
lines(x, SEdiret, col="red",lty=2)
#membuat legenda
legend(1,150000,legend=c("SEeblup","SEdirect"), col=c("blue","red"),
       pch=c("o","*"),lty=c(1,2), ncol=1)

Referensi

Battese, G. E., Harter, R. M., & Fuller, W. A. (1988). An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83(401), 28–36. https://doi.org/10.1080/01621459.1988.10478561
J.N.K R, Molina I. 2015. Small Area Estimation. Second Edi. Hoboken, New Jersey: John Wiley & Sons, Inc. https://www.wiley.com/en-us/Small+Area+Estimation%2C+2nd+Edition-p-9781118735787

Empirical Best Linear Unbiassed Prediction (EBLUP)

in Maximum Likelihood Estimation

Achmad Fauzan

12 April 2025

1 Motivasi EBLUP

2 Contoh 1

3 Contoh 2

No	Variabel	Keterangan
1	Y	Rata-rata pengeluaran rumah tangga per kapita per bulan di desa se kabupaten Jember
2	\(X_1\)	Persentase keluarga yang menggantungkan hidupnya pada pertanian
3	\(X_2\)	Jumlah keluarga yang menerima askeskin dalam satu tahun terakhir
4	\(X_3\)	Jumlah keluarga pengguna listrik PLN
5	\(X_4\)	Jumlah keluarga yang mengenyam pendidikan SD, SMP, SMA, dan PT
6	\(X_5\)	Jumlah keluarga yang tinggal di pemukiman kumuh
7	\(X_6\)	Jumlah keluarga yang memiliki Surat Keterangan Tidak Mampu (SKTM) dan satu tahun terakhir
8	\(X_7\)	Jumlah keluarga yang pernah belajar di lembaga pendidikan dan keterampilan
9	\(X_8\)	Jumlah keluarga yang memiliki anggota keluarga sebagai Tenaga Kerja Indonesia