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.700836 2.9841553  7.84397932 11.710340 13.715867 15.752670
#> theta[2]  11.002184 1.9177401  7.21284169  9.753640 11.032156 12.277713
#> theta[3]  14.293842 2.5135478  9.39279508 12.581547 14.338298 15.952506
#> theta[4]  11.603413 1.3746126  8.94809116 10.677383 11.598682 12.523812
#> theta[5]  14.765726 0.7757108 13.26068276 14.264171 14.773048 15.267743
#> theta[6]  10.647656 1.2441374  8.19825323  9.805757 10.657267 11.456757
#> theta[7]  14.302558 0.9181777 12.50821974 13.657989 14.309011 14.967690
#> theta[8]  13.463353 1.7860722  9.81117055 12.295585 13.488433 14.704016
#> theta[9]  16.195471 0.9162524 14.37998332 15.606028 16.195125 16.791448
#> theta[10] 14.868310 0.8778040 13.15198983 14.286343 14.864006 15.464710
#> theta[11] 14.279204 0.7006086 12.90568991 13.824211 14.276854 14.761533
#> theta[12] 15.797201 1.9585923 12.19814520 14.441074 15.753379 17.125426
#> theta[13]  9.264908 0.5819150  8.08774108  8.867728  9.268531  9.665873
#> theta[14]  9.884987 0.9729736  8.06028220  9.242607  9.861275 10.523704
#> theta[15] 15.556360 0.7073490 14.15213699 15.084257 15.566087 16.028095
#> theta[16] 12.592962 1.5091298  9.74180976 11.573265 12.603995 13.618402
#> theta[17]  2.240940 0.6934159  0.86549126  1.772983  2.258542  2.699123
#> theta[18]  7.182423 1.1477621  4.92911569  6.391046  7.198753  7.950405
#> theta[19]  2.140899 0.7315419  0.72577852  1.660936  2.146606  2.623694
#> theta[20]  6.092369 1.8735669  2.48732960  4.871362  6.069239  7.343218
#> theta[21]  4.382645 0.7847955  2.87506109  3.856958  4.391940  4.922319
#> theta[22] 14.763542 1.1659052 12.46168373 14.010233 14.774807 15.523308
#> theta[23] 13.979654 0.4891300 13.04876754 13.653168 13.988402 14.302781
#> theta[24] 17.621004 1.1157127 15.46274578 16.853739 17.614291 18.393386
#> theta[25] 14.467239 1.3136979 11.87221136 13.582232 14.467912 15.355930
#> theta[26] 10.350990 1.9962971  6.50560259  9.033527 10.298760 11.738736
#> theta[27]  9.600212 2.7409839  4.30801645  7.745921  9.564851 11.451550
#> theta[28] 14.992532 0.9124515 13.21804478 14.355291 15.012871 15.611041
#> theta[29] 13.228416 1.6398461 10.00886112 12.128460 13.275009 14.362726
#> theta[30] 10.089158 1.7100308  6.70987803  8.938736 10.070834 11.199461
#> theta[31] 14.105885 2.2385186 10.01221700 12.515254 14.084827 15.702305
#> theta[32] 14.598148 0.9202476 12.78277735 13.961146 14.612962 15.238612
#> theta[33]  9.861159 3.2494574  3.48290948  7.811813  9.918704 11.985244
#> theta[34] 13.802577 1.4163881 10.91915536 12.887296 13.793869 14.761708
#> theta[35] 14.622883 2.0498503 10.68981085 13.224925 14.627310 15.997240
#> theta[36] 13.460671 1.9468090  9.72255160 12.115582 13.453638 14.800826
#> theta[37] 13.970652 0.8394202 12.34230781 13.416130 13.980232 14.544799
#> theta[38] 14.271862 0.9245124 12.39262358 13.659733 14.249844 14.859097
#> theta[39] 12.474313 1.4042908  9.71209645 11.522691 12.435372 13.443709
#> theta[40]  9.172924 2.3438900  4.64532172  7.635095  9.161632 10.782794
#> theta[41] 18.451202 0.8225654 16.88824932 17.906831 18.442864 18.973087
#> theta[42] 12.517591 0.9958116 10.57104796 11.854765 12.521124 13.189190
#> theta[43] 16.441635 3.3713821  9.98643467 14.216332 16.474472 18.763984
#> theta[44] 16.995192 1.4685020 14.09231644 16.006534 16.943888 18.013495
#> theta[45] 18.522659 1.4152881 15.71690807 17.598067 18.521091 19.472287
#> theta[46] 13.694495 1.9680612  9.88535342 12.388310 13.741360 15.026926
#> theta[47] 12.459937 2.6884725  6.91258690 10.771005 12.471745 14.249369
#> theta[48] 11.381127 0.6365273 10.10601276 10.954261 11.384943 11.817928
#> theta[49]  6.576709 1.3810977  3.77871944  5.686474  6.587215  7.532944
#> theta[50] 15.914073 1.5315639 12.94313296 14.870424 15.905990 16.945538
#> theta[51]  7.746773 1.3947616  5.04641247  6.777415  7.732308  8.722022
#> theta[52] 26.608486 0.8481536 24.95193100 26.042762 26.597180 27.189558
#> theta[53] 18.568659 1.6743363 15.33341380 17.381816 18.593495 19.737058
#> theta[54] 17.919667 1.2735884 15.49658614 16.997518 17.931536 18.777829
#> theta[55] 11.987320 0.7457585 10.58052842 11.464200 11.974263 12.476034
#> theta[56] 11.480232 0.5366446 10.43217795 11.123793 11.462312 11.839456
#> theta[57]  8.974027 0.6585201  7.64928852  8.539617  8.969026  9.421485
#> theta[58] 13.748271 1.4014663 11.02307384 12.801677 13.748821 14.706119
#> theta[59] 16.499646 0.9392728 14.70273982 15.864520 16.492476 17.138120
#> theta[60] 15.170610 0.8331571 13.54068204 14.589362 15.154640 15.720865
#> theta[61]  5.562194 2.0192344  1.66756416  4.242036  5.544260  6.838406
#> theta[62]  5.921018 0.8152810  4.33439861  5.360361  5.929732  6.466860
#> theta[63]  5.325055 0.9431702  3.46060466  4.715807  5.324508  5.956944
#> theta[64] 12.320606 1.5036587  9.29910576 11.343397 12.351762 13.328397
#> theta[65] 12.346010 2.4354655  7.51980840 10.751590 12.333615 13.997261
#> theta[66]  9.869536 1.2491205  7.46501914  9.002784  9.866286 10.739826
#> theta[67] 16.784198 1.0484244 14.71772304 16.067718 16.782611 17.482436
#> theta[68] 11.121854 1.2869124  8.60766738 10.278632 11.117700 11.995835
#> theta[69]  8.155602 0.8728289  6.47051964  7.569918  8.142711  8.765172
#> theta[70] 11.757205 0.8402124 10.12742932 11.197752 11.762027 12.317114
#> theta[71] 13.727313 0.7261703 12.26337797 13.253340 13.738274 14.189806
#> theta[72] 10.764833 0.7209311  9.37593513 10.251159 10.767605 11.254039
#> theta[73]  3.847679 1.4875354  1.08681135  2.813547  3.828363  4.834862
#> theta[74] 10.377174 1.4879010  7.43336517  9.381432 10.340932 11.381611
#> theta[75]  5.924570 1.3060877  3.34781244  5.066262  5.884776  6.789324
#> theta[76]  9.439930 0.7974144  7.84695331  8.894490  9.423539  9.996444
#> theta[77] 13.939418 0.7524845 12.45122022 13.453085 13.958963 14.449581
#> theta[78]  9.609245 1.9763412  5.63167388  8.296957  9.664640 10.951302
#> theta[79] 10.448898 1.1452686  8.17984749  9.692578 10.434544 11.212619
#> theta[80]  5.879242 0.7363014  4.41888910  5.383451  5.884098  6.382580
#> theta[81]  6.485499 0.7315356  5.08918573  5.975724  6.475850  6.979795
#> theta[82] 14.517446 2.2752274 10.02092391 13.056631 14.520009 15.982429
#> theta[83] 10.732904 1.8414765  7.06382094  9.476700 10.709408 11.970493
#> theta[84] 10.023458 3.0588079  4.11494364  8.005588  9.959581 12.103175
#> theta[85] 13.359665 1.5444237 10.36065047 12.330022 13.377616 14.400086
#> theta[86] 18.556646 1.4559635 15.79455279 17.548890 18.559535 19.536169
#> theta[87] 11.505079 0.7779917  9.95716038 10.998972 11.502114 12.019273
#> theta[88]  7.998143 1.8828039  4.34509426  6.756058  7.987460  9.217200
#> theta[89]  7.736669 1.4034570  4.97786882  6.828128  7.738249  8.692154
#> theta[90]  4.549741 2.3030862  0.09235355  2.980458  4.591112  6.141857
#>               97.5%
#> theta[1]  19.553260
#> theta[2]  14.691170
#> theta[3]  19.261200
#> theta[4]  14.276343
#> theta[5]  16.319966
#> theta[6]  13.098977
#> theta[7]  16.118226
#> theta[8]  16.911840
#> theta[9]  18.047082
#> theta[10] 16.617872
#> theta[11] 15.639207
#> theta[12] 19.743498
#> theta[13] 10.341023
#> theta[14] 11.742808
#> theta[15] 16.931937
#> theta[16] 15.616784
#> theta[17]  3.568188
#> theta[18]  9.388288
#> theta[19]  3.581800
#> theta[20]  9.788045
#> theta[21]  5.849933
#> theta[22] 17.059409
#> theta[23] 14.942191
#> theta[24] 19.834414
#> theta[25] 17.066379
#> theta[26] 14.192439
#> theta[27] 15.013204
#> theta[28] 16.780795
#> theta[29] 16.388945
#> theta[30] 13.410506
#> theta[31] 18.559874
#> theta[32] 16.351543
#> theta[33] 16.266627
#> theta[34] 16.638363
#> theta[35] 18.626196
#> theta[36] 17.222022
#> theta[37] 15.592394
#> theta[38] 16.122887
#> theta[39] 15.332106
#> theta[40] 13.730763
#> theta[41] 20.176509
#> theta[42] 14.419135
#> theta[43] 22.797035
#> theta[44] 19.933279
#> theta[45] 21.252088
#> theta[46] 17.349979
#> theta[47] 17.584733
#> theta[48] 12.607343
#> theta[49]  9.159185
#> theta[50] 18.953635
#> theta[51] 10.463506
#> theta[52] 28.243657
#> theta[53] 21.832715
#> theta[54] 20.394616
#> theta[55] 13.475695
#> theta[56] 12.533681
#> theta[57] 10.260468
#> theta[58] 16.385521
#> theta[59] 18.321176
#> theta[60] 16.822831
#> theta[61]  9.688778
#> theta[62]  7.500279
#> theta[63]  7.156373
#> theta[64] 15.204075
#> theta[65] 16.976851
#> theta[66] 12.288227
#> theta[67] 18.805221
#> theta[68] 13.627819
#> theta[69]  9.877345
#> theta[70] 13.379226
#> theta[71] 15.145427
#> theta[72] 12.175670
#> theta[73]  6.763252
#> theta[74] 13.297443
#> theta[75]  8.565210
#> theta[76] 10.977161
#> theta[77] 15.372398
#> theta[78] 13.468731
#> theta[79] 12.676088
#> theta[80]  7.296614
#> theta[81]  7.977884
#> theta[82] 19.134967
#> theta[83] 14.330841
#> theta[84] 15.933701
#> theta[85] 16.371435
#> theta[86] 21.435820
#> theta[87] 13.026198
#> theta[88] 11.784134
#> theta[89] 10.573122
#> theta[90]  8.980101

Area Estimatio

model$Est_area
#>         Mean        SD      2.5%       25%       50%       75%     97.5%
#> 1  13.174194 1.8757970  9.548572 11.906090 13.156644 14.464995 16.987577
#> 2  12.447131 0.6608199 11.183743 11.991859 12.440342 12.905846 13.711261
#> 3  14.480509 0.7697081 12.912427 13.965633 14.492558 15.001342 15.942461
#> 4  14.872766 0.6683304 13.576497 14.422683 14.859455 15.316530 16.240171
#> 5  11.280297 0.4504848 10.379203 10.979672 11.286902 11.577315 12.152841
#> 6   7.803630 0.7432633  6.342944  7.315838  7.813511  8.287237  9.329164
#> 7   4.504830 0.8690714  2.810087  3.922831  4.489988  5.100613  6.160734
#> 8  15.881963 0.6406304 14.642258 15.463884 15.886584 16.324258 17.113437
#> 9  11.383365 1.4085046  8.601455 10.406262 11.389193 12.340254 14.120491
#> 10 12.303727 0.9620901 10.418536 11.638255 12.298204 12.947526 14.230625
#> 11 12.563215 1.6785112  9.370077 11.416515 12.541644 13.669699 15.979627
#> 12 13.974230 1.1285062 11.848715 13.185825 13.945363 14.767715 16.190132
#> 13 13.453780 0.6983896 12.100951 12.960998 13.465686 13.930321 14.833792
#> 14 13.266242 0.9487562 11.426317 12.607216 13.273275 13.910274 15.106527
#> 15 17.267055 1.6176639 14.120503 16.191423 17.269713 18.402143 20.340883
#> 16 12.642532 1.3658930  9.939696 11.713622 12.674476 13.577839 15.291727
#> 17 10.428405 0.8844264  8.693297  9.826081 10.446626 11.026856 12.132657
#> 18 20.814389 0.7952688 19.246545 20.269274 20.820103 21.351741 22.380848
#> 19 10.404270 0.4047848  9.602860 10.139229 10.407243 10.684260 11.186211
#> 20 15.242160 0.6028092 14.071268 14.843170 15.242522 15.645692 16.381621
#> 21  5.572848 0.8592665  3.911412  5.020063  5.566820  6.135076  7.312850
#> 22 11.539506 1.0687694  9.472401 10.827886 11.552108 12.267741 13.631946
#> 23 12.072095 0.6353944 10.789805 11.652326 12.079158 12.512064 13.311277
#> 24 12.128729 0.4773417 11.219069 11.814361 12.123462 12.449714 13.050944
#> 25  6.463888 0.8689617  4.757401  5.892280  6.473107  7.064715  8.173768
#> 26 11.371467 0.7127147  9.964253 10.906070 11.384509 11.860617 12.777731
#> 27  7.460122 0.5221459  6.476078  7.113395  7.457077  7.824783  8.479261
#> 28 12.158085 1.5714062  9.016493 11.108547 12.136820 13.179602 15.311834
#> 29 14.729116 0.7954282 13.194907 14.189566 14.723236 15.262712 16.300551
#> 30  6.539990 1.2589493  4.018163  5.695882  6.546468  7.412784  8.933607

Coefficient Estimation

model$coefficient
#>               Mean         SD       2.5%        25%        50%        75%
#> beta[0] -0.3784285 0.64703757 -1.6256174 -0.8166177 -0.3862564 0.06250906
#> beta[1]  0.2118237 0.60754539 -0.9587698 -0.1940662  0.2202596 0.61795121
#> beta[2]  1.2290059 0.06518217  1.1036651  1.1846098  1.2304438 1.27184298
#>             97.5%
#> beta[0] 0.8938012
#> beta[1] 1.3864935
#> beta[2] 1.3566824

Random effect variance estimation

model$refVar
#>   var_area  var_sub
#> 1 8.159361 8.813393

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.188   Min.   : 0.2392  
#>  1st Qu.: 6.533   1st Qu.: 0.7050  
#>  Median :11.566   Median : 1.6810  
#>  Mean   :13.713   Mean   : 2.3821  
#>  3rd Qu.:17.680   3rd Qu.: 3.3408  
#>  Max.   :50.620   Max.   :11.3662

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.821   Min.   :0.1639  
#>  1st Qu.: 5.209   1st Qu.:0.4392  
#>  Median : 7.486   Median :0.6855  
#>  Mean   : 8.622   Mean   :1.0154  
#>  3rd Qu.:11.981   3rd Qu.:1.2407  
#>  Max.   :19.292   Max.   :3.5186

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.188  
#>  1st Qu.:  6.642   1st Qu.: 6.533  
#>  Median : 12.764   Median :11.566  
#>  Mean   : 25.141   Mean   :13.713  
#>  3rd Qu.: 21.041   3rd Qu.:17.680  
#>  Max.   :680.378   Max.   :50.620
  • 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.159361 8.813393