An Application to HB Twofold Normal Model On sampel dataset

FIRST STEPS: Load package and load the data

library(saeHB.twofold)
data("dataTwofold")

STEPS 2: Fitting HB Model

model=NormalTF(y~x1+x2,vardir="vardir",area = "codearea",weight = "w",iter.mcmc = 100000,thin=50,burn.in = 1000,data=dataTwofold)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 90
#>    Unobserved stochastic nodes: 125
#>    Total graph size: 953
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 90
#>    Unobserved stochastic nodes: 125
#>    Total graph size: 953
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 90
#>    Unobserved stochastic nodes: 125
#>    Total graph size: 953
#> 
#> Initializing model
#> Warning in par(opar): argument 1 does not name a graphical parameter

STEP 3 Extract mean estimation

Subarea Estimation

model$Est_sub
#>                Mean        SD        2.5%       25%       50%       75%
#> theta[1]  13.594925 3.0140631  7.73686315 11.551874 13.559450 15.631783
#> theta[2]  11.064887 1.9050634  7.35792715  9.756531 11.006045 12.295356
#> theta[3]  14.143834 2.5191976  9.30302236 12.508716 14.165294 15.795477
#> theta[4]  11.588684 1.2991965  8.98706214 10.748213 11.599420 12.471958
#> theta[5]  14.765808 0.7725453 13.27936210 14.251547 14.775313 15.270731
#> theta[6]  10.685405 1.2332284  8.21738673  9.873609 10.693293 11.532595
#> theta[7]  14.252884 0.9086451 12.48089503 13.594362 14.244441 14.891069
#> theta[8]  13.465864 1.8599347  9.83134450 12.217462 13.470692 14.710690
#> theta[9]  16.165675 0.8917237 14.47721455 15.565701 16.170923 16.729571
#> theta[10] 14.842212 0.8951357 13.10008741 14.215020 14.836125 15.470635
#> theta[11] 14.280482 0.7139409 12.89204160 13.794373 14.286636 14.750315
#> theta[12] 15.842296 1.9815892 11.93600404 14.496921 15.841736 17.133793
#> theta[13]  9.276981 0.5757624  8.15296051  8.904147  9.259437  9.656904
#> theta[14]  9.914109 0.9606422  8.12190041  9.242204  9.901735 10.580929
#> theta[15] 15.520237 0.7126655 14.19215616 15.009948 15.520819 16.007885
#> theta[16] 12.650745 1.5113664  9.65435560 11.640640 12.633881 13.661343
#> theta[17]  2.245791 0.6930789  0.84966215  1.791284  2.238840  2.697205
#> theta[18]  7.182724 1.1458057  4.94961612  6.436011  7.177194  7.973712
#> theta[19]  2.143186 0.7328289  0.72955215  1.646969  2.134490  2.636649
#> theta[20]  6.073198 1.8321595  2.44113689  4.891935  6.078508  7.323520
#> theta[21]  4.440928 0.7722255  2.91570134  3.922980  4.415796  4.966334
#> theta[22] 14.729809 1.2234791 12.36484878 13.903162 14.693355 15.553985
#> theta[23] 13.962236 0.4771211 13.02582636 13.640802 13.954050 14.297394
#> theta[24] 17.668737 1.1113312 15.51322365 16.931086 17.677727 18.420879
#> theta[25] 14.481489 1.3659974 11.85551208 13.557308 14.483160 15.393528
#> theta[26] 10.330787 1.9440861  6.37564484  9.065498 10.335846 11.637550
#> theta[27]  9.467229 2.7887581  4.10369802  7.676445  9.459300 11.252248
#> theta[28] 14.983357 0.9008644 13.27143359 14.374554 15.002214 15.574075
#> theta[29] 13.133026 1.6868601  9.83401593 11.986960 13.106057 14.270921
#> theta[30] 10.195117 1.6979623  6.88156540  9.059268 10.217965 11.359241
#> theta[31] 14.094618 2.2291707  9.84621770 12.556145 14.038233 15.696938
#> theta[32] 14.564141 0.9429044 12.70928309 13.924778 14.577002 15.221126
#> theta[33] 10.099279 3.0888263  3.88189109  8.073014 10.108282 12.189431
#> theta[34] 13.753739 1.3884788 10.97645095 12.845085 13.729050 14.661599
#> theta[35] 14.620853 1.9943240 10.60437335 13.295121 14.617990 15.950824
#> theta[36] 13.440315 2.0042289  9.65681505 12.084639 13.432032 14.811116
#> theta[37] 13.985349 0.8373531 12.36677450 13.417616 13.979874 14.530412
#> theta[38] 14.254148 0.9105075 12.44237487 13.628614 14.283507 14.856019
#> theta[39] 12.520193 1.4182394  9.69464389 11.557072 12.574364 13.471578
#> theta[40]  9.228416 2.3506490  4.55059067  7.657218  9.213859 10.802326
#> theta[41] 18.413693 0.8173686 16.79566899 17.860746 18.426607 18.963054
#> theta[42] 12.494478 1.0095185 10.55091561 11.827602 12.476862 13.163182
#> theta[43] 16.453004 3.2902988 10.13177402 14.278573 16.445045 18.680615
#> theta[44] 16.985550 1.4711446 14.21896543 15.994025 16.941956 18.005003
#> theta[45] 18.526444 1.3515378 15.84903643 17.605213 18.515324 19.468375
#> theta[46] 13.629115 1.9591559  9.86315388 12.270123 13.644670 14.931376
#> theta[47] 12.369762 2.6944616  6.97998343 10.533556 12.327385 14.187493
#> theta[48] 11.380950 0.6720117 10.04950916 10.913343 11.373483 11.850325
#> theta[49]  6.720785 1.4005713  3.92130412  5.803455  6.712906  7.622371
#> theta[50] 15.991670 1.5450380 12.99566662 14.962120 15.981333 16.994847
#> theta[51]  7.791652 1.4050370  5.09960193  6.818275  7.822791  8.703789
#> theta[52] 26.602922 0.8525628 24.97821757 26.037090 26.597917 27.193020
#> theta[53] 18.524322 1.6635328 15.17485031 17.448693 18.591206 19.649551
#> theta[54] 17.944526 1.2464942 15.40767246 17.121816 17.961286 18.789026
#> theta[55] 11.969594 0.7042348 10.59637989 11.498212 11.954251 12.453746
#> theta[56] 11.485084 0.5251993 10.46155087 11.134862 11.479730 11.838593
#> theta[57]  8.999329 0.6634485  7.68275325  8.569431  8.992873  9.433523
#> theta[58] 13.821943 1.4175039 11.04636139 12.866163 13.832513 14.781652
#> theta[59] 16.483252 0.9239204 14.58815788 15.875678 16.480842 17.105509
#> theta[60] 15.179707 0.8204032 13.57797236 14.639971 15.179128 15.721794
#> theta[61]  5.622955 1.9781039  1.71321841  4.325171  5.671322  6.898275
#> theta[62]  5.924883 0.8100535  4.30883760  5.381724  5.901289  6.499494
#> theta[63]  5.340135 0.9265505  3.54997602  4.711355  5.302034  5.993711
#> theta[64] 12.197445 1.5573168  9.20703797 11.142945 12.197317 13.263971
#> theta[65] 12.366848 2.4432263  7.41606499 10.745833 12.427157 14.028886
#> theta[66]  9.886264 1.2524121  7.44682304  9.056517  9.928099 10.746564
#> theta[67] 16.803160 1.0586993 14.67456579 16.079396 16.796378 17.498837
#> theta[68] 11.156523 1.3003548  8.45446038 10.330106 11.165389 12.039924
#> theta[69]  8.166392 0.8559034  6.47222791  7.605533  8.143536  8.772598
#> theta[70] 11.759499 0.8508413 10.21513395 11.189965 11.732145 12.305514
#> theta[71] 13.713015 0.7547404 12.20199419 13.208503 13.694420 14.221959
#> theta[72] 10.734262 0.7118026  9.35656141 10.260615 10.737245 11.197560
#> theta[73]  3.829813 1.4882190  0.91836296  2.809947  3.865795  4.864632
#> theta[74] 10.427344 1.5269173  7.46777712  9.400282 10.444643 11.427381
#> theta[75]  5.968607 1.3117129  3.37312282  5.068742  5.963021  6.892261
#> theta[76]  9.453906 0.8020740  7.87767017  8.910760  9.454904 10.014344
#> theta[77] 13.915757 0.7685914 12.39952428 13.420836 13.921933 14.420964
#> theta[78]  9.570049 1.9158717  5.87124577  8.307479  9.577342 10.859396
#> theta[79] 10.452481 1.1567403  8.24793396  9.675661 10.460652 11.260498
#> theta[80]  5.913930 0.7241750  4.52284219  5.428540  5.898262  6.403777
#> theta[81]  6.512096 0.7464968  5.01840520  6.020531  6.517719  7.012687
#> theta[82] 14.513518 2.2519577 10.16802963 12.935743 14.522494 15.964739
#> theta[83] 10.680166 1.8405165  7.19023951  9.447423 10.634482 11.909999
#> theta[84] 10.063014 3.0624393  3.89703447  8.104000 10.052307 12.082477
#> theta[85] 13.417926 1.6403937 10.02832163 12.334080 13.445474 14.550774
#> theta[86] 18.563076 1.4528671 15.73697403 17.551715 18.567203 19.513985
#> theta[87] 11.474234 0.8068654  9.86852291 10.946006 11.447724 12.002438
#> theta[88]  7.941102 1.8936991  4.33357536  6.666516  7.959537  9.212012
#> theta[89]  7.700753 1.4410456  4.77410324  6.744063  7.710097  8.661038
#> theta[90]  4.482209 2.3010720  0.01801621  2.875045  4.478305  5.950563
#>               97.5%
#> theta[1]  19.559197
#> theta[2]  14.856442
#> theta[3]  19.318089
#> theta[4]  14.024365
#> theta[5]  16.269898
#> theta[6]  13.099843
#> theta[7]  15.992347
#> theta[8]  16.990561
#> theta[9]  17.864972
#> theta[10] 16.593871
#> theta[11] 15.703884
#> theta[12] 19.768989
#> theta[13] 10.439705
#> theta[14] 11.772837
#> theta[15] 16.863365
#> theta[16] 15.631477
#> theta[17]  3.566384
#> theta[18]  9.493284
#> theta[19]  3.571691
#> theta[20]  9.584369
#> theta[21]  5.917274
#> theta[22] 17.137336
#> theta[23] 14.858395
#> theta[24] 19.868061
#> theta[25] 17.265784
#> theta[26] 14.293813
#> theta[27] 14.917854
#> theta[28] 16.737735
#> theta[29] 16.442006
#> theta[30] 13.538537
#> theta[31] 18.324460
#> theta[32] 16.368409
#> theta[33] 16.204646
#> theta[34] 16.581959
#> theta[35] 18.600225
#> theta[36] 17.217238
#> theta[37] 15.652501
#> theta[38] 15.977705
#> theta[39] 15.218313
#> theta[40] 13.799354
#> theta[41] 20.015090
#> theta[42] 14.465053
#> theta[43] 22.843410
#> theta[44] 19.881296
#> theta[45] 21.120094
#> theta[46] 17.787116
#> theta[47] 17.689560
#> theta[48] 12.661778
#> theta[49]  9.516662
#> theta[50] 18.955563
#> theta[51] 10.584787
#> theta[52] 28.294958
#> theta[53] 21.716692
#> theta[54] 20.396017
#> theta[55] 13.356967
#> theta[56] 12.484065
#> theta[57] 10.260430
#> theta[58] 16.500204
#> theta[59] 18.248132
#> theta[60] 16.826243
#> theta[61]  9.551614
#> theta[62]  7.474333
#> theta[63]  7.151661
#> theta[64] 15.159575
#> theta[65] 17.146081
#> theta[66] 12.253958
#> theta[67] 18.910577
#> theta[68] 13.693866
#> theta[69]  9.837791
#> theta[70] 13.476756
#> theta[71] 15.251085
#> theta[72] 12.123975
#> theta[73]  6.693362
#> theta[74] 13.601872
#> theta[75]  8.540874
#> theta[76] 10.954572
#> theta[77] 15.414114
#> theta[78] 13.248390
#> theta[79] 12.679222
#> theta[80]  7.362659
#> theta[81]  7.954839
#> theta[82] 18.988695
#> theta[83] 14.292940
#> theta[84] 16.314887
#> theta[85] 16.563549
#> theta[86] 21.347980
#> theta[87] 13.121585
#> theta[88] 11.632028
#> theta[89] 10.511383
#> theta[90]  8.982347

Area Estimatio

model$Est_area
#>         Mean        SD      2.5%       25%       50%       75%     97.5%
#> 1  13.099276 1.8850409  9.462958 11.793300 13.129487 14.346726 16.956476
#> 2  12.457517 0.6515855 11.150841 12.022501 12.455950 12.892502 13.730511
#> 3  14.452526 0.7633072 12.930624 13.932443 14.446680 14.941051 15.982173
#> 4  14.875075 0.6784071 13.549130 14.426863 14.853017 15.341079 16.249243
#> 5  11.284189 0.4456603 10.422869 10.969697 11.284890 11.596549 12.124735
#> 6   7.827009 0.7487982  6.307533  7.345965  7.838038  8.330070  9.274511
#> 7   4.516273 0.8307645  2.878086  3.945435  4.522484  5.082251  6.107967
#> 8  15.890341 0.6416436 14.607750 15.448032 15.893309 16.319424 17.152263
#> 9  11.333789 1.4047479  8.572614 10.426547 11.354692 12.238701 14.051188
#> 10 12.322630 0.9721767 10.443832 11.665316 12.317999 12.983821 14.216820
#> 11 12.643947 1.6214277  9.365978 11.526015 12.652210 13.732703 15.834397
#> 12 13.949870 1.1327413 11.801542 13.182352 13.939223 14.720917 16.105355
#> 13 13.471964 0.6995338 12.098402 12.987252 13.491312 13.950516 14.809763
#> 14 13.265514 0.9597172 11.407371 12.630559 13.257948 13.921583 15.101482
#> 15 17.270554 1.5783011 14.243248 16.192500 17.269420 18.352425 20.473195
#> 16 12.584368 1.3549792  9.992088 11.672637 12.544796 13.477735 15.270399
#> 17 10.511800 0.8997008  8.762217  9.903676 10.492546 11.120618 12.323349
#> 18 20.810601 0.7786697 19.201619 20.313700 20.815785 21.342293 22.292728
#> 19 10.412981 0.3998483  9.629327 10.144399 10.409999 10.678137 11.207021
#> 20 15.258885 0.6004720 14.075205 14.867168 15.267637 15.657189 16.452805
#> 21  5.600738 0.8284109  3.907860  5.067858  5.604660  6.156470  7.242019
#> 22 11.502115 1.1002907  9.292076 10.759208 11.529555 12.253517 13.635545
#> 23 12.093536 0.6476149 10.821965 11.649497 12.088973 12.538122 13.321288
#> 24 12.117668 0.4790957 11.190923 11.797504 12.097557 12.440302 13.082102
#> 25  6.487406 0.8996504  4.727583  5.843034  6.505914  7.104391  8.200421
#> 26 11.354133 0.7043550  9.934102 10.902321 11.345158 11.835203 12.722291
#> 27  7.482920 0.5332027  6.442763  7.119192  7.499681  7.840200  8.530109
#> 28 12.149125 1.5567836  9.129199 11.105523 12.156175 13.170711 15.353275
#> 29 14.739689 0.8328503 13.074206 14.172567 14.742696 15.297527 16.355481
#> 30  6.484850 1.2703475  4.065159  5.596096  6.441222  7.330158  9.003678

Coefficient Estimation

model$coefficient
#>               Mean         SD       2.5%        25%        50%       75%
#> beta[0] -0.2702159 0.66139593 -1.5978795 -0.7241061 -0.2670289 0.1898495
#> beta[1]  0.2348534 0.61285485 -0.9305834 -0.1714123  0.2283585 0.6613784
#> beta[2]  1.2170637 0.06871058  1.0831428  1.1714959  1.2159849 1.2615295
#>             97.5%
#> beta[0] 0.9653253
#> beta[1] 1.4304783
#> beta[2] 1.3564975

Random effect variance estimation

model$refVar
#>   var_area  var_sub
#> 1 8.092291 8.793771

STEP 3 : Extract CV and MSE Subarea

CV $CV(\theta \hat)=\frac{SD(\theta \hat)}{(\theta \hat)} \times 100$
MSE MSE = V(θ)̂

CV=(model$Est_sub$SD)/(model$Est_sub$Mean)*100
MSE=model$Est_sub$SD^2
summary(cbind(CV,MSE))
#>        CV              MSE         
#>  Min.   : 3.205   Min.   : 0.2276  
#>  1st Qu.: 6.513   1st Qu.: 0.7069  
#>  Median :11.502   Median : 1.6894  
#>  Mean   :13.692   Mean   : 2.3730  
#>  3rd Qu.:17.379   3rd Qu.: 3.3798  
#>  Max.   :51.338   Max.   :10.8261

STEP 4 Extract CV and MSE of Area

CV2=(model$Est_area$SD)/(model$Est_area$Mean)*100
MSE2=model$Est_area$SD^2
summary(cbind(CV2,MSE2))
#>       CV2              MSE2       
#>  Min.   : 3.742   Min.   :0.1599  
#>  1st Qu.: 5.202   1st Qu.:0.4335  
#>  Median : 7.562   Median :0.6882  
#>  Mean   : 8.599   Mean   :1.0084  
#>  3rd Qu.:11.988   3rd Qu.:1.2650  
#>  Max.   :19.589   Max.   :3.5534

STEP 5 : You can also compare the CV between subarea direct estimator and HB Twofold estimator

dirCV=sqrt(dataTwofold$vardir)/(dataTwofold$y)*100
summary(cbind(dirCV,CV))
#>      dirCV               CV        
#>  Min.   :  3.187   Min.   : 3.205  
#>  1st Qu.:  6.642   1st Qu.: 6.513  
#>  Median : 12.764   Median :11.502  
#>  Mean   : 25.141   Mean   :13.692  
#>  3rd Qu.: 21.041   3rd Qu.:17.379  
#>  Max.   :680.378   Max.   :51.338

Look that! ,CV on our model is less than CV on direct estimator.

boxplot(cbind(dirCV,CV),ylim=c(0,50))

model$refVar
#>   var_area  var_sub
#> 1 8.092291 8.793771