EBLUP-SAE
Empirical Best Linear Unbiased Prediction- Small Area Estimation
Package
Package yang perlu digunakan dalam analisis ini yaitu : sae.
## Warning: package 'sae' was built under R version 4.0.5## Loading required package: MASS## Warning: package 'MASS' was built under R version 4.0.4## Loading required package: lme4## Loading required package: MatrixData
Data yang digunakan berasal dari data survey jagung dan kacang kedelai, dan data satelit di 12 Kabupaten di Iowa. Data ini tersedia di R pada package sae. Akan dilakukan pendugaan area kecil terhadap data ini.
## 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## 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.05(Xmean <- data.frame(cornsoybeanmeans[, c("CountyIndex","MeanCornPixPerSeg","MeanSoyBeansPixPerSeg")]))## 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.05## 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 556#Call function that returns EBLUP estimates of the means of corn crop area and also parametric bootstrap MSE estimates with B=200 bootstrap replicates for all the counties, without removing any data point. We consider the default fitting method, REML.
set.seed(123)
BHF <- pbmseBHF(CornHec ~ CornPix + SoyBeansPix, dom = County,
meanxpop = Xmean, popnsize = Popn, B = 200, data = cornsoybean) #y = CorHec , x1=CornPix, x2=SoyBeansPix##
## 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 ?isSingular## b = 8
## b = 9
## b = 10## boundary (singular) fit: see ?isSingular## b = 11
## b = 12
## b = 13
## b = 14
## b = 15
## b = 16
## b = 17## boundary (singular) fit: see ?isSingular## b = 18
## b = 19## boundary (singular) fit: see ?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 ?isSingular## b = 30
## b = 31
## b = 32## boundary (singular) fit: see ?isSingular## b = 33## boundary (singular) fit: see ?isSingular## b = 34
## b = 35
## b = 36
## b = 37
## b = 38
## b = 39
## b = 40## boundary (singular) fit: see ?isSingular## b = 41
## b = 42
## b = 43
## b = 44
## b = 45
## b = 46
## b = 47## boundary (singular) fit: see ?isSingular## b = 48
## b = 49
## b = 50
## b = 51
## b = 52
## b = 53
## b = 54## boundary (singular) fit: see ?isSingular## b = 55
## b = 56
## b = 57
## b = 58
## b = 59## boundary (singular) fit: see ?isSingular## b = 60
## b = 61
## b = 62
## b = 63
## b = 64
## b = 65## boundary (singular) fit: see ?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 ?isSingular## b = 76
## b = 77## boundary (singular) fit: see ?isSingular## b = 78
## b = 79## boundary (singular) fit: see ?isSingular## b = 80
## b = 81## boundary (singular) fit: see ?isSingular## b = 82
## b = 83
## b = 84
## b = 85## boundary (singular) fit: see ?isSingular## b = 86
## b = 87
## b = 88
## b = 89
## b = 90
## b = 91
## b = 92## boundary (singular) fit: see ?isSingular## b = 93## boundary (singular) fit: see ?isSingular## b = 94
## b = 95## boundary (singular) fit: see ?isSingular## b = 96## boundary (singular) fit: see ?isSingular## b = 97## boundary (singular) fit: see ?isSingular## b = 98
## b = 99
## b = 100
## b = 101
## b = 102
## b = 103
## b = 104
## b = 105## boundary (singular) fit: see ?isSingular## b = 106## boundary (singular) fit: see ?isSingular## b = 107## boundary (singular) fit: see ?isSingular## b = 108## boundary (singular) fit: see ?isSingular## b = 109## boundary (singular) fit: see ?isSingular## b = 110
## b = 111
## b = 112
## b = 113## boundary (singular) fit: see ?isSingular## b = 114
## b = 115## boundary (singular) fit: see ?isSingular## b = 116
## b = 117
## b = 118
## b = 119
## b = 120## boundary (singular) fit: see ?isSingular## b = 121
## b = 122
## b = 123
## b = 124
## b = 125## boundary (singular) fit: see ?isSingular## b = 126## boundary (singular) fit: see ?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 ?isSingular## b = 138
## b = 139## boundary (singular) fit: see ?isSingular## b = 140
## b = 141
## b = 142
## b = 143
## b = 144## boundary (singular) fit: see ?isSingular## b = 145## boundary (singular) fit: see ?isSingular## b = 146## boundary (singular) fit: see ?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 ?isSingular## b = 157
## b = 158
## b = 159
## b = 160## boundary (singular) fit: see ?isSingular## b = 161## boundary (singular) fit: see ?isSingular## b = 162
## b = 163
## b = 164
## b = 165
## b = 166
## b = 167
## b = 168## boundary (singular) fit: see ?isSingular## b = 169
## b = 170
## b = 171
## b = 172
## b = 173
## b = 174
## b = 175## boundary (singular) fit: see ?isSingular## b = 176## boundary (singular) fit: see ?isSingular## b = 177
## b = 178## boundary (singular) fit: see ?isSingular## b = 179
## b = 180
## b = 181
## b = 182
## b = 183
## b = 184
## b = 185## boundary (singular) fit: see ?isSingular## b = 186## boundary (singular) fit: see ?isSingular## b = 187
## b = 188
## b = 189
## b = 190## boundary (singular) fit: see ?isSingular## b = 191## boundary (singular) fit: see ?isSingular## b = 192
## b = 193
## b = 194
## b = 195
## b = 196## boundary (singular) fit: see ?isSingular## b = 197
## b = 198
## b = 199
## b = 200#Menghitung sqrtmse
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 9.219720
## 2 Hamilton 1 123.5274 8.640516
## 3 Worth 1 113.0343 8.653565
## 4 Humboldt 2 114.9901 8.020315
## 5 Franklin 3 137.2660 8.319550
## 6 Pocahontas 3 108.9807 7.131785
## 7 Winnebago 3 116.4839 6.852536
## 8 Wright 3 122.7711 7.029138
## 9 Webster 4 111.5648 6.787478
## 10 Hancock 5 124.1565 6.245703
## 11 Kossuth 5 112.4626 6.157168
## 12 Hardin 6 131.2515 5.955844Dapat 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.
Referensi
Battese GE, Harter RM, Fuller WA. 1988. An error-components model for prediction of county crop areas using survey and satellite data. J Am Stat Assoc. 83(401):28–36. doi:10.1080/01621459.1988.10478561.
J.N.K R, Molina I. 2015. Small Area Estimation. Second Edi. Hoboken, New Jersey: John Wiley & Sons, Inc.
IPB University, nadirabelinda@apps.ipb.ac.id↩︎