saeHB-Twofold-Beta

STEP 1: Load package and load the data

library(saeHB.TF.beta)
data("dataBeta")

STEPS 2: Fitting HB Model

model <- betaTF(y~X1+X2,area="codearea",weight="w",iter.mcmc = 10000, burn.in = 3000, iter.update = 5, thin = 10, data=dataBeta)
#> 
#> SAMPLING FOR MODEL 'saeHB_TF_beta' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000133 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.33 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 3001 / 10000 [ 30%]  (Sampling)
#> Chain 1: Iteration: 4000 / 10000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 5000 / 10000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 4.239 seconds (Warm-up)
#> Chain 1:                13.986 seconds (Sampling)
#> Chain 1:                18.225 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'saeHB_TF_beta' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.67 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 3001 / 10000 [ 30%]  (Sampling)
#> Chain 1: Iteration: 4000 / 10000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 5000 / 10000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.67 seconds (Warm-up)
#> Chain 1:                12.25 seconds (Sampling)
#> Chain 1:                15.92 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'saeHB_TF_beta' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.64 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 3001 / 10000 [ 30%]  (Sampling)
#> Chain 1: Iteration: 4000 / 10000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 5000 / 10000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.339 seconds (Warm-up)
#> Chain 1:                9.313 seconds (Sampling)
#> Chain 1:                12.652 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'saeHB_TF_beta' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.0002 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 3001 / 10000 [ 30%]  (Sampling)
#> Chain 1: Iteration: 4000 / 10000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 5000 / 10000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.067 seconds (Warm-up)
#> Chain 1:                8.547 seconds (Sampling)
#> Chain 1:                11.614 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'saeHB_TF_beta' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.5e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.65 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 3001 / 10000 [ 30%]  (Sampling)
#> Chain 1: Iteration: 4000 / 10000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 5000 / 10000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.004 seconds (Warm-up)
#> Chain 1:                6.211 seconds (Sampling)
#> Chain 1:                9.215 seconds (Total)
#> Chain 1:

STEP 3 Extract mean estimation

Subarea Estimation

model$Est_sub
#>             Mean         SD       2.5%       25%       50%       75%     97.5%
#> mu[1]  0.9053967 0.08195363 0.67851291 0.8818738 0.9299119 0.9603876 0.9873144
#> mu[2]  0.8926761 0.08760498 0.67347240 0.8612844 0.9154443 0.9519162 0.9888949
#> mu[3]  0.8868311 0.08887440 0.64658461 0.8512343 0.9104661 0.9473585 0.9866026
#> mu[4]  0.8868978 0.08701021 0.64662824 0.8527519 0.9117412 0.9504478 0.9830285
#> mu[5]  0.6271048 0.18100766 0.25176080 0.5057484 0.6404679 0.7655775 0.9301059
#> mu[6]  0.8216952 0.12028181 0.52548477 0.7697127 0.8461044 0.9063645 0.9688645
#> mu[7]  0.8948383 0.07950368 0.68129007 0.8603116 0.9158148 0.9522523 0.9866266
#> mu[8]  0.8624380 0.10544918 0.58550770 0.8165252 0.8922498 0.9370547 0.9830706
#> mu[9]  0.8799658 0.09982128 0.60810802 0.8445298 0.9072289 0.9476080 0.9844358
#> mu[10] 0.9123285 0.07908197 0.70933065 0.8852371 0.9348038 0.9624647 0.9908243
#> mu[11] 0.8720506 0.09865071 0.62345252 0.8278944 0.8967816 0.9399233 0.9846169
#> mu[12] 0.9053396 0.07518627 0.70027774 0.8769720 0.9276158 0.9589379 0.9876333
#> mu[13] 0.6971537 0.16602980 0.34038058 0.5873610 0.7245334 0.8371638 0.9402407
#> mu[14] 0.8188202 0.13481296 0.44575324 0.7671898 0.8571600 0.9104273 0.9692113
#> mu[15] 0.4465486 0.21414551 0.10179592 0.2822292 0.4195753 0.6075303 0.8796106
#> mu[16] 0.8327054 0.11914243 0.52509691 0.7753740 0.8629299 0.9182501 0.9780020
#> mu[17] 0.9153811 0.07127910 0.74311794 0.8903543 0.9354156 0.9639012 0.9890464
#> mu[18] 0.7870216 0.14488772 0.43461564 0.7075130 0.8188011 0.8969234 0.9718825
#> mu[19] 0.3788355 0.21580093 0.07319142 0.2032868 0.3426698 0.5343229 0.8426713
#> mu[20] 0.8151413 0.12778213 0.50495844 0.7459182 0.8474112 0.9077456 0.9701964
#> mu[21] 0.5964160 0.18134710 0.20502766 0.4747901 0.6073351 0.7375049 0.8969731
#> mu[22] 0.8024685 0.12713643 0.49016011 0.7406594 0.8290650 0.8967925 0.9723601
#> mu[23] 0.7374890 0.16106601 0.34383560 0.6384648 0.7597792 0.8628978 0.9563236
#> mu[24] 0.6259265 0.17771714 0.24697717 0.4966795 0.6394002 0.7650359 0.9147318
#> mu[25] 0.4808298 0.20884347 0.12687534 0.3067006 0.4657572 0.6481988 0.8781021
#> mu[26] 0.4268665 0.19737120 0.12287057 0.2783191 0.4060867 0.5622995 0.8280090
#> mu[27] 0.6189552 0.18567270 0.21853089 0.4954805 0.6423776 0.7596318 0.9142918
#> mu[28] 0.9296569 0.05702640 0.77986863 0.9095508 0.9459108 0.9689667 0.9900320
#> mu[29] 0.8437936 0.11985869 0.53266660 0.7926927 0.8766417 0.9279530 0.9840234
#> mu[30] 0.8589059 0.10546324 0.57483154 0.8114975 0.8926066 0.9341297 0.9792554

Area Estimation

model$Est_area
#>         Mean         SD      2.5%       25%       50%       75%     97.5%
#> 1  1.2177115 0.08536626 1.0181561 1.1820636 1.2372036 1.2765637 1.3239864
#> 2  0.9708202 0.10644154 0.7296944 0.9083033 0.9816517 1.0510719 1.1291948
#> 3  1.2879826 0.10589976 1.0218098 1.2365061 1.3076122 1.3671666 1.4216273
#> 4  0.9760982 0.07245221 0.7960258 0.9429232 0.9915813 1.0288044 1.0688764
#> 5  0.9770270 0.16822147 0.6173117 0.8589678 0.9907493 1.0999273 1.2592528
#> 6  1.4803831 0.13373837 1.1445366 1.4061231 1.5037906 1.5816653 1.6666454
#> 7  0.8182944 0.18003381 0.4717322 0.6963216 0.8192384 0.9475808 1.1436643
#> 8  0.8426272 0.14515393 0.5350700 0.7483065 0.8512957 0.9512469 1.0822125
#> 9  0.6995611 0.18666905 0.3323952 0.5709749 0.6960454 0.8297694 1.0774089
#> 10 0.7984735 0.06038606 0.6443066 0.7691049 0.8108283 0.8424015 0.8804124

Coefficient

model$coefficient
#>           Mean        SD       2.5%        25%       50%       75%     97.5%
#> b[0] 0.8392249 0.2572229  0.3474709 0.66877717 0.8322050 1.0062950 1.3329864
#> b[1] 0.3031667 0.3551812 -0.4624135 0.07607504 0.3010123 0.5428577 0.9941086
#> b[2] 1.2698733 0.4122647  0.4746837 0.99754388 1.2794708 1.5544243 2.0308225

Random Effect

model$area_randeff
#>             f_mean
#> f[1]   0.255401039
#> f[2]   0.008793188
#> f[3]   0.539419782
#> f[4]   0.817424298
#> f[5]  -0.401040333
#> f[6]   0.528203138
#> f[7]  -0.765338417
#> f[8]  -0.365873029
#> f[9]  -1.014203212
#> f[10]  0.386653135
model$sub_randeff
#>            u_mean
#> u[1]   0.18075235
#> u[2]  -0.05199439
#> u[3]   0.24953733
#> u[4]   0.21384398
#> u[5]  -0.54239478
#> u[6]   0.27412731
#> u[7]  -0.07659570
#> u[8]   0.19446461
#> u[9]   0.33980887
#> u[10]  0.59489277
#> u[11]  0.06774146
#> u[12]  0.23690097
#> u[13] -0.19702884
#> u[14]  0.67010373
#> u[15] -0.90804435
#> u[16]  0.04379546
#> u[17]  0.94831480
#> u[18] -0.42631628
#> u[19] -1.11326722
#> u[20]  0.26048696
#> u[21]  0.04943089
#> u[22]  0.36841040
#> u[23] -0.45954833
#> u[24] -0.32354465
#> u[25] -0.81566715
#> u[26] -0.63627854
#> u[27]  0.41708341
#> u[28]  0.51890820
#> u[29] -0.48195210
#> u[30]  0.33018644

Random Effect Variance

model$refVar
#>      b_var    a_var
#> 1 1.051256 1.022215

STEP 4 Check Convergence via Plot MCMC

Trace Plot, Density Plot, ACF Plot, R-Hat Plot

model$plot
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

#> 
#> [[4]]

STEP 5 : Extract CV and MSE

Subarea

CV_sub <- (model$Est_sub$SD)/(model$Est_sub$Mean)*100
MSE_sub <- model$Est_sub$SD^2
summary(cbind(CV_sub,MSE_sub))
#>      CV_sub          MSE_sub        
#>  Min.   : 6.134   Min.   :0.003252  
#>  1st Qu.: 9.866   1st Qu.:0.007731  
#>  Median :14.473   Median :0.014417  
#>  Mean   :19.770   Mean   :0.019109  
#>  3rd Qu.:27.248   3rd Qu.:0.030579  
#>  Max.   :56.964   Max.   :0.046570

Area

CV_area <- (model$Est_area$SD)/(model$Est_area$Mean)*100
MSE_area <- model$Est_area$SD^2
summary(cbind(CV_area,MSE_area))
#>     CV_area          MSE_area       
#>  Min.   : 7.010   Min.   :0.003646  
#>  1st Qu.: 7.728   1st Qu.:0.008269  
#>  Median : 9.999   Median :0.014608  
#>  Mean   :13.334   Mean   :0.017324  
#>  3rd Qu.:17.224   3rd Qu.:0.026491  
#>  Max.   :26.684   Max.   :0.034845