An Aplication to SAE HB ME under Beta Distribution On Sample Data

case 1 : Auxiliary variables only contains variable with error in aux variable

Load Package and Load Data Set

library(saeHB.ME.beta)
data("dataHBMEbeta")

Fitting HB Model

example <- meHBbeta(Y~x1+x2, var.x = c("v.x1","v.x2"),
                   iter.update = 3, iter.mcmc = 10000,
                   thin = 3, burn.in = 1000, data = dataHBMEbeta)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 126
#>    Total graph size: 623
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 126
#>    Total graph size: 623
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 126
#>    Total graph size: 623
#> 
#> Initializing model

Ekstract Mean Estimation

Small Area Mean Estimation

example$Est
#>             mean         sd      2.5%       25%       50%       75%     97.5%
#> mu[1]  0.8416853 0.07503160 0.6663454 0.8021204 0.8539008 0.8957161 0.9512338
#> mu[2]  0.8463504 0.08074311 0.6518505 0.8053343 0.8611295 0.9047494 0.9569837
#> mu[3]  0.8439023 0.07723618 0.6550671 0.8034900 0.8586887 0.8991897 0.9504220
#> mu[4]  0.8235074 0.08840573 0.6013875 0.7798147 0.8405103 0.8869287 0.9461643
#> mu[5]  0.8459879 0.07871477 0.6571124 0.8041918 0.8598413 0.9037009 0.9549178
#> mu[6]  0.8381432 0.08121686 0.6410625 0.7947501 0.8526293 0.8967996 0.9529408
#> mu[7]  0.7641579 0.12265984 0.4619847 0.7009210 0.7881012 0.8536115 0.9314142
#> mu[8]  0.8478323 0.07618529 0.6669401 0.8065519 0.8602505 0.9041729 0.9564811
#> mu[9]  0.8426264 0.08422705 0.6367777 0.7998750 0.8583432 0.9037904 0.9593601
#> mu[10] 0.8437598 0.07836709 0.6511212 0.8035820 0.8577890 0.9000326 0.9516863
#> mu[11] 0.8444968 0.08173540 0.6435679 0.8035870 0.8609307 0.9034221 0.9542632
#> mu[12] 0.8358571 0.08391067 0.6276648 0.7947853 0.8519023 0.8950722 0.9510740
#> mu[13] 0.8467580 0.07764394 0.6620718 0.8042278 0.8606256 0.9032737 0.9554820
#> mu[14] 0.7951709 0.09998640 0.5523233 0.7409157 0.8117109 0.8687643 0.9377560
#> mu[15] 0.8404791 0.07985322 0.6497380 0.7985050 0.8539948 0.8975110 0.9516097
#> mu[16] 0.8452365 0.07863805 0.6510475 0.8039895 0.8592398 0.9031062 0.9542445
#> mu[17] 0.7635485 0.11793018 0.4725584 0.6988790 0.7870598 0.8496561 0.9270834
#> mu[18] 0.8419619 0.08750373 0.6243736 0.8007166 0.8592568 0.9030673 0.9581212
#> mu[19] 0.8558909 0.07420725 0.6796130 0.8189064 0.8705806 0.9096122 0.9545350
#> mu[20] 0.8333575 0.08310627 0.6308405 0.7912695 0.8483918 0.8938964 0.9500743
#> mu[21] 0.8433135 0.08358542 0.6392663 0.8045128 0.8603984 0.9035820 0.9538114
#> mu[22] 0.7981174 0.10398068 0.5472636 0.7430299 0.8170609 0.8743005 0.9451216
#> mu[23] 0.8425901 0.08423087 0.6275140 0.7988084 0.8600944 0.9035623 0.9564442
#> mu[24] 0.8036584 0.10240068 0.5453360 0.7494526 0.8223033 0.8773301 0.9435467
#> mu[25] 0.8437857 0.07852212 0.6534682 0.8022027 0.8591593 0.9019915 0.9530249
#> mu[26] 0.8248564 0.08662463 0.6053740 0.7792868 0.8415287 0.8880432 0.9457300
#> mu[27] 0.7754028 0.11291888 0.4985139 0.7185032 0.7958109 0.8567496 0.9330853
#> mu[28] 0.8456545 0.08025105 0.6439812 0.8069383 0.8608143 0.9028653 0.9553957
#> mu[29] 0.7903936 0.10536840 0.5342441 0.7352144 0.8110449 0.8683767 0.9377924
#> mu[30] 0.8442913 0.09857787 0.5885866 0.7996503 0.8684811 0.9131264 0.9695456

Coefficient Estimation

example$coefficient
#>            Mean        SD       2.5%         25%        50%       75%     97.5%
#> b[0] 1.65485883 0.2397138  1.1721340  1.49785315 1.65687373 1.8181941 2.1259435
#> b[1] 0.06602099 0.2130744 -0.3334683 -0.08100469 0.06366662 0.2099927 0.4784728
#> b[2] 0.02207541 0.1244913 -0.2189939 -0.06144032 0.02164405 0.1048711 0.2743714

Random Effect Variance Estimation

example$refvar
#> [1] 0.3233804

Extract MSE

MSE_HBMEbeta=example$Est$sd^2

Extract RSE

RSE_HBMEbeta=sqrt(MSE_HBMEbeta)/example$Est$mean*100

You can compare with Direct Estimation

Extract Direct Estimation

Y_direct=dataHBMEbeta[,1]
MSE_direct=dataHBMEbeta[,6]
RSE_direct=sqrt(MSE_direct)/Y_direct*100

Comparing Y

Y_HBMEbeta=example$Est$mean
Y=as.data.frame(cbind(Y_direct,Y_HBMEbeta))
summary(Y)
#>     Y_direct          Y_HBMEbeta    
#>  Min.   :1.50e-06   Min.   :0.7635  
#>  1st Qu.:9.99e-01   1st Qu.:0.8238  
#>  Median :1.00e+00   Median :0.8423  
#>  Mean   :8.64e-01   Mean   :0.8284  
#>  3rd Qu.:1.00e+00   3rd Qu.:0.8444  
#>  Max.   :1.00e+00   Max.   :0.8559

Comparing Mean Squared Error (MSE)

MSE=as.data.frame(cbind(MSE_direct,MSE_HBMEbeta))
summary(MSE)
#>    MSE_direct         MSE_HBMEbeta     
#>  Min.   :1.500e-07   Min.   :0.005507  
#>  1st Qu.:3.306e-05   1st Qu.:0.006187  
#>  Median :5.868e-04   Median :0.006947  
#>  Mean   :7.635e-03   Mean   :0.007938  
#>  3rd Qu.:8.967e-03   3rd Qu.:0.009242  
#>  Max.   :5.569e-02   Max.   :0.015045

Comparing Relative Standard Error (RSE)

RSE=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta))
summary(RSE)
#>    RSE_direct       RSE_HBMEbeta   
#>  Min.   :      0   Min.   : 8.670  
#>  1st Qu.:      1   1st Qu.: 9.305  
#>  Median :      2   Median : 9.942  
#>  Mean   : 175957   Mean   :10.698  
#>  3rd Qu.:     11   3rd Qu.:11.441  
#>  Max.   :5247828   Max.   :16.052

case 2: Auxiliary variables contains variable with error and without error

Fitting HB Model

example_mix <- meHBbeta(Y~x1+x2+x3, var.x = c("v.x1","v.x2"),
                   iter.update = 3, iter.mcmc = 10000,
                   thin = 3, burn.in = 1000, data = dataHBMEbeta)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 127
#>    Total graph size: 717
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 127
#>    Total graph size: 717
#> 
#> Initializing model
#> 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 30
#>    Unobserved stochastic nodes: 127
#>    Total graph size: 717
#> 
#> Initializing model

Ekstract Mean Estimation

Small Area Mean Estimation

example_mix$Est
#>             mean         sd      2.5%       25%       50%       75%     97.5%
#> mu[1]  0.8438472 0.07722785 0.6564269 0.8042641 0.8577980 0.8992938 0.9501498
#> mu[2]  0.8596662 0.07557709 0.6720641 0.8220169 0.8742556 0.9140819 0.9627313
#> mu[3]  0.8547806 0.07612989 0.6633313 0.8153258 0.8690362 0.9092639 0.9604253
#> mu[4]  0.8370814 0.08439967 0.6322903 0.7934238 0.8527133 0.8983265 0.9545385
#> mu[5]  0.8782910 0.07386179 0.6834716 0.8472002 0.8946201 0.9293734 0.9704964
#> mu[6]  0.8377267 0.08138082 0.6426118 0.7955398 0.8518593 0.8979894 0.9526109
#> mu[7]  0.7247963 0.13053092 0.4234166 0.6509727 0.7474068 0.8196420 0.9169745
#> mu[8]  0.8488618 0.07823629 0.6593455 0.8087150 0.8641079 0.9054945 0.9579902
#> mu[9]  0.8244094 0.09328407 0.5984369 0.7774757 0.8433142 0.8912295 0.9512582
#> mu[10] 0.8445639 0.08054587 0.6527087 0.8038852 0.8589208 0.9024319 0.9557523
#> mu[11] 0.8387194 0.08236559 0.6321406 0.7950637 0.8525188 0.8991819 0.9552008
#> mu[12] 0.8729207 0.07410988 0.6828150 0.8400825 0.8887964 0.9251204 0.9682803
#> mu[13] 0.8403636 0.08041833 0.6425800 0.7990467 0.8552089 0.8981045 0.9554823
#> mu[14] 0.7719146 0.10374423 0.5191747 0.7164453 0.7909995 0.8465631 0.9245720
#> mu[15] 0.8543727 0.07433915 0.6753214 0.8160795 0.8680057 0.9086660 0.9568454
#> mu[16] 0.8187067 0.08823105 0.6162447 0.7696732 0.8328894 0.8839035 0.9466283
#> mu[17] 0.8196221 0.10742766 0.5431860 0.7701667 0.8435198 0.8967238 0.9555676
#> mu[18] 0.8465896 0.08407211 0.6426450 0.8037892 0.8629280 0.9060966 0.9605002
#> mu[19] 0.8341665 0.08116591 0.6390407 0.7904435 0.8484006 0.8938307 0.9489397
#> mu[20] 0.8423202 0.08403403 0.6385869 0.8003616 0.8595833 0.9043118 0.9537410
#> mu[21] 0.8768050 0.07499803 0.6905004 0.8411294 0.8935940 0.9301278 0.9713527
#> mu[22] 0.8207075 0.09980845 0.5702020 0.7693380 0.8423226 0.8935895 0.9540973
#> mu[23] 0.8409029 0.08530114 0.6303807 0.7987527 0.8567405 0.9009653 0.9575608
#> mu[24] 0.7943102 0.10319248 0.5435801 0.7385560 0.8151732 0.8688209 0.9415578
#> mu[25] 0.8682000 0.07214065 0.6828236 0.8338007 0.8823259 0.9201328 0.9646815
#> mu[26] 0.7945181 0.09280405 0.5733252 0.7447108 0.8084101 0.8635545 0.9317437
#> mu[27] 0.7341388 0.12619341 0.4257459 0.6677606 0.7557034 0.8275831 0.9156738
#> mu[28] 0.8723977 0.07573922 0.6785626 0.8404803 0.8902613 0.9250844 0.9666099
#> mu[29] 0.7750740 0.10883178 0.5151701 0.7169573 0.7930869 0.8546596 0.9313368
#> mu[30] 0.8345906 0.10640752 0.5588260 0.7850818 0.8569105 0.9091802 0.9717777

Coefficient Estimation

example_mix$coefficient
#>             Mean        SD       2.5%         25%        50%        75%
#> b[0] 1.321988106 0.2966870  0.7379480  1.12115821 1.32322839 1.52435384
#> b[1] 0.046397710 0.2236702 -0.3877839 -0.10260848 0.04411156 0.19258886
#> b[2] 0.008715332 0.1296287 -0.2407654 -0.07631908 0.01056678 0.09321969
#> b[3] 0.750751276 0.4451693 -0.1283277  0.45004538 0.74731428 1.05556146
#>          97.5%
#> b[0] 1.9026168
#> b[1] 0.4931409
#> b[2] 0.2645530
#> b[3] 1.6151043

Random Effect Variance Estimation

example_mix$refvar
#> [1] 0.3065507

Extract MSE

MSE_HBMEbeta_mix=example_mix$Est$sd^2

Extract RSE

RSE_HBMEbeta_mix=sqrt(MSE_HBMEbeta_mix)/example_mix$Est$mean*100

You can compare with Direct Estimation

Comparing Y

Y_HBMEbeta_mix=example_mix$Est$mean
Y_mix=as.data.frame(cbind(Y_direct,Y_HBMEbeta_mix))
summary(Y)
#>     Y_direct          Y_HBMEbeta    
#>  Min.   :1.50e-06   Min.   :0.7635  
#>  1st Qu.:9.99e-01   1st Qu.:0.8238  
#>  Median :1.00e+00   Median :0.8423  
#>  Mean   :8.64e-01   Mean   :0.8284  
#>  3rd Qu.:1.00e+00   3rd Qu.:0.8444  
#>  Max.   :1.00e+00   Max.   :0.8559

Comparing Mean Squared Error (MSE)

MSE_mix=as.data.frame(cbind(MSE_direct,MSE_HBMEbeta_mix))
summary(MSE_mix)
#>    MSE_direct        MSE_HBMEbeta_mix  
#>  Min.   :1.500e-07   Min.   :0.005204  
#>  1st Qu.:3.306e-05   1st Qu.:0.005838  
#>  Median :5.868e-04   Median :0.006923  
#>  Mean   :7.635e-03   Mean   :0.008075  
#>  3rd Qu.:8.967e-03   3rd Qu.:0.009647  
#>  Max.   :5.569e-02   Max.   :0.017038

Comparing Relative Standard Error (RSE)

RSE_mix=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta_mix))
summary(RSE_mix)
#>    RSE_direct      RSE_HBMEbeta_mix
#>  Min.   :      0   Min.   : 8.309  
#>  1st Qu.:      1   1st Qu.: 8.968  
#>  Median :      2   Median : 9.876  
#>  Mean   : 175957   Mean   :10.773  
#>  3rd Qu.:     11   3rd Qu.:12.041  
#>  Max.   :5247828   Max.   :18.009
note : you can use dataHBMEbetaNS for using dataset with non-sampled area