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.8436743 0.07735009 0.6550221 0.8047131 0.8580827 0.8987907 0.9506191
#> mu[2]  0.8458687 0.08072727 0.6402200 0.8043438 0.8608484 0.9046814 0.9559296
#> mu[3]  0.8435473 0.07798622 0.6481176 0.8040780 0.8591482 0.8986058 0.9507517
#> mu[4]  0.8230587 0.09084435 0.5897036 0.7781604 0.8399052 0.8869473 0.9487441
#> mu[5]  0.8477206 0.07935916 0.6528529 0.8064490 0.8638750 0.9051100 0.9570607
#> mu[6]  0.8401538 0.08230128 0.6374021 0.7991808 0.8561279 0.8996262 0.9528110
#> mu[7]  0.7652276 0.12149889 0.4530404 0.6998132 0.7897098 0.8548229 0.9316360
#> mu[8]  0.8496584 0.07731011 0.6603683 0.8112215 0.8642435 0.9052572 0.9547779
#> mu[9]  0.8434602 0.08458266 0.6343257 0.8031404 0.8595747 0.9039116 0.9579195
#> mu[10] 0.8448981 0.07804508 0.6548369 0.8066568 0.8587257 0.9013021 0.9524943
#> mu[11] 0.8473281 0.07736603 0.6611758 0.8062899 0.8624545 0.9026265 0.9547883
#> mu[12] 0.8362070 0.08356159 0.6375215 0.7933280 0.8517217 0.8964099 0.9525618
#> mu[13] 0.8457941 0.07937621 0.6461194 0.8041810 0.8607892 0.9032192 0.9527868
#> mu[14] 0.7947126 0.09863120 0.5500934 0.7428136 0.8107521 0.8654935 0.9370334
#> mu[15] 0.8433908 0.07839378 0.6546895 0.8054030 0.8576562 0.8994305 0.9541302
#> mu[16] 0.8459563 0.07627208 0.6606242 0.8053661 0.8604155 0.9004900 0.9534746
#> mu[17] 0.7611595 0.11897981 0.4804497 0.6939994 0.7832698 0.8507481 0.9270664
#> mu[18] 0.8414341 0.08525273 0.6326457 0.7985570 0.8586181 0.9027947 0.9573566
#> mu[19] 0.8596145 0.07196700 0.6821247 0.8231956 0.8721013 0.9116026 0.9583572
#> mu[20] 0.8353138 0.08195356 0.6394918 0.7909845 0.8497625 0.8957614 0.9504922
#> mu[21] 0.8458095 0.08035121 0.6498301 0.8044477 0.8613048 0.9035721 0.9562269
#> mu[22] 0.8021732 0.10199827 0.5468752 0.7520820 0.8219365 0.8732023 0.9460473
#> mu[23] 0.8453049 0.08308474 0.6374279 0.8024493 0.8609705 0.9040802 0.9569121
#> mu[24] 0.8095434 0.09625605 0.5676367 0.7569789 0.8270031 0.8805171 0.9431141
#> mu[25] 0.8464838 0.07536123 0.6607136 0.8076304 0.8604844 0.9016398 0.9516352
#> mu[26] 0.8253905 0.08528126 0.6253685 0.7790831 0.8406748 0.8875188 0.9475830
#> mu[27] 0.7729684 0.11589687 0.4774464 0.7114459 0.7955380 0.8575615 0.9343527
#> mu[28] 0.8472263 0.07982049 0.6520825 0.8082864 0.8627476 0.9042278 0.9578021
#> mu[29] 0.7931947 0.10378910 0.5377861 0.7393610 0.8114102 0.8685412 0.9409732
#> mu[30] 0.8413018 0.09936571 0.5880837 0.7936364 0.8620292 0.9126177 0.9683289

Coefficient Estimation

example$coefficient
#>            Mean        SD       2.5%         25%        50%       75%     97.5%
#> b[0] 1.65221293 0.2372866  1.1964501  1.48983226 1.65441969 1.8133777 2.1162262
#> b[1] 0.07899627 0.2105544 -0.3414484 -0.06060007 0.07837086 0.2240071 0.4824158
#> b[2] 0.02813679 0.1247927 -0.2127049 -0.05458379 0.02709273 0.1096369 0.2770789

Random Effect Variance Estimation

example$refvar
#> [1] 0.3254072

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.7612  
#>  1st Qu.:9.99e-01   1st Qu.:0.8236  
#>  Median :1.00e+00   Median :0.8434  
#>  Mean   :8.64e-01   Mean   :0.8296  
#>  3rd Qu.:1.00e+00   3rd Qu.:0.8459  
#>  Max.   :1.00e+00   Max.   :0.8596

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.005179  
#>  1st Qu.:3.306e-05   1st Qu.:0.006105  
#>  Median :5.868e-04   Median :0.006745  
#>  Mean   :7.635e-03   Mean   :0.007820  
#>  3rd Qu.:8.967e-03   3rd Qu.:0.009012  
#>  Max.   :5.569e-02   Max.   :0.014762

Comparing Relative Standard Error (RSE)

RSE=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta))
summary(RSE)
#>    RSE_direct       RSE_HBMEbeta   
#>  Min.   :      0   Min.   : 8.372  
#>  1st Qu.:      1   1st Qu.: 9.258  
#>  Median :      2   Median : 9.804  
#>  Mean   : 175957   Mean   :10.602  
#>  3rd Qu.:     11   3rd Qu.:11.618  
#>  Max.   :5247828   Max.   :15.877

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.8432703 0.07811595 0.6525169 0.8019284 0.8569977 0.9009389 0.9547623
#> mu[2]  0.8609552 0.07845569 0.6644535 0.8229216 0.8782319 0.9171474 0.9626691
#> mu[3]  0.8582950 0.07667694 0.6651080 0.8208157 0.8736234 0.9125798 0.9596885
#> mu[4]  0.8404908 0.08511013 0.6303771 0.7966421 0.8583990 0.9014441 0.9562038
#> mu[5]  0.8795967 0.07549319 0.6901568 0.8461060 0.8979673 0.9333064 0.9736563
#> mu[6]  0.8373425 0.08495690 0.6291810 0.7943367 0.8547155 0.8979123 0.9554784
#> mu[7]  0.7200354 0.13256484 0.4105285 0.6471242 0.7448583 0.8192456 0.9110731
#> mu[8]  0.8524130 0.07690692 0.6588372 0.8122020 0.8665680 0.9075190 0.9586833
#> mu[9]  0.8239808 0.09114622 0.6089862 0.7744438 0.8406071 0.8928057 0.9532461
#> mu[10] 0.8453634 0.08036090 0.6457588 0.8054194 0.8596590 0.9035337 0.9552729
#> mu[11] 0.8357960 0.08607041 0.6189846 0.7912129 0.8511413 0.8974571 0.9539811
#> mu[12] 0.8734183 0.07741717 0.6706549 0.8403516 0.8911870 0.9274080 0.9686387
#> mu[13] 0.8358400 0.08611697 0.6375670 0.7909366 0.8502625 0.8986072 0.9560941
#> mu[14] 0.7687870 0.10712688 0.5216712 0.7057305 0.7854219 0.8474828 0.9269003
#> mu[15] 0.8545845 0.07475907 0.6785083 0.8139381 0.8684851 0.9097672 0.9597363
#> mu[16] 0.8172271 0.09052377 0.5974310 0.7696432 0.8330749 0.8824850 0.9481497
#> mu[17] 0.8170406 0.11539631 0.5194404 0.7619438 0.8438489 0.9017401 0.9587462
#> mu[18] 0.8462358 0.08442389 0.6387833 0.8063511 0.8625411 0.9053067 0.9604072
#> mu[19] 0.8341133 0.08326185 0.6330237 0.7921379 0.8468643 0.8939717 0.9509302
#> mu[20] 0.8467810 0.08055947 0.6440359 0.8051035 0.8625578 0.9056571 0.9578906
#> mu[21] 0.8803107 0.07130247 0.7015990 0.8477693 0.8961253 0.9312760 0.9725240
#> mu[22] 0.8191636 0.10147909 0.5698528 0.7693331 0.8399783 0.8937765 0.9539431
#> mu[23] 0.8415033 0.08503327 0.6313088 0.7983574 0.8570644 0.9014187 0.9576047
#> mu[24] 0.7914729 0.10560461 0.5294044 0.7338186 0.8099070 0.8678384 0.9390144
#> mu[25] 0.8714204 0.07067913 0.6945880 0.8350759 0.8863060 0.9230522 0.9653513
#> mu[26] 0.7909124 0.09821208 0.5655852 0.7328084 0.8058891 0.8636746 0.9357595
#> mu[27] 0.7339148 0.12461890 0.4437852 0.6654220 0.7528177 0.8254328 0.9182404
#> mu[28] 0.8745893 0.07245089 0.6933386 0.8391909 0.8908521 0.9253647 0.9678016
#> mu[29] 0.7713764 0.11240701 0.5031873 0.7094338 0.7910266 0.8547998 0.9305015
#> mu[30] 0.8394143 0.10085203 0.5904817 0.7905107 0.8616663 0.9117590 0.9715047

Coefficient Estimation

example_mix$coefficient
#>            Mean        SD       2.5%        25%        50%       75%     97.5%
#> b[0] 1.28787585 0.2931234  0.7272828  1.0945180 1.27868589 1.4896316 1.8713818
#> b[1] 0.04349063 0.2188143 -0.3986533 -0.1027667 0.04044306 0.1836323 0.4829281
#> b[2] 0.02825011 0.1257353 -0.2167400 -0.0597437 0.03009467 0.1155645 0.2737182
#> b[3] 0.80327699 0.4695010 -0.1220506  0.4877225 0.81588628 1.1202939 1.6742491

Random Effect Variance Estimation

example_mix$refvar
#> [1] 0.3216946

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.7612  
#>  1st Qu.:9.99e-01   1st Qu.:0.8236  
#>  Median :1.00e+00   Median :0.8434  
#>  Mean   :8.64e-01   Mean   :0.8296  
#>  3rd Qu.:1.00e+00   3rd Qu.:0.8459  
#>  Max.   :1.00e+00   Max.   :0.8596

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.004996  
#>  1st Qu.:3.306e-05   1st Qu.:0.006021  
#>  Median :5.868e-04   Median :0.007224  
#>  Mean   :7.635e-03   Mean   :0.008283  
#>  3rd Qu.:8.967e-03   3rd Qu.:0.010040  
#>  Max.   :5.569e-02   Max.   :0.017573

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.100  
#>  1st Qu.:      1   1st Qu.: 9.045  
#>  Median :      2   Median :10.116  
#>  Mean   : 175957   Mean   :10.910  
#>  3rd Qu.:     11   3rd Qu.:12.295  
#>  Max.   :5247828   Max.   :18.411
note : you can use dataHBMEbetaNS for using dataset with non-sampled area