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.741504 3.0856670  7.8010353 11.659458 13.628460 15.787649
#> theta[2]  11.092691 1.9180724  7.1987404  9.820287 11.076013 12.362363
#> theta[3]  14.390827 2.5807087  9.4263030 12.665893 14.378202 16.108445
#> theta[4]  11.616381 1.3407315  8.9751695 10.686310 11.649095 12.520909
#> theta[5]  14.738931 0.7763137 13.2796075 14.198006 14.736296 15.290461
#> theta[6]  10.612000 1.2665430  8.0427400  9.763193 10.593674 11.487157
#> theta[7]  14.292080 0.8951624 12.5591964 13.676577 14.291456 14.902012
#> theta[8]  13.400855 1.8274872  9.8102010 12.174250 13.411870 14.579576
#> theta[9]  16.206453 0.9155210 14.3750868 15.618583 16.207820 16.794910
#> theta[10] 14.823813 0.8775305 13.1284323 14.235533 14.831045 15.382206
#> theta[11] 14.311051 0.6918043 12.9617292 13.841279 14.308388 14.789031
#> theta[12] 15.809711 1.9796102 11.9443534 14.464636 15.784353 17.236204
#> theta[13]  9.298998 0.5743344  8.1466978  8.919245  9.317439  9.682927
#> theta[14]  9.950196 0.9681224  7.9964212  9.312942  9.968393 10.584745
#> theta[15] 15.554675 0.6902498 14.2569966 15.088441 15.541135 15.985033
#> theta[16] 12.611936 1.5818383  9.5868703 11.526827 12.556409 13.671945
#> theta[17]  2.196914 0.7177179  0.7356928  1.714375  2.197089  2.691536
#> theta[18]  7.199177 1.1583580  4.8340849  6.448731  7.229045  7.959109
#> theta[19]  2.110965 0.7401162  0.6279908  1.632817  2.106957  2.604813
#> theta[20]  6.127710 1.8474213  2.6039379  4.841596  6.125601  7.390807
#> theta[21]  4.432889 0.7471209  3.0421213  3.928373  4.423203  4.919864
#> theta[22] 14.731722 1.1448974 12.4479135 13.949484 14.719196 15.536869
#> theta[23] 13.968729 0.4663789 13.0723168 13.646933 13.983192 14.281234
#> theta[24] 17.581238 1.0904061 15.4822527 16.859843 17.588733 18.294606
#> theta[25] 14.421590 1.3653716 11.7664966 13.495209 14.431883 15.349361
#> theta[26] 10.293635 1.9819774  6.4557988  8.971032 10.257156 11.638256
#> theta[27]  9.587082 2.8040436  4.0455986  7.751201  9.629002 11.486004
#> theta[28] 14.961907 0.8858421 13.2460682 14.368029 14.957921 15.567672
#> theta[29] 13.149861 1.6264879 10.0691727 12.022490 13.137077 14.256480
#> theta[30] 10.058354 1.6669207  6.9290321  8.887520 10.040820 11.205201
#> theta[31] 14.181591 2.2317677  9.8054831 12.628145 14.146691 15.722344
#> theta[32] 14.546117 0.9322131 12.7579346 13.897288 14.570186 15.184316
#> theta[33]  9.903317 3.1597719  3.4715728  7.861396  9.858348 11.904885
#> theta[34] 13.829367 1.3910253 11.1708362 12.905833 13.822698 14.789266
#> theta[35] 14.590212 2.0240139 10.4946809 13.255110 14.597462 15.926169
#> theta[36] 13.482963 1.9611633  9.6321472 12.150202 13.490641 14.806365
#> theta[37] 14.000376 0.8334945 12.3410026 13.453484 14.005062 14.558568
#> theta[38] 14.269839 0.8853946 12.4413704 13.667886 14.239397 14.882568
#> theta[39] 12.507644 1.4089969  9.6965009 11.600640 12.506736 13.412383
#> theta[40]  9.181506 2.3807535  4.4869077  7.585705  9.197934 10.760113
#> theta[41] 18.445731 0.8272858 16.8229265 17.885651 18.436509 19.008748
#> theta[42] 12.467545 1.0135991 10.5182864 11.807416 12.439559 13.159509
#> theta[43] 16.522281 3.3168856 10.1014102 14.301455 16.596986 18.746549
#> theta[44] 16.971693 1.4635211 14.1430622 15.962318 17.000347 17.967128
#> theta[45] 18.562263 1.4438475 15.7648396 17.598033 18.548465 19.508032
#> theta[46] 13.663864 1.9287392  9.8171269 12.381917 13.713335 14.965290
#> theta[47] 12.409930 2.7075959  7.0391760 10.626483 12.407064 14.212290
#> theta[48] 11.394637 0.6561391 10.1079715 10.938157 11.414267 11.846956
#> theta[49]  6.659018 1.4036437  3.9463892  5.752088  6.668480  7.604093
#> theta[50] 15.988940 1.5442126 12.9361201 14.949975 16.000873 17.056198
#> theta[51]  7.826956 1.4054536  5.1440886  6.851027  7.839901  8.790994
#> theta[52] 26.572399 0.8431722 24.8730613 26.017482 26.566171 27.149396
#> theta[53] 18.510158 1.6919785 15.2329200 17.377127 18.491017 19.640350
#> theta[54] 17.963218 1.2641465 15.5514710 17.121944 17.944516 18.830777
#> theta[55] 12.009766 0.7173436 10.6078349 11.539437 12.011572 12.467499
#> theta[56] 11.480137 0.5294561 10.4190093 11.120376 11.480707 11.847100
#> theta[57]  8.987959 0.6613076  7.6625364  8.560964  8.993705  9.430382
#> theta[58] 13.797735 1.3943081 11.0673117 12.850670 13.833256 14.751014
#> theta[59] 16.501714 0.9341327 14.6860785 15.892548 16.468423 17.118164
#> theta[60] 15.168645 0.8156395 13.6175660 14.618071 15.140668 15.719714
#> theta[61]  5.563405 2.0129117  1.6492211  4.201345  5.593858  6.883907
#> theta[62]  5.951582 0.8066634  4.4002582  5.402388  5.936851  6.494711
#> theta[63]  5.339710 0.9288898  3.5870637  4.700898  5.339310  5.983513
#> theta[64] 12.265493 1.5649799  9.2136633 11.177799 12.303381 13.311575
#> theta[65] 12.364880 2.4287056  7.5618266 10.720560 12.336434 13.989128
#> theta[66]  9.891618 1.2686746  7.3352581  9.028408  9.920222 10.733844
#> theta[67] 16.781051 1.0500601 14.7874640 16.067549 16.765845 17.484125
#> theta[68] 11.122807 1.3174622  8.6125753 10.236462 11.092071 12.045382
#> theta[69]  8.135012 0.8556041  6.4453440  7.564122  8.146733  8.691202
#> theta[70] 11.756617 0.8551771 10.1280694 11.170586 11.770651 12.339956
#> theta[71] 13.703290 0.7530676 12.2464640 13.169604 13.726747 14.211451
#> theta[72] 10.758303 0.7224027  9.2890558 10.274157 10.756603 11.222581
#> theta[73]  3.903756 1.5148677  1.0237087  2.923638  3.865361  4.928876
#> theta[74] 10.415687 1.5178687  7.3164983  9.434598 10.430529 11.469050
#> theta[75]  5.967400 1.2927477  3.4238148  5.093783  5.982464  6.816727
#> theta[76]  9.439124 0.7673781  7.9216811  8.915954  9.436281  9.952396
#> theta[77] 13.927077 0.7474939 12.4704197 13.411831 13.921043 14.428489
#> theta[78]  9.757071 1.9246293  5.8423982  8.476305  9.774556 10.997492
#> theta[79] 10.404234 1.1833064  8.2148893  9.580308 10.372055 11.241866
#> theta[80]  5.875624 0.7340487  4.4989059  5.378228  5.850410  6.351004
#> theta[81]  6.482501 0.7085809  5.0672843  5.994592  6.492989  6.955405
#> theta[82] 14.410758 2.2247430  9.8853769 12.972243 14.408887 15.887631
#> theta[83] 10.635753 1.8612518  7.0391150  9.348748 10.666094 11.914285
#> theta[84] 10.001386 3.0433782  4.1723398  7.957520  9.929638 11.996383
#> theta[85] 13.302238 1.6148062 10.2492096 12.167669 13.299020 14.392525
#> theta[86] 18.591262 1.4382957 15.7332995 17.578191 18.610121 19.531044
#> theta[87] 11.449633 0.8033221  9.8616914 10.908095 11.468799 11.992462
#> theta[88]  7.964436 1.8868242  4.0918727  6.757296  7.967302  9.295286
#> theta[89]  7.681279 1.4029638  5.0110224  6.707305  7.672076  8.635890
#> theta[90]  4.523632 2.3913413 -0.1444548  2.844499  4.515239  6.158400
#>               97.5%
#> theta[1]  20.015861
#> theta[2]  14.897712
#> theta[3]  19.347867
#> theta[4]  14.153514
#> theta[5]  16.216615
#> theta[6]  13.065719
#> theta[7]  16.046213
#> theta[8]  17.081288
#> theta[9]  18.036070
#> theta[10] 16.596235
#> theta[11] 15.644146
#> theta[12] 19.723922
#> theta[13] 10.388192
#> theta[14] 11.897317
#> theta[15] 16.984344
#> theta[16] 15.701274
#> theta[17]  3.584280
#> theta[18]  9.474577
#> theta[19]  3.626581
#> theta[20]  9.786767
#> theta[21]  5.921857
#> theta[22] 16.944627
#> theta[23] 14.887307
#> theta[24] 19.753551
#> theta[25] 17.076075
#> theta[26] 14.169452
#> theta[27] 15.005441
#> theta[28] 16.706840
#> theta[29] 16.371411
#> theta[30] 13.377479
#> theta[31] 18.652941
#> theta[32] 16.326672
#> theta[33] 16.323539
#> theta[34] 16.479561
#> theta[35] 18.620701
#> theta[36] 17.308179
#> theta[37] 15.596444
#> theta[38] 16.036171
#> theta[39] 15.321049
#> theta[40] 13.888526
#> theta[41] 20.061178
#> theta[42] 14.475080
#> theta[43] 22.992079
#> theta[44] 19.780602
#> theta[45] 21.497620
#> theta[46] 17.524941
#> theta[47] 17.829480
#> theta[48] 12.658452
#> theta[49]  9.432133
#> theta[50] 18.964178
#> theta[51] 10.476779
#> theta[52] 28.215289
#> theta[53] 21.756685
#> theta[54] 20.422197
#> theta[55] 13.447903
#> theta[56] 12.456897
#> theta[57] 10.265050
#> theta[58] 16.414720
#> theta[59] 18.359245
#> theta[60] 16.815917
#> theta[61]  9.608259
#> theta[62]  7.523168
#> theta[63]  7.201795
#> theta[64] 15.184530
#> theta[65] 17.099401
#> theta[66] 12.301053
#> theta[67] 18.851151
#> theta[68] 13.691623
#> theta[69]  9.849734
#> theta[70] 13.364485
#> theta[71] 15.189320
#> theta[72] 12.161693
#> theta[73]  6.969618
#> theta[74] 13.344254
#> theta[75]  8.504967
#> theta[76] 10.924388
#> theta[77] 15.354051
#> theta[78] 13.452206
#> theta[79] 12.661736
#> theta[80]  7.379711
#> theta[81]  7.857724
#> theta[82] 18.850880
#> theta[83] 14.260728
#> theta[84] 15.982645
#> theta[85] 16.500090
#> theta[86] 21.480337
#> theta[87] 13.014733
#> theta[88] 11.387332
#> theta[89] 10.454630
#> theta[90]  9.165464

Area Estimatio

model$Est_area
#>         Mean        SD      2.5%       25%       50%       75%     97.5%
#> 1  13.242848 1.9396623  9.546461 11.911779 13.213304 14.528353 17.091161
#> 2  12.426957 0.6904870 11.062640 11.978513 12.415188 12.894155 13.760647
#> 3  14.458236 0.7673141 13.020126 13.933154 14.459347 14.975624 15.972032
#> 4  14.872075 0.6673939 13.553206 14.427664 14.862498 15.328910 16.153920
#> 5  11.314462 0.4371514 10.472102 11.028518 11.329109 11.584075 12.143375
#> 6   7.803667 0.7540292  6.343227  7.297956  7.796631  8.320330  9.265336
#> 7   4.527867 0.8523857  2.872556  3.944267  4.525766  5.122739  6.202536
#> 8  15.851982 0.6278229 14.622670 15.417892 15.847975 16.268400 17.080973
#> 9  11.345494 1.3753111  8.675782 10.426557 11.340472 12.268617 14.040803
#> 10 12.259857 0.9441913 10.464042 11.598917 12.253974 12.900808 14.173578
#> 11 12.597674 1.6704112  9.297845 11.514585 12.584774 13.663293 15.928169
#> 12 13.979619 1.1223177 11.829874 13.212717 14.015906 14.718718 16.180917
#> 13 13.476104 0.6805209 12.102727 13.012681 13.474426 13.938956 14.791554
#> 14 13.250575 0.9599723 11.380285 12.606016 13.257960 13.895077 15.131747
#> 15 17.307562 1.6252821 14.131704 16.189273 17.318940 18.427987 20.388657
#> 16 12.615811 1.3321746 10.000563 11.691563 12.648289 13.504182 15.212228
#> 17 10.507208 0.8888326  8.781417  9.876870 10.498204 11.134173 12.179840
#> 18 20.804826 0.7917722 19.250232 20.283235 20.819011 21.317708 22.359175
#> 19 10.416968 0.4009617  9.620886 10.148577 10.429693 10.677361 11.152057
#> 20 15.253864 0.5904001 14.116197 14.866325 15.236685 15.654661 16.414436
#> 21  5.587228 0.8415768  3.926810  5.027346  5.593024  6.156714  7.231207
#> 22 11.530172 1.0865787  9.460675 10.788108 11.531992 12.283105 13.722067
#> 23 12.064544 0.6483701 10.782049 11.622620 12.056756 12.494264 13.347357
#> 24 12.119319 0.4748259 11.173822 11.808419 12.116283 12.435712 13.058219
#> 25  6.510241 0.8809165  4.820643  5.923279  6.493985  7.102820  8.233856
#> 26 11.409352 0.6993754 10.057983 10.953762 11.411093 11.876456 12.837572
#> 27  7.444032 0.5347580  6.427998  7.070266  7.441987  7.825312  8.474654
#> 28 12.075691 1.5681859  8.881538 11.020570 12.062012 13.135378 15.123180
#> 29 14.706958 0.7991713 13.220845 14.152136 14.697842 15.226685 16.323904
#> 30  6.503061 1.3063727  3.870801  5.642313  6.528773  7.405501  8.950579

Coefficient Estimation

model$coefficient
#>               Mean         SD       2.5%        25%        50%        75%
#> beta[0] -0.4072923 0.62307565 -1.6015664 -0.8293815 -0.3999044 0.01528482
#> beta[1]  0.2502772 0.58513628 -0.9274328 -0.1229403  0.2368556 0.64960189
#> beta[2]  1.2276135 0.06511588  1.1022560  1.1844025  1.2244511 1.27212477
#>             97.5%
#> beta[0] 0.8124954
#> beta[1] 1.3920988
#> beta[2] 1.3600645

Random effect variance estimation

model$refVar
#>   var_area  var_sub
#> 1 8.090813 8.795565

STEP 3 : Extract CV and MSE Subarea

  • CV \(CV(\theta \hat)=\frac{SD(\theta \hat)}{(\theta \hat)} \times 100\)
  • MSE \(MSE= V(\theta\hat)\)
CV=(model$Est_sub$SD)/(model$Est_sub$Mean)*100
MSE=model$Est_sub$SD^2
summary(cbind(CV,MSE))
#>        CV              MSE         
#>  Min.   : 3.173   Min.   : 0.2175  
#>  1st Qu.: 6.485   1st Qu.: 0.6988  
#>  Median :11.693   Median : 1.7035  
#>  Mean   :13.750   Mean   : 2.3935  
#>  3rd Qu.:17.474   3rd Qu.: 3.3947  
#>  Max.   :52.863   Max.   :11.0017

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.806   Min.   :0.1608  
#>  1st Qu.: 5.114   1st Qu.:0.4498  
#>  Median : 7.473   Median :0.6735  
#>  Mean   : 8.626   Mean   :1.0201  
#>  3rd Qu.:11.731   3rd Qu.:1.2399  
#>  Max.   :20.089   Max.   :3.7623

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.173  
#>  1st Qu.:  6.642   1st Qu.: 6.485  
#>  Median : 12.763   Median :11.693  
#>  Mean   : 25.141   Mean   :13.750  
#>  3rd Qu.: 21.041   3rd Qu.:17.474  
#>  Max.   :680.378   Max.   :52.863
  • 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.090813 8.795565