library(sae)## Loading required package: MASS## Loading required package: lme4## Loading required package: Matrixlibrary(lme4)
library(knitr)
library(kableExtra)# Load data set for segments (units within domains)
data(cornsoybean)
cornsoybean## County CornHec SoyBeansHec CornPix SoyBeansPix
## 1 1 165.76 8.09 374 55
## 2 2 96.32 106.03 209 218
## 3 3 76.08 103.60 253 250
## 4 4 185.35 6.47 432 96
## 5 4 116.43 63.82 367 178
## 6 5 162.08 43.50 361 137
## 7 5 152.04 71.43 288 206
## 8 5 161.75 42.49 369 165
## 9 6 92.88 105.26 206 218
## 10 6 149.94 76.49 316 221
## 11 6 64.75 174.34 145 338
## 12 7 127.07 95.67 355 128
## 13 7 133.55 76.57 295 147
## 14 7 77.70 93.48 223 204
## 15 8 206.39 37.84 459 77
## 16 8 108.33 131.12 290 217
## 17 8 118.17 124.44 307 258
## 18 9 99.96 144.15 252 303
## 19 9 140.43 103.60 293 221
## 20 9 98.95 88.59 206 222
## 21 9 131.04 115.58 302 274
## 22 10 114.12 99.15 313 190
## 23 10 100.60 124.56 246 270
## 24 10 127.88 110.88 353 172
## 25 10 116.90 109.14 271 228
## 26 10 87.41 143.66 237 297
## 27 11 93.48 91.05 221 167
## 28 11 121.00 132.33 369 191
## 29 11 109.91 143.14 343 249
## 30 11 122.66 104.13 342 182
## 31 11 104.21 118.57 294 179
## 32 12 88.59 102.59 220 262
## 33 12 88.59 29.46 340 87
## 34 12 165.35 69.28 355 160
## 35 12 104.00 99.15 261 221
## 36 12 88.63 143.66 187 345
## 37 12 153.70 94.49 350 190# Load data set for counties
data(cornsoybeanmeans)
cornsoybeanmeans## CountyIndex CountyName SampSegments PopnSegments MeanCornPixPerSeg
## 1 1 CerroGordo 1 545 295.29
## 2 2 Hamilton 1 566 300.40
## 3 3 Worth 1 394 289.60
## 4 4 Humboldt 2 424 290.74
## 5 5 Franklin 3 564 318.21
## 6 6 Pocahontas 3 570 257.17
## 7 7 Winnebago 3 402 291.77
## 8 8 Wright 3 567 301.26
## 9 9 Webster 4 687 262.17
## 10 10 Hancock 5 569 314.28
## 11 11 Kossuth 5 965 298.65
## 12 12 Hardin 6 556 325.99
## MeanSoyBeansPixPerSeg
## 1 189.70
## 2 196.65
## 3 205.28
## 4 220.22
## 5 188.06
## 6 247.13
## 7 185.37
## 8 221.36
## 9 247.09
## 10 198.66
## 11 204.61
## 12 177.05attach(cornsoybeanmeans)
# Construct data frame with county means of auxiliary variables for domains. First column must include the county code
Xmean <- data.frame(CountyIndex, MeanCornPixPerSeg, MeanSoyBeansPixPerSeg)
Xmean## CountyIndex MeanCornPixPerSeg MeanSoyBeansPixPerSeg
## 1 1 295.29 189.70
## 2 2 300.40 196.65
## 3 3 289.60 205.28
## 4 4 290.74 220.22
## 5 5 318.21 188.06
## 6 6 257.17 247.13
## 7 7 291.77 185.37
## 8 8 301.26 221.36
## 9 9 262.17 247.09
## 10 10 314.28 198.66
## 11 11 298.65 204.61
## 12 12 325.99 177.05Popn <- data.frame(CountyIndex, PopnSegments)
Popn## CountyIndex PopnSegments
## 1 1 545
## 2 2 566
## 3 3 394
## 4 4 424
## 5 5 564
## 6 6 570
## 7 7 402
## 8 8 567
## 9 9 687
## 10 10 569
## 11 11 965
## 12 12 556resultCorn <- eblupBHF(CornHec ~ CornPix + SoyBeansPix,
dom=County,
meanxpop=Xmean,
popnsize=Popn,
data=cornsoybean)
resultCorn$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 6resultCorn$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
## Xs(Intercept) 17.96398 30.97450 0.580
## XsCornPix 0.36634 0.06496 5.640
## XsSoyBeansPix -0.03036 0.06758 -0.449
##
## Correlation of Fixed Effects:
## Xs(In) XsCrnP
## XsCornPix -0.947
## XsSoyBensPx -0.897 0.734
##
## $fixed
## Xs(Intercept) XsCornPix XsSoyBeansPix
## 17.9639791 0.3663352 -0.0303638
##
## $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
##
## $errorvar
## [1] 297.7128
##
## $refvar
## [1] 63.3149
##
## $loglike
## [1] -161.0058
##
## $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.0396194detach(cornsoybeanmeans)
set.seed(123)
BHFCorn <- pbmseBHF(CornHec ~ CornPix + SoyBeansPix,
dom = County,
meanxpop = Xmean,
popnsize = Popn,
B = 200,
data = cornsoybean)##
## Bootstrap procedure with B = 200 iterations starts.
## b = 1
## b = 2
## b = 3
## b = 4
## b = 5
## b = 6
## b = 7## boundary (singular) fit: see help('isSingular')## b = 8
## b = 9
## b = 10## boundary (singular) fit: see help('isSingular')## b = 11
## b = 12
## b = 13
## b = 14
## b = 15
## b = 16
## b = 17## boundary (singular) fit: see help('isSingular')## b = 18
## b = 19## boundary (singular) fit: see help('isSingular')## b = 20
## b = 21
## b = 22
## b = 23
## b = 24
## b = 25
## b = 26
## b = 27
## b = 28
## b = 29## boundary (singular) fit: see help('isSingular')## b = 30
## b = 31
## b = 32## boundary (singular) fit: see help('isSingular')## b = 33## boundary (singular) fit: see help('isSingular')## b = 34
## b = 35
## 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
## b = 46
## b = 47## boundary (singular) fit: see help('isSingular')## b = 48
## b = 49
## b = 50
## b = 51
## b = 52
## b = 53
## b = 54## boundary (singular) fit: see help('isSingular')## b = 55
## b = 56
## b = 57
## b = 58
## b = 59## boundary (singular) fit: see help('isSingular')## b = 60
## b = 61
## b = 62
## b = 63
## b = 64
## b = 65## boundary (singular) fit: see help('isSingular')## b = 66
## b = 67
## b = 68
## b = 69
## b = 70
## b = 71
## b = 72
## b = 73
## b = 74
## b = 75## boundary (singular) fit: see help('isSingular')## b = 76
## b = 77## boundary (singular) fit: see help('isSingular')## b = 78
## b = 79## boundary (singular) fit: see help('isSingular')## b = 80
## b = 81## boundary (singular) fit: see help('isSingular')## b = 82
## b = 83
## b = 84
## b = 85## boundary (singular) fit: see help('isSingular')## b = 86
## b = 87
## b = 88
## b = 89
## b = 90
## b = 91
## b = 92## boundary (singular) fit: see help('isSingular')## b = 93## boundary (singular) fit: see help('isSingular')## b = 94
## b = 95## boundary (singular) fit: see help('isSingular')## b = 96## boundary (singular) fit: see help('isSingular')## b = 97## boundary (singular) fit: see help('isSingular')## b = 98
## b = 99
## b = 100
## b = 101
## b = 102
## b = 103
## b = 104
## b = 105## boundary (singular) fit: see help('isSingular')## b = 106## boundary (singular) fit: see help('isSingular')## b = 107## boundary (singular) fit: see help('isSingular')## b = 108## boundary (singular) fit: see help('isSingular')## b = 109## boundary (singular) fit: see help('isSingular')## b = 110
## b = 111
## b = 112
## b = 113## boundary (singular) fit: see help('isSingular')## b = 114
## b = 115## boundary (singular) fit: see help('isSingular')## b = 116
## b = 117
## b = 118
## b = 119
## b = 120## boundary (singular) fit: see help('isSingular')## b = 121
## b = 122
## b = 123
## b = 124
## b = 125## boundary (singular) fit: see help('isSingular')## b = 126## boundary (singular) fit: see help('isSingular')## b = 127
## b = 128
## b = 129
## b = 130
## b = 131
## b = 132
## b = 133
## b = 134
## b = 135
## b = 136
## b = 137## boundary (singular) fit: see help('isSingular')## b = 138
## b = 139## boundary (singular) fit: see help('isSingular')## b = 140
## b = 141
## b = 142
## b = 143
## b = 144## boundary (singular) fit: see help('isSingular')## b = 145## boundary (singular) fit: see help('isSingular')## b = 146## boundary (singular) fit: see help('isSingular')## b = 147
## b = 148
## b = 149
## b = 150
## b = 151
## b = 152
## b = 153
## b = 154
## b = 155
## b = 156## boundary (singular) fit: see help('isSingular')## b = 157
## b = 158
## b = 159
## b = 160## boundary (singular) fit: see help('isSingular')## b = 161## boundary (singular) fit: see help('isSingular')## b = 162
## b = 163
## b = 164
## b = 165
## b = 166
## b = 167
## b = 168## boundary (singular) fit: see help('isSingular')## b = 169
## b = 170
## b = 171
## b = 172
## b = 173
## b = 174
## b = 175## boundary (singular) fit: see help('isSingular')## b = 176## boundary (singular) fit: see help('isSingular')## b = 177
## b = 178## boundary (singular) fit: see help('isSingular')## b = 179
## b = 180
## b = 181
## b = 182
## b = 183
## b = 184
## b = 185## boundary (singular) fit: see help('isSingular')## b = 186## boundary (singular) fit: see help('isSingular')## b = 187
## b = 188
## b = 189
## b = 190## boundary (singular) fit: see help('isSingular')## b = 191## boundary (singular) fit: see help('isSingular')## b = 192
## b = 193
## b = 194
## b = 195
## b = 196## boundary (singular) fit: see help('isSingular')## b = 197
## b = 198
## b = 199
## b = 200mseCorn <- BHFCorn$mse
mseCorn## domain mse
## 1 1 85.00324
## 2 2 74.65850
## 3 3 74.88418
## 4 4 64.32547
## 5 5 69.21492
## 6 6 50.86237
## 7 7 46.95725
## 8 8 49.40878
## 9 9 46.06986
## 10 10 39.00881
## 11 11 37.91072
## 12 12 35.47207residCorn <- BHFCorn$est$fit$residuals
cat("\nResidual dari Corn :\n")##
## Residual dari Corn :residCorn## [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.0396194sigmav2estCorn <- BHFCorn$est$fit$refvar
cat("\nEstimasi Sigma(v)^2 dari Corn =",sigmav2estCorn)##
## Estimasi Sigma(v)^2 dari Corn = 63.3149sigmae2estCorn <- BHFCorn$est$fit$errorvar
cat("\nEstimasi Sigma(e)^2 dari Corn =",sigmae2estCorn)##
## Estimasi Sigma(e)^2 dari Corn = 297.7128## ============================================================
## Persiapan Data
## ============================================================
attach(cornsoybeanmeans)
Xmean <- data.frame(CountyIndex, MeanCornPixPerSeg, MeanSoyBeansPixPerSeg)
Xmean## CountyIndex MeanCornPixPerSeg MeanSoyBeansPixPerSeg
## 1 1 295.29 189.70
## 2 2 300.40 196.65
## 3 3 289.60 205.28
## 4 4 290.74 220.22
## 5 5 318.21 188.06
## 6 6 257.17 247.13
## 7 7 291.77 185.37
## 8 8 301.26 221.36
## 9 9 262.17 247.09
## 10 10 314.28 198.66
## 11 11 298.65 204.61
## 12 12 325.99 177.05Popn <- data.frame(CountyIndex, PopnSegments)
Popn## CountyIndex PopnSegments
## 1 1 545
## 2 2 566
## 3 3 394
## 4 4 424
## 5 5 564
## 6 6 570
## 7 7 402
## 8 8 567
## 9 9 687
## 10 10 569
## 11 11 965
## 12 12 556detach(cornsoybeanmeans)
# Variabel dari data sampel
y <- cornsoybean$CornHec
X <- cbind(1, cornsoybean$CornPix, cornsoybean$SoyBeansPix)
dom <- cornsoybean$County
# Informasi area
areas <- sort(unique(dom))
m <- length(areas)
ni <- as.vector(table(dom))
names(ni) <- areas
cat("Jumlah sampel per area (ni):\n")## Jumlah sampel per area (ni):print(ni)## 1 2 3 4 5 6 7 8 9 10 11 12
## 1 1 1 2 3 3 3 3 4 5 5 6cat("Total sampel:", sum(ni), "\n\n")## Total sampel: 37# Rata-rata X populasi
Xbar_pop <- Xmean[order(Xmean$CountyIndex), ]
# Rata-rata X dan Y sampel per area
ybar_ia <- tapply(y, dom, mean)
Xbar_ia <- aggregate(
data.frame(CornPix = cornsoybean$CornPix,
SoyBeansPix = cornsoybean$SoyBeansPix),
by = list(dom = dom),
FUN = mean
)
## ============================================================
## LANGKAH 1: Estimasi Komponen Varians (sigma_v^2, sigma_e^2)
## Menggunakan REML via lme4
## Persamaan (10) dari slide: Vi = sigma_v^2 * 1_ni*1_ni' + sigma_e^2 * I_ni
## ============================================================
data_model <- data.frame(
y = y,
Xs1 = cornsoybean$CornPix,
Xs2 = cornsoybean$SoyBeansPix,
dom = factor(dom)
)
fit_lme <- lmer(y ~ Xs1 + Xs2 + (1 | dom),
data = data_model,
REML = TRUE)
summary(fit_lme)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ Xs1 + Xs2 + (1 | dom)
## Data: data_model
##
## 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
## (Intercept) 17.96398 30.97450 0.580
## Xs1 0.36634 0.06496 5.640
## Xs2 -0.03036 0.06758 -0.449
##
## Correlation of Fixed Effects:
## (Intr) Xs1
## Xs1 -0.947
## Xs2 -0.897 0.734# Ekstrak komponen varians
vc <- as.data.frame(VarCorr(fit_lme))
sigma_v2 <- vc$vcov[vc$grp == "dom"] # sigma_v^2 (refvar)
sigma_e2 <- vc$vcov[vc$grp == "Residual"] # sigma_e^2 (errorvar)
cat("=== Estimasi Komponen Varians ===\n")## === Estimasi Komponen Varians ===cat("Estimasi Sigma(v)^2 =", sigma_v2, "\n")## Estimasi Sigma(v)^2 = 63.3149cat("Estimasi Sigma(e)^2 =", sigma_e2, "\n\n")## Estimasi Sigma(e)^2 = 297.7128## ============================================================
## LANGKAH 2: Estimasi gamma_hat_i (bobot shrinkage)
## Persamaan dari slide:
## gamma_hat_i = sigma_v^2 / (sigma_v^2 + sigma_e^2 / a_i.)
## a_ij = k_ij^(-2) = 1 (karena k_ij = 1)
## a_i. = sum_j(a_ij) = n_i
## ============================================================
a_i <- ni # a_i. = n_i karena k_ij = 1
gamma_hat <- sigma_v2 / (sigma_v2 + sigma_e2 / a_i)
cat("=== Gamma_hat per Area ===\n")## === Gamma_hat per Area ===print(round(data.frame(
Area = areas,
ni = ni,
gamma_hat = gamma_hat
), 6))## Area ni gamma_hat
## 1 1 1 0.175374
## 2 2 1 0.175374
## 3 3 1 0.175374
## 4 4 2 0.298414
## 5 5 3 0.389504
## 6 6 3 0.389504
## 7 7 3 0.389504
## 8 8 3 0.389504
## 9 9 4 0.459659
## 10 10 5 0.515352
## 11 11 5 0.515352
## 12 12 6 0.560638cat("\n")## ============================================================
## LANGKAH 3: Estimasi beta_hat (Fixed Effects / BLUE)
## Persamaan (13): beta_tilde = (sum X_i' V_i^-1 X_i)^-1 (sum X_i' V_i^-1 y_i)
## Diambil langsung dari hasil lmer
## ============================================================
beta_hat <- as.vector(fixef(fit_lme))
names(beta_hat) <- c("Intercept", "CornPix", "SoyBeansPix")
cat("=== Beta_hat (Fixed Effects) ===\n")## === Beta_hat (Fixed Effects) ===print(beta_hat)## Intercept CornPix SoyBeansPix
## 17.9639791 0.3663352 -0.0303638cat("\n")## ============================================================
## LANGKAH 4: Hitung EBLUP untuk mu_i
## Persamaan (17) dari slide:
## mu_hat_i = gamma_hat_i * [ybar_ia + (Xbar_i - xbar_ia)' beta_hat]
## + (1 - gamma_hat_i) * Xbar_i' * beta_hat
## ============================================================
eblup_corn_manual <- numeric(m)
names(eblup_corn_manual) <- areas
for (i in seq_along(areas)) {
area_id <- areas[i]
# ybar_ia: rata-rata y sampel area ke-i
ybar_i <- ybar_ia[as.character(area_id)]
# xbar_ia: rata-rata X sampel area ke-i
row_ia <- Xbar_ia[Xbar_ia$dom == area_id, ]
xbar_ia_i <- c(row_ia$CornPix, row_ia$SoyBeansPix)
# Xbar_i: rata-rata X populasi area ke-i
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
# Regresi sintetik: Xbar_i' * beta_hat
synth_i <- beta_hat["Intercept"] + sum(Xbar_i * beta_hat[-1])
# Regresi survey: ybar_ia + (Xbar_i - xbar_ia)' * beta_hat
survey_i <- ybar_i + sum((Xbar_i - xbar_ia_i) * beta_hat[-1])
# EBLUP (Pers. 17)
eblup_corn_manual[i] <- gamma_hat[i] * survey_i +
(1 - gamma_hat[i]) * synth_i
}
cat("=== EBLUP Manual untuk Corn (All Counties) ===\n")## === EBLUP Manual untuk Corn (All Counties) ===print(data.frame(
domain = areas,
eblup = round(eblup_corn_manual, 4),
sampsize = ni
))## domain eblup sampsize
## 1 1 122.5637 1
## 2 2 123.5152 1
## 3 3 113.0907 1
## 4 4 115.0207 2
## 5 5 137.1962 3
## 6 6 108.9454 3
## 7 7 116.5155 3
## 8 8 122.7615 3
## 9 9 111.5303 4
## 10 10 124.1803 5
## 11 11 112.5047 5
## 12 12 131.2579 6## ============================================================
## LANGKAH 5: Hitung Random Effects v_hat_i
## v_hat_i = gamma_i * (ybar_ia - xbar_ia' * beta_hat)
## ============================================================
v_hat <- numeric(m)
names(v_hat) <- areas
for (i in seq_along(areas)) {
area_id <- areas[i]
ybar_i <- ybar_ia[as.character(area_id)]
row_ia <- Xbar_ia[Xbar_ia$dom == area_id, ]
xbar_ia_i <- c(row_ia$CornPix, row_ia$SoyBeansPix)
xbar_pred <- beta_hat["Intercept"] + sum(xbar_ia_i * beta_hat[-1])
v_hat[i] <- gamma_hat[i] * (ybar_i - xbar_pred)
}
cat("\n=== Random Effects Manual (v_hat_i) ===\n")##
## === Random Effects Manual (v_hat_i) ===print(round(v_hat, 6))## 1 2 3 4 5 6 7 8
## 2.184574 1.475118 -4.730863 -2.764825 8.370915 4.274827 -2.705540 1.156682
## 9 10 11 12
## 5.026852 -2.883398 -8.652532 -0.751808cat("\n=== Random Effects dari lme4 (pembanding) ===\n")##
## === Random Effects dari lme4 (pembanding) ===print(round(as.numeric(ranef(fit_lme)$dom[, 1]), 6))## [1] 2.184574 1.475118 -4.730863 -2.764825 8.370915 4.274827 -2.705540
## [8] 1.156682 5.026852 -2.883398 -8.652532 -0.751808## ============================================================
## LANGKAH 6: Hitung Residual Manual
## resid_ij = y_ij - (X_ij * beta_hat + v_hat_i)
## ============================================================
v_hat_ij <- rep(v_hat, times = ni)
fitted_manual <- as.numeric(X %*% beta_hat) + v_hat_ij
resid_manual <- y - fitted_manual
cat("\n=== Residual Manual ===\n")##
## === Residual Manual ===print(round(resid_manual, 6))## 1 2 3 4 4 5 5
## 10.272080 6.936148 -22.244981 14.808951 -27.809428 7.657928 26.455502
## 5 6 6 6 7 7 7
## 5.247432 1.795444 18.649660 -0.344451 -14.350880 14.686146 -13.056981
## 8 8 8 9 9 9 9
## 21.459481 -10.438934 -5.581717 -6.147079 16.813345 7.234874 5.735610
## 10 10 10 10 10 11 11
## -9.854387 3.599177 -11.294344 9.465517 -5.473983 8.279221 -17.689662
## 11 11 11 12 12 12 12
## -17.493846 -6.411885 -7.368885 -1.260607 -50.534499 22.947029 -2.115267
## 12 12
## 13.388650 14.039619cat("\n=== Residual dari lme4 (pembanding) ===\n")##
## === Residual dari lme4 (pembanding) ===print(round(as.numeric(residuals(fit_lme)), 6))## [1] 10.272080 6.936148 -22.244981 14.808951 -27.809428 7.657928
## [7] 26.455502 5.247432 1.795444 18.649660 -0.344451 -14.350880
## [13] 14.686146 -13.056981 21.459481 -10.438934 -5.581717 -6.147079
## [19] 16.813345 7.234874 5.735610 -9.854387 3.599177 -11.294344
## [25] 9.465517 -5.473983 8.279221 -17.689662 -17.493846 -6.411885
## [31] -7.368885 -1.260607 -50.534499 22.947029 -2.115267 13.388650
## [37] 14.039619## ============================================================
## LANGKAH 8: Bootstrap MSE Manual (B = 200)
## Mengikuti 8 langkah dari slide Pendugaan MSE
## ============================================================
set.seed(123)
B <- 200
mse_boot <- numeric(m)
names(mse_boot) <- areas
for (b in 1:B) {
# Langkah 2: Bangkitkan v_i* ~ N(0, sigma_v2) dan e_ij* ~ N(0, sigma_e2)
v_star <- rnorm(m, 0, sqrt(sigma_v2))
e_star <- rnorm(sum(ni), 0, sqrt(sigma_e2))
# Langkah 3: Data bootstrap y_ij* = X_ij * beta_hat + v_i* + e_ij*
v_star_ij <- rep(v_star, times = ni)
y_star <- as.numeric(X %*% beta_hat) + v_star_ij + e_star
# Langkah 4: Parameter bootstrap mu_i*(b) = Xbar_i' * beta_hat + v_i*
mu_star <- numeric(m)
for (i in seq_along(areas)) {
area_id <- areas[i]
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
mu_star[i] <- beta_hat["Intercept"] +
sum(Xbar_i * beta_hat[-1]) +
v_star[i]
}
# Langkah 5: Estimasi ulang parameter dari data bootstrap
data_boot <- data.frame(
y_star = y_star,
Xs1 = cornsoybean$CornPix,
Xs2 = cornsoybean$SoyBeansPix,
dom = factor(dom)
)
fit_boot <- tryCatch(
lmer(y_star ~ Xs1 + Xs2 + (1 | dom),
data = data_boot, REML = TRUE),
error = function(e) NULL
)
if (is.null(fit_boot)) next
vc_boot <- as.data.frame(VarCorr(fit_boot))
sigma_v2_b <- max(vc_boot$vcov[vc_boot$grp == "dom"], 0)
sigma_e2_b <- vc_boot$vcov[vc_boot$grp == "Residual"]
beta_hat_b <- as.vector(fixef(fit_boot))
gamma_b <- sigma_v2_b / (sigma_v2_b + sigma_e2_b / a_i)
# Rata-rata y* dan X sampel per area
ybar_boot <- tapply(y_star, dom, mean)
Xbar_ia_b <- aggregate(
data.frame(CornPix = cornsoybean$CornPix,
SoyBeansPix = cornsoybean$SoyBeansPix),
by = list(dom = dom),
FUN = mean
)
# Langkah 6: EBLUP bootstrap mu_hat_i*(b)
mu_hat_star <- numeric(m)
for (i in seq_along(areas)) {
area_id <- areas[i]
ybar_b_i <- ybar_boot[as.character(area_id)]
row_ia_b <- Xbar_ia_b[Xbar_ia_b$dom == area_id, ]
xbar_ia_b <- c(row_ia_b$CornPix, row_ia_b$SoyBeansPix)
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
synth_b <- beta_hat_b[1] + sum(Xbar_i * beta_hat_b[-1])
survey_b <- ybar_b_i + sum((Xbar_i - xbar_ia_b) * beta_hat_b[-1])
mu_hat_star[i] <- gamma_b[i] * survey_b +
(1 - gamma_b[i]) * synth_b
}
# Langkah 8: Akumulasi MSE
mse_boot <- mse_boot + (mu_hat_star - mu_star)^2
}## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')# Rata-rata MSE bootstrap
mse_boot <- mse_boot / B
RMSE <- sqrt(mse_boot)
## ============================================================
## LANGKAH 9: Ringkasan Hasil MSE Manual
## ============================================================
mseCorn_manual <- data.frame(
domain = areas,
eblup = round(eblup_corn_manual, 4),
mse = round(mse_boot, 4),
rmse = round(RMSE, 4),
cv_pct = round(100 * RMSE / eblup_corn_manual, 2)
)
cat("\n=== MSE Bootstrap Manual (B=200) ===\n")##
## === MSE Bootstrap Manual (B=200) ===print(mseCorn_manual)## domain eblup mse rmse cv_pct
## 1 1 122.5637 73.7317 8.5867 7.01
## 2 2 123.5152 82.1643 9.0645 7.34
## 3 3 113.0907 79.8922 8.9382 7.90
## 4 4 115.0207 62.0760 7.8788 6.85
## 5 5 137.1962 61.8099 7.8619 5.73
## 6 6 108.9454 50.9765 7.1398 6.55
## 7 7 116.5155 57.7378 7.5985 6.52
## 8 8 122.7615 56.9339 7.5455 6.15
## 9 9 111.5303 41.1925 6.4181 5.75
## 10 10 124.1803 44.9904 6.7075 5.40
## 11 11 112.5047 42.3023 6.5040 5.78
## 12 12 131.2579 45.6489 6.7564 5.15# Residual, sigma_v^2, sigma_e^2 (konsisten dengan alur awal)
residCorn_manual <- resid_manual
sigmav2estCorn_manual <- sigma_v2
sigmae2estCorn_manual <- sigma_e2
cat("\nResidual Manual:\n")##
## Residual Manual:print(round(residCorn_manual, 6))## 1 2 3 4 4 5 5
## 10.272080 6.936148 -22.244981 14.808951 -27.809428 7.657928 26.455502
## 5 6 6 6 7 7 7
## 5.247432 1.795444 18.649660 -0.344451 -14.350880 14.686146 -13.056981
## 8 8 8 9 9 9 9
## 21.459481 -10.438934 -5.581717 -6.147079 16.813345 7.234874 5.735610
## 10 10 10 10 10 11 11
## -9.854387 3.599177 -11.294344 9.465517 -5.473983 8.279221 -17.689662
## 11 11 11 12 12 12 12
## -17.493846 -6.411885 -7.368885 -1.260607 -50.534499 22.947029 -2.115267
## 12 12
## 13.388650 14.039619cat("\nEstimasi Sigma(v)^2 =", sigmav2estCorn_manual)##
## Estimasi Sigma(v)^2 = 63.3149cat("\nEstimasi Sigma(e)^2 =", sigmae2estCorn_manual, "\n")##
## Estimasi Sigma(e)^2 = 297.7128## ============================================================
## LANGKAH 10: Verifikasi dengan pbmseBHF()
## ============================================================
cat("\n--Perbandingan MSE Manual vs pbmseBHF() dari Estimasi Corn--\n")##
## --Perbandingan MSE Manual vs pbmseBHF() dari Estimasi Corn--compare <- data.frame(
Domain = areas,
EBLUP_Manual = round(eblup_corn_manual, 4),
EBLUP_Fungsi = round(resultCorn$eblup$eblup, 4),
MSE_Manual = round(mse_boot, 4),
MSE_Fungsi = round(BHFCorn$mse$mse, 4)
)
kable(compare,
caption = "Perbandingan EBLUP dan MSE Manual vs pbmseBHF() dari Estimasi Corn",
col.names = c("Domain", "EBLUP Manual", "EBLUP Fungsi",
"MSE Manual", "MSE Fungsi"),
align = c("c","r","r","r","r"),
booktabs = TRUE) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "bordered", "condensed"),
full_width = FALSE,
position = "center"
) %>%
row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")| Domain | EBLUP Manual | EBLUP Fungsi | MSE Manual | MSE Fungsi |
|---|---|---|---|---|
| 1 | 122.5637 | 122.5825 | 73.7317 | 85.0032 |
| 2 | 123.5152 | 123.5274 | 82.1643 | 74.6585 |
| 3 | 113.0907 | 113.0343 | 79.8922 | 74.8842 |
| 4 | 115.0207 | 114.9901 | 62.0760 | 64.3255 |
| 5 | 137.1962 | 137.2660 | 61.8099 | 69.2149 |
| 6 | 108.9454 | 108.9807 | 50.9765 | 50.8624 |
| 7 | 116.5155 | 116.4839 | 57.7378 | 46.9572 |
| 8 | 122.7615 | 122.7711 | 56.9339 | 49.4088 |
| 9 | 111.5303 | 111.5648 | 41.1925 | 46.0699 |
| 10 | 124.1803 | 124.1565 | 44.9904 | 39.0088 |
| 11 | 112.5047 | 112.4626 | 42.3023 | 37.9107 |
| 12 | 131.2579 | 131.2515 | 45.6489 | 35.4721 |
attach(cornsoybeanmeans)
# Compute EBLUPs of county means of SOY BEANS crop areas for a SUBSET of COUNTIES using ML method
domains <- c(1, 5, 10)
resultBean <- eblupBHF(SoyBeansHec ~ CornPix + SoyBeansPix,
dom=County,
selectdom=domains,
meanxpop=Xmean,
popnsize=Popn,
method="ML",
data=cornsoybean)
resultBean$eblup## domain eblup sampsize
## 1 1 78.63690 1
## 2 5 66.26077 3
## 3 10 100.61124 5resultBean$fit## $summary
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: ys ~ -1 + Xs + (1 | dom)
##
## AIC BIC logLik -2*log(L) df.resid
## 323.2 331.2 -156.6 313.2 32
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9871 -0.5487 -0.1716 0.4525 1.8864
##
## Random effects:
## Groups Name Variance Std.Dev.
## dom (Intercept) 219.3 14.81
## Residual 170.3 13.05
## Number of obs: 37, groups: dom, 12
##
## Fixed effects:
## Estimate Std. Error t value
## Xs(Intercept) -16.34566 25.46384 -0.642
## XsCornPix 0.02807 0.05348 0.525
## XsSoyBeansPix 0.49679 0.05537 8.972
##
## Correlation of Fixed Effects:
## Xs(In) XsCrnP
## XsCornPix -0.935
## XsSoyBensPx -0.878 0.718
##
## $fixed
## Xs(Intercept) XsCornPix XsSoyBeansPix
## -16.34565637 0.02806513 0.49678511
##
## $random
## (Intercept)
## 1 -7.5341959
## 2 4.6221614
## 3 -6.3898792
## 4 -20.0120589
## 5 -19.7224359
## 6 0.1073773
## 7 13.8480550
## 8 10.2146568
## 9 -3.9717549
## 10 9.4323906
## 11 25.2796992
## 12 -5.8740155
##
## $errorvar
## [1] 170.2856
##
## $refvar
## [1] 219.322
##
## $loglike
## [1] -156.5848
##
## $residuals
## [1] -5.8496865 3.5887301 -4.9612183 -16.9877903 1.4500644 1.3770215
## [7] -2.9223964 -13.7674825 7.4177096 -25.9298098 18.5954698 24.6159873
## [13] -2.2390220 -11.6250838 -7.1633474 21.3097445 -6.2155520 6.8691120
## [19] 5.9048204 -7.1602986 -8.6973764 2.8897105 -9.5627343 22.4392373
## [25] -4.8193881 -3.6233460 -7.0495488 18.1539697 0.8801269 -4.8172059
## [31] 12.4602756 -11.5223541 -1.0827760 2.0509344 4.2551650 -10.7593686
## [37] 12.4977069detach(cornsoybeanmeans)
set.seed(123)
BHFBean <- pbmseBHF(SoyBeansHec ~ CornPix + SoyBeansPix,
dom = County,
meanxpop = Xmean,
popnsize = Popn,
B = 200,
method="ML",
data = cornsoybean)##
## Bootstrap procedure with B = 200 iterations starts.
## b = 1
## b = 2
## b = 3
## b = 4
## b = 5
## b = 6
## b = 7## boundary (singular) fit: see help('isSingular')## b = 8
## b = 9
## b = 10
## b = 11
## b = 12
## b = 13
## b = 14
## b = 15
## b = 16
## b = 17## boundary (singular) fit: see help('isSingular')## b = 18
## b = 19
## b = 20
## b = 21
## b = 22
## b = 23
## b = 24
## b = 25
## b = 26
## b = 27
## b = 28
## b = 29
## b = 30
## b = 31
## b = 32
## b = 33
## b = 34
## b = 35
## 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
## b = 46
## b = 47
## b = 48
## b = 49
## b = 50
## b = 51
## b = 52
## b = 53
## b = 54
## b = 55
## b = 56
## b = 57
## b = 58
## b = 59
## b = 60
## b = 61
## b = 62
## b = 63
## b = 64
## b = 65
## b = 66
## b = 67
## b = 68
## b = 69
## b = 70
## b = 71
## b = 72
## b = 73
## b = 74
## b = 75
## b = 76
## b = 77
## b = 78
## b = 79
## b = 80
## b = 81
## b = 82
## b = 83
## b = 84
## b = 85
## b = 86
## b = 87
## b = 88
## b = 89
## b = 90
## b = 91
## b = 92
## b = 93
## b = 94
## b = 95
## b = 96
## b = 97
## b = 98
## b = 99
## b = 100
## b = 101
## b = 102
## b = 103
## b = 104
## b = 105
## b = 106
## b = 107
## b = 108## boundary (singular) fit: see help('isSingular')## b = 109
## b = 110
## b = 111
## b = 112
## b = 113
## b = 114
## b = 115
## b = 116
## b = 117
## b = 118
## b = 119
## b = 120
## b = 121
## b = 122
## b = 123
## b = 124
## b = 125
## b = 126## boundary (singular) fit: see help('isSingular')## b = 127
## b = 128
## b = 129
## b = 130
## b = 131
## b = 132
## b = 133
## b = 134
## b = 135
## b = 136
## b = 137
## b = 138
## b = 139
## b = 140
## b = 141
## b = 142
## b = 143
## b = 144
## b = 145
## b = 146
## b = 147
## b = 148
## b = 149
## b = 150
## b = 151
## b = 152
## b = 153
## b = 154
## b = 155
## b = 156## boundary (singular) fit: see help('isSingular')## b = 157
## b = 158
## b = 159
## b = 160
## b = 161
## b = 162
## b = 163
## b = 164
## b = 165
## b = 166
## b = 167
## b = 168
## b = 169
## b = 170
## b = 171
## b = 172
## b = 173
## b = 174
## b = 175
## b = 176
## b = 177
## b = 178
## b = 179
## b = 180
## b = 181
## b = 182
## b = 183
## b = 184
## b = 185
## b = 186
## b = 187
## b = 188
## b = 189
## b = 190
## b = 191
## b = 192
## b = 193
## b = 194
## b = 195
## b = 196
## b = 197
## b = 198
## b = 199
## b = 200mseBean <- BHFBean$mse
mseBean## domain mse
## 1 1 145.03062
## 2 2 122.08925
## 3 3 115.27881
## 4 4 78.08470
## 5 5 64.62653
## 6 6 47.72946
## 7 7 46.94591
## 8 8 51.73991
## 9 9 37.32158
## 10 10 30.01241
## 11 11 30.28057
## 12 12 28.85243residBean <- BHFBean$est$fit$residuals
cat("\nResidual dari Bean untuk All Counties :\n")##
## Residual dari Bean untuk All Counties :residBean## [1] -5.8496865 3.5887301 -4.9612183 -16.9877903 1.4500644 1.3770215
## [7] -2.9223964 -13.7674825 7.4177096 -25.9298098 18.5954698 24.6159873
## [13] -2.2390220 -11.6250838 -7.1633474 21.3097445 -6.2155520 6.8691120
## [19] 5.9048204 -7.1602986 -8.6973764 2.8897105 -9.5627343 22.4392373
## [25] -4.8193881 -3.6233460 -7.0495488 18.1539697 0.8801269 -4.8172059
## [31] 12.4602756 -11.5223541 -1.0827760 2.0509344 4.2551650 -10.7593686
## [37] 12.4977069sigmav2estBean <- BHFBean$est$fit$refvar
cat("\nEstimasi Sigma(v)^2 dari Bean =",sigmav2estBean)##
## Estimasi Sigma(v)^2 dari Bean = 219.322sigmae2estBean <- BHFBean$est$fit$errorvar
cat("\nEstimasi Sigma(e)^2 dari Bean =",sigmae2estBean)##
## Estimasi Sigma(e)^2 dari Bean = 170.2856# Ambil untuk domain tertentu
mseBean1 <- BHFBean$mse[domains, ]
mseBean1## domain mse
## 1 1 145.03062
## 5 5 64.62653
## 10 10 30.01241residBean1 <- BHFBean$est$fit$residuals[domains]
cat("\nResidual dari Bean untuk Subset Counties (domain = 1, 5, 10) :\n")##
## Residual dari Bean untuk Subset Counties (domain = 1, 5, 10) :residBean1## [1] -5.849687 1.450064 -25.929810## ============================================================
## Pendugaan dari Soy Beans untuk Subset of Counties
## menggunakan ML Method - MANUAL
## ============================================================
## Persiapan Data
y_bean <- cornsoybean$SoyBeansHec
X_bean <- cbind(1, cornsoybean$CornPix, cornsoybean$SoyBeansPix)
dom <- cornsoybean$County
areas <- sort(unique(dom))
m <- length(areas)
ni <- as.vector(table(dom)); names(ni) <- areas
Xbar_pop <- Xmean[order(Xmean$CountyIndex), ]
ybar_ia_bean <- tapply(y_bean, dom, mean)
Xbar_ia_bean <- aggregate(
data.frame(CornPix = cornsoybean$CornPix,
SoyBeansPix = cornsoybean$SoyBeansPix),
by = list(dom = dom),
FUN = mean
)
## ============================================================
## LANGKAH 1: Estimasi Komponen Varians (sigma_v^2, sigma_e^2)
## Menggunakan ML (bukan REML) sesuai method="ML"
## ============================================================
data_model_bean <- data.frame(
y = y_bean,
Xs1 = cornsoybean$CornPix,
Xs2 = cornsoybean$SoyBeansPix,
dom = factor(dom)
)
fit_lme_bean <- lmer(y ~ Xs1 + Xs2 + (1 | dom),
data = data_model_bean,
REML = FALSE) # FALSE = ML method
summary(fit_lme_bean)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ Xs1 + Xs2 + (1 | dom)
## Data: data_model_bean
##
## AIC BIC logLik -2*log(L) df.resid
## 323.2 331.2 -156.6 313.2 32
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9871 -0.5487 -0.1716 0.4525 1.8864
##
## Random effects:
## Groups Name Variance Std.Dev.
## dom (Intercept) 219.3 14.81
## Residual 170.3 13.05
## Number of obs: 37, groups: dom, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -16.34566 25.46384 -0.642
## Xs1 0.02807 0.05348 0.525
## Xs2 0.49679 0.05537 8.972
##
## Correlation of Fixed Effects:
## (Intr) Xs1
## Xs1 -0.935
## Xs2 -0.878 0.718vc_bean <- as.data.frame(VarCorr(fit_lme_bean))
sigma_v2_bean <- vc_bean$vcov[vc_bean$grp == "dom"]
sigma_e2_bean <- vc_bean$vcov[vc_bean$grp == "Residual"]
cat("=== Estimasi Komponen Varians (ML) ===\n")## === Estimasi Komponen Varians (ML) ===cat("Estimasi Sigma(v)^2 dari Bean =", sigma_v2_bean, "\n")## Estimasi Sigma(v)^2 dari Bean = 219.322cat("Estimasi Sigma(e)^2 dari Bean =", sigma_e2_bean, "\n\n")## Estimasi Sigma(e)^2 dari Bean = 170.2856## ============================================================
## LANGKAH 2: Hitung gamma_hat_i untuk SEMUA area
## gamma_hat_i = sigma_v^2 / (sigma_v^2 + sigma_e^2 / n_i)
## ============================================================
a_i_bean <- ni
gamma_hat_bean <- sigma_v2_bean / (sigma_v2_bean + sigma_e2_bean / a_i_bean)
cat("=== Gamma_hat per Area (semua) ===\n")## === Gamma_hat per Area (semua) ===print(data.frame(
Area = areas,
ni = ni,
gamma_hat = round(gamma_hat_bean, 6)
), row.names = FALSE)## Area ni gamma_hat
## 1 1 0.562931
## 2 1 0.562931
## 3 1 0.562931
## 4 2 0.720353
## 5 3 0.794404
## 6 3 0.794404
## 7 3 0.794404
## 8 3 0.794404
## 9 4 0.837448
## 10 5 0.865588
## 11 5 0.865588
## 12 6 0.885424cat("\n")## ============================================================
## LANGKAH 3: Fixed Effects beta_hat (ML)
## ============================================================
beta_hat_bean <- as.vector(fixef(fit_lme_bean))
names(beta_hat_bean) <- c("Intercept", "CornPix", "SoyBeansPix")
cat("=== Beta_hat dari Bean (ML) ===\n")## === Beta_hat dari Bean (ML) ===print(beta_hat_bean)## Intercept CornPix SoyBeansPix
## -16.34565637 0.02806513 0.49678511cat("\n")## ============================================================
## LANGKAH 4: Hitung EBLUP untuk SEMUA area
## Persamaan (17):
## mu_hat_i = gamma_i * [ybar_ia + (Xbar_i - xbar_ia)' beta_hat]
## + (1 - gamma_i) * Xbar_i' * beta_hat
## ============================================================
eblup_bean_all <- numeric(m)
names(eblup_bean_all) <- areas
for (i in seq_along(areas)) {
area_id <- areas[i]
ybar_i <- ybar_ia_bean[as.character(area_id)]
row_ia <- Xbar_ia_bean[Xbar_ia_bean$dom == area_id, ]
xbar_ia_i <- c(row_ia$CornPix, row_ia$SoyBeansPix)
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
synth_i <- beta_hat_bean["Intercept"] + sum(Xbar_i * beta_hat_bean[-1])
survey_i <- ybar_i + sum((Xbar_i - xbar_ia_i) * beta_hat_bean[-1])
eblup_bean_all[i] <- gamma_hat_bean[i] * survey_i +
(1 - gamma_hat_bean[i]) * synth_i
}
## Ambil hanya untuk subset domains = c(1, 5, 10)
domains <- c(1, 5, 10)
idx_domains <- which(areas %in% domains)
eblup_bean_subset <- eblup_bean_all[idx_domains]
cat("=== EBLUP Manual Bean - Semua Area ===\n")## === EBLUP Manual Bean - Semua Area ===print(data.frame(
domain = areas,
eblup = round(eblup_bean_all, 4),
sampsize = ni
), row.names = FALSE)## domain eblup sampsize
## 1 78.6476 1
## 2 94.4001 1
## 3 87.3722 1
## 4 81.2040 2
## 5 66.2879 3
## 6 113.7497 3
## 7 97.7800 3
## 8 112.2923 3
## 9 109.7911 4
## 10 100.5984 5
## 11 118.9629 5
## 12 74.8851 6cat("\n=== EBLUP Manual Bean - Subset (domain = 1, 5, 10) ===\n")##
## === EBLUP Manual Bean - Subset (domain = 1, 5, 10) ===print(data.frame(
domain = areas[idx_domains],
eblup = round(eblup_bean_subset, 4),
sampsize = ni[idx_domains]
), row.names = FALSE)## domain eblup sampsize
## 1 78.6476 1
## 5 66.2879 3
## 10 100.5984 5## ============================================================
## LANGKAH 5: Random Effects v_hat_i (semua area)
## v_hat_i = gamma_i * (ybar_ia - xbar_ia' * beta_hat)
## ============================================================
v_hat_bean <- numeric(m)
names(v_hat_bean) <- areas
for (i in seq_along(areas)) {
area_id <- areas[i]
ybar_i <- ybar_ia_bean[as.character(area_id)]
row_ia <- Xbar_ia_bean[Xbar_ia_bean$dom == area_id, ]
xbar_ia_i <- c(row_ia$CornPix, row_ia$SoyBeansPix)
xbar_pred <- beta_hat_bean["Intercept"] + sum(xbar_ia_i * beta_hat_bean[-1])
v_hat_bean[i] <- gamma_hat_bean[i] * (ybar_i - xbar_pred)
}
cat("\n=== Random Effects Manual Bean (v_hat_i) ===\n")##
## === Random Effects Manual Bean (v_hat_i) ===print(round(v_hat_bean, 6))## 1 2 3 4 5 6 7
## -7.534196 4.622161 -6.389879 -20.012059 -19.722436 0.107377 13.848055
## 8 9 10 11 12
## 10.214657 -3.971755 9.432391 25.279699 -5.874015cat("\n=== Random Effects dari lme4 (pembanding) ===\n")##
## === Random Effects dari lme4 (pembanding) ===print(round(as.numeric(ranef(fit_lme_bean)$dom[, 1]), 6))## [1] -7.534196 4.622161 -6.389879 -20.012059 -19.722436 0.107377
## [7] 13.848055 10.214657 -3.971755 9.432391 25.279699 -5.874015## ============================================================
## LANGKAH 6: Residual Manual (semua observasi)
## resid_ij = y_ij - (X_ij * beta_hat + v_hat_i)
## ============================================================
v_hat_ij_bean <- rep(v_hat_bean, times = ni)
fitted_bean <- as.numeric(X_bean %*% beta_hat_bean) + v_hat_ij_bean
resid_bean_manual <- y_bean - fitted_bean
cat("\n=== Residual Manual Bean (All Counties) ===\n")##
## === Residual Manual Bean (All Counties) ===print(round(resid_bean_manual, 6))## 1 2 3 4 4 5 5
## -5.849687 3.588730 -4.961218 -16.987790 1.450064 1.377022 -2.922396
## 5 6 6 6 7 7 7
## -13.767482 7.417710 -25.929810 18.595470 24.615987 -2.239022 -11.625084
## 8 8 8 9 9 9 9
## -7.163347 21.309744 -6.215552 6.869112 5.904820 -7.160299 -8.697376
## 10 10 10 10 10 11 11
## 2.889711 -9.562734 22.439237 -4.819388 -3.623346 -7.049549 18.153970
## 11 11 11 12 12 12 12
## 0.880127 -4.817206 12.460276 -11.522354 -1.082776 2.050934 4.255165
## 12 12
## -10.759369 12.497707## Residual untuk subset domain
obs_subset <- which(dom %in% domains)
cat("\n=== Residual Manual Bean (Subset domain = 1, 5, 10) ===\n")##
## === Residual Manual Bean (Subset domain = 1, 5, 10) ===print(round(resid_bean_manual[obs_subset], 6))## 1 5 5 5 10 10 10
## -5.849687 1.377022 -2.922396 -13.767482 2.889711 -9.562734 22.439237
## 10 10
## -4.819388 -3.623346## ============================================================
## LANGKAH 8: Bootstrap MSE Manual (B = 200) - Semua Area
## ============================================================
set.seed(123)
B <- 200
mse_boot_bean <- numeric(m)
names(mse_boot_bean) <- areas
for (b in 1:B) {
# Langkah 2: Bangkitkan v_i* dan e_ij*
v_star <- rnorm(m, 0, sqrt(sigma_v2_bean))
e_star <- rnorm(sum(ni), 0, sqrt(sigma_e2_bean))
# Langkah 3: y_ij* = X_ij * beta_hat + v_i* + e_ij*
v_star_ij <- rep(v_star, times = ni)
y_star <- as.numeric(X_bean %*% beta_hat_bean) + v_star_ij + e_star
# Langkah 4: mu_i*(b) = Xbar_i' * beta_hat + v_i*
mu_star <- numeric(m)
for (i in seq_along(areas)) {
area_id <- areas[i]
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
mu_star[i] <- beta_hat_bean["Intercept"] +
sum(Xbar_i * beta_hat_bean[-1]) +
v_star[i]
}
# Langkah 5: Estimasi ulang parameter dari data bootstrap (ML)
data_boot <- data.frame(
y_star = y_star,
Xs1 = cornsoybean$CornPix,
Xs2 = cornsoybean$SoyBeansPix,
dom = factor(dom)
)
fit_boot <- tryCatch(
lmer(y_star ~ Xs1 + Xs2 + (1 | dom),
data = data_boot, REML = FALSE), # ML method
error = function(e) NULL
)
if (is.null(fit_boot)) next
vc_boot <- as.data.frame(VarCorr(fit_boot))
sigma_v2_b <- max(vc_boot$vcov[vc_boot$grp == "dom"], 0)
sigma_e2_b <- vc_boot$vcov[vc_boot$grp == "Residual"]
beta_hat_b <- as.vector(fixef(fit_boot))
gamma_b <- sigma_v2_b / (sigma_v2_b + sigma_e2_b / a_i_bean)
ybar_boot <- tapply(y_star, dom, mean)
Xbar_ia_b <- aggregate(
data.frame(CornPix = cornsoybean$CornPix,
SoyBeansPix = cornsoybean$SoyBeansPix),
by = list(dom = dom),
FUN = mean
)
# Langkah 6: EBLUP bootstrap
mu_hat_star <- numeric(m)
for (i in seq_along(areas)) {
area_id <- areas[i]
ybar_b_i <- ybar_boot[as.character(area_id)]
row_ia_b <- Xbar_ia_b[Xbar_ia_b$dom == area_id, ]
xbar_ia_b <- c(row_ia_b$CornPix, row_ia_b$SoyBeansPix)
row_pop <- Xbar_pop[Xbar_pop$CountyIndex == area_id, ]
Xbar_i <- c(row_pop$MeanCornPixPerSeg, row_pop$MeanSoyBeansPixPerSeg)
synth_b <- beta_hat_b[1] + sum(Xbar_i * beta_hat_b[-1])
survey_b <- ybar_b_i + sum((Xbar_i - xbar_ia_b) * beta_hat_b[-1])
mu_hat_star[i] <- gamma_b[i] * survey_b +
(1 - gamma_b[i]) * synth_b
}
# Langkah 8: Akumulasi MSE
mse_boot_bean <- mse_boot_bean + (mu_hat_star - mu_star)^2
}## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')mse_boot_bean <- mse_boot_bean / B
RMSE_bean <- sqrt(mse_boot_bean)
## ============================================================
## LANGKAH 9: Ringkasan Hasil
## ============================================================
# Semua area
cat("\n=== MSE Bootstrap Manual Bean - Semua Area ===\n")##
## === MSE Bootstrap Manual Bean - Semua Area ===mse_all <- data.frame(
domain = areas,
eblup = round(eblup_bean_all, 4),
mse = round(mse_boot_bean, 4),
rmse = round(RMSE_bean, 4),
cv_pct = round(100 * RMSE_bean / abs(eblup_bean_all), 2)
)
print(mse_all, row.names = FALSE)## domain eblup mse rmse cv_pct
## 1 78.6476 119.4393 10.9288 13.90
## 2 94.4001 121.9107 11.0413 11.70
## 3 87.3722 122.2808 11.0581 12.66
## 4 81.2040 71.1901 8.4374 10.39
## 5 66.2879 51.4310 7.1715 10.82
## 6 113.7497 47.3275 6.8795 6.05
## 7 97.7800 56.1159 7.4911 7.66
## 8 112.2923 60.6782 7.7896 6.94
## 9 109.7911 34.5603 5.8788 5.35
## 10 100.5984 33.3578 5.7756 5.74
## 11 118.9629 35.5731 5.9643 5.01
## 12 74.8851 34.7854 5.8979 7.88# Subset domain = c(10, 1, 5)
mse_subset <- mse_all[mse_all$domain %in% domains, ]
cat("\n=== MSE Bootstrap Manual Bean - Subset (domain = 1, 5, 10) ===\n")##
## === MSE Bootstrap Manual Bean - Subset (domain = 1, 5, 10) ===print(mse_subset, row.names = FALSE)## domain eblup mse rmse cv_pct
## 1 78.6476 119.4393 10.9288 13.90
## 5 66.2879 51.4310 7.1715 10.82
## 10 100.5984 33.3578 5.7756 5.74# Residual dan komponen varians
residBean_manual <- resid_bean_manual
sigmav2estBean_manual <- sigma_v2_bean
sigmae2estBean_manual <- sigma_e2_bean
cat("\nResidual dari Bean (All Counties) :\n")##
## Residual dari Bean (All Counties) :print(round(residBean_manual, 6))## 1 2 3 4 4 5 5
## -5.849687 3.588730 -4.961218 -16.987790 1.450064 1.377022 -2.922396
## 5 6 6 6 7 7 7
## -13.767482 7.417710 -25.929810 18.595470 24.615987 -2.239022 -11.625084
## 8 8 8 9 9 9 9
## -7.163347 21.309744 -6.215552 6.869112 5.904820 -7.160299 -8.697376
## 10 10 10 10 10 11 11
## 2.889711 -9.562734 22.439237 -4.819388 -3.623346 -7.049549 18.153970
## 11 11 11 12 12 12 12
## 0.880127 -4.817206 12.460276 -11.522354 -1.082776 2.050934 4.255165
## 12 12
## -10.759369 12.497707cat("\nResidual dari Bean untuk Subset Counties (domain = 1, 5, 10) :\n")##
## Residual dari Bean untuk Subset Counties (domain = 1, 5, 10) :print(round(residBean_manual[obs_subset], 6))## 1 5 5 5 10 10 10
## -5.849687 1.377022 -2.922396 -13.767482 2.889711 -9.562734 22.439237
## 10 10
## -4.819388 -3.623346cat("\nEstimasi Sigma(v)^2 dari Bean =", sigmav2estBean_manual)##
## Estimasi Sigma(v)^2 dari Bean = 219.322cat("\nEstimasi Sigma(e)^2 dari Bean =", sigmae2estBean_manual, "\n")##
## Estimasi Sigma(e)^2 dari Bean = 170.2856## ============================================================
## LANGKAH 10: Verifikasi dengan pbmseBHF()
## ============================================================
## Subset Counties (domain = 10, 1, 5)\n")
compare_bean_subset <- data.frame(
Domain = areas[idx_domains],
EBLUP_Manual = round(eblup_bean_subset, 4),
EBLUP_Fungsi = round(resultBean$eblup$eblup, 4),
MSE_Manual = round(mse_boot_bean[idx_domains], 4),
MSE_Fungsi = round(mseBean1$mse, 4)
)
kable(compare_bean_subset,
caption = "Perbandingan EBLUP dan MSE Manual vs pbmseBHF() dari Estimasi Corn - Subset (Domain = 1, 5, 10)",
col.names = c("Domain", "EBLUP Manual", "EBLUP Fungsi",
"MSE Manual", "MSE Fungsi"),
align = c("c","r","r","r","r"),
booktabs = TRUE) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "bordered", "condensed"),
full_width = FALSE,
position = "center"
) %>%
row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")| Domain | EBLUP Manual | EBLUP Fungsi | MSE Manual | MSE Fungsi | |
|---|---|---|---|---|---|
| 1 | 1 | 78.6476 | 78.6369 | 119.4393 | 145.0306 |
| 5 | 5 | 66.2879 | 66.2608 | 51.4310 | 64.6265 |
| 10 | 10 | 100.5984 | 100.6112 | 33.3578 | 30.0124 |