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.9014769 0.06936093 0.7257684 0.8712883 0.9169440 0.9504468 0.9833955
#> mu[2]  0.7775968 0.13950892 0.4513797 0.7005513 0.8074728 0.8816158 0.9657572
#> mu[3]  0.9348819 0.05563493 0.7766253 0.9158085 0.9500579 0.9721796 0.9920648
#> mu[4]  0.9188865 0.06285304 0.7546586 0.8916753 0.9356428 0.9635376 0.9892233
#> mu[5]  0.9157478 0.06786241 0.7288942 0.8919420 0.9336556 0.9617780 0.9883086
#> mu[6]  0.8526465 0.10555796 0.5896977 0.7975672 0.8791645 0.9301823 0.9826674
#> mu[7]  0.5458412 0.22093259 0.1450091 0.3699476 0.5472179 0.7222873 0.9273907
#> mu[8]  0.9467127 0.04730011 0.8235115 0.9315511 0.9603042 0.9774117 0.9931057
#> mu[9]  0.8666614 0.10384973 0.5945812 0.8210583 0.8949309 0.9395383 0.9816303
#> mu[10] 0.9366700 0.04960332 0.8050297 0.9159127 0.9498583 0.9714215 0.9909851
#> mu[11] 0.9490847 0.04797719 0.8107559 0.9361210 0.9635982 0.9792311 0.9939132
#> mu[12] 0.8170334 0.11422137 0.5422671 0.7568341 0.8417220 0.9028191 0.9660620
#> mu[13] 0.7725224 0.14360364 0.4342217 0.6839850 0.8000838 0.8842707 0.9714839
#> mu[14] 0.9028851 0.07205969 0.7146339 0.8715317 0.9220976 0.9545249 0.9853101
#> mu[15] 0.9452388 0.04312261 0.8335708 0.9288603 0.9567450 0.9745581 0.9931357
#> mu[16] 0.8756351 0.08611002 0.6515605 0.8363840 0.8956041 0.9387777 0.9797624
#> mu[17] 0.6557450 0.17674322 0.2813359 0.5319164 0.6746834 0.7954139 0.9329857
#> mu[18] 0.9547193 0.04380568 0.8369104 0.9428281 0.9671736 0.9817878 0.9949269
#> mu[19] 0.9221948 0.05842717 0.7661582 0.8997354 0.9383235 0.9627091 0.9878729
#> mu[20] 0.9570032 0.03762601 0.8576019 0.9460195 0.9675407 0.9812432 0.9948962
#> mu[21] 0.9389192 0.05726275 0.7792374 0.9216674 0.9561066 0.9762862 0.9936317
#> mu[22] 0.8491903 0.10426431 0.5807587 0.7984674 0.8739790 0.9241017 0.9783677
#> mu[23] 0.9051896 0.08250940 0.6762629 0.8779907 0.9274266 0.9613380 0.9904097
#> mu[24] 0.7722720 0.12811637 0.4664481 0.7002112 0.7948103 0.8681949 0.9536831
#> mu[25] 0.8865441 0.07620690 0.6927788 0.8479387 0.9032242 0.9422322 0.9801423
#> mu[26] 0.8489853 0.09643323 0.6058780 0.8018204 0.8711678 0.9187309 0.9729361
#> mu[27] 0.8140757 0.15547147 0.3948924 0.7441326 0.8622698 0.9313268 0.9825537
#> mu[28] 0.9599913 0.03960918 0.8548818 0.9491588 0.9716947 0.9848246 0.9961929
#> mu[29] 0.8035007 0.12402434 0.4943832 0.7419471 0.8293812 0.8951268 0.9675836
#> mu[30] 0.8901128 0.10501580 0.6015574 0.8547233 0.9222994 0.9617059 0.9926407

Coefficient Estimation

example$coefficient
#>           Mean        SD       2.5%       25%       50%       75%     97.5%
#> b[0] 1.0401620 0.3219482  0.4076461 0.8241162 1.0390895 1.2533070 1.6679399
#> b[1] 0.2917453 0.2789920 -0.2376792 0.1008333 0.2864888 0.4744689 0.8579734
#> b[2] 0.3983480 0.1120675  0.1861012 0.3205009 0.3958005 0.4716675 0.6276406

Random Effect Variance Estimation

example$refvar
#> [1] 0.5370229

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.   :0.1110   Min.   :0.5458  
#>  1st Qu.:0.8404   1st Qu.:0.8250  
#>  Median :0.9941   Median :0.8958  
#>  Mean   :0.8696   Mean   :0.8673  
#>  3rd Qu.:1.0000   3rd Qu.:0.9362  
#>  Max.   :1.0000   Max.   :0.9600

Comparing Mean Squared Error (MSE)

MSE=as.data.frame(cbind(MSE_direct,MSE_HBMEbeta))
summary(MSE)
#>    MSE_direct        MSE_HBMEbeta     
#>  Min.   :0.001290   Min.   :0.001416  
#>  1st Qu.:0.007948   1st Qu.:0.003141  
#>  Median :0.018572   Median :0.006308  
#>  Mean   :0.021281   Mean   :0.010147  
#>  3rd Qu.:0.033778   3rd Qu.:0.012571  
#>  Max.   :0.053116   Max.   :0.048811

Comparing Relative Standard Error (RSE)

RSE=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta))
summary(RSE)
#>    RSE_direct       RSE_HBMEbeta   
#>  Min.   :  3.592   Min.   : 3.932  
#>  1st Qu.:  9.021   1st Qu.: 5.988  
#>  Median : 13.658   Median : 8.856  
#>  Mean   : 24.129   Mean   :11.242  
#>  3rd Qu.: 22.290   3rd Qu.:13.580  
#>  Max.   :184.551   Max.   :40.476

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.8970133 0.06810979 0.7323970 0.8648256 0.9126631 0.9455287 0.9811981
#> mu[2]  0.8070847 0.12408488 0.4964117 0.7392668 0.8366825 0.8998387 0.9674536
#> mu[3]  0.9320252 0.05276613 0.7911582 0.9119322 0.9460987 0.9684988 0.9901196
#> mu[4]  0.9144207 0.06533837 0.7405387 0.8874852 0.9313511 0.9599225 0.9881461
#> mu[5]  0.9285016 0.05643053 0.7746122 0.9055382 0.9445479 0.9675252 0.9894984
#> mu[6]  0.8496556 0.10276880 0.5842047 0.8031884 0.8724716 0.9234793 0.9754258
#> mu[7]  0.5647885 0.20937895 0.1688913 0.4058302 0.5765506 0.7359121 0.9087996
#> mu[8]  0.9385665 0.05080471 0.7955667 0.9222046 0.9522873 0.9719726 0.9912799
#> mu[9]  0.8516288 0.10382330 0.5868327 0.8014782 0.8761902 0.9276669 0.9770651
#> mu[10] 0.9286563 0.05506942 0.7829271 0.9078152 0.9436160 0.9666592 0.9891093
#> mu[11] 0.9375954 0.05130428 0.7995507 0.9167948 0.9531937 0.9735819 0.9918736
#> mu[12] 0.8600248 0.09710993 0.6022564 0.8167239 0.8834003 0.9302885 0.9763832
#> mu[13] 0.7809472 0.13372358 0.4546189 0.7089040 0.8066999 0.8802375 0.9658072
#> mu[14] 0.8764190 0.08430121 0.6642307 0.8366529 0.8950809 0.9369373 0.9806060
#> mu[15] 0.9398863 0.04540734 0.8191448 0.9206904 0.9526822 0.9715086 0.9906909
#> mu[16] 0.8521331 0.09609209 0.6076881 0.8034450 0.8749071 0.9218514 0.9730482
#> mu[17] 0.7580418 0.16184038 0.3762273 0.6631549 0.7957762 0.8847765 0.9641101
#> mu[18] 0.9471518 0.04404329 0.8240718 0.9306686 0.9594113 0.9776933 0.9930924
#> mu[19] 0.8981253 0.07021686 0.7135316 0.8657403 0.9158889 0.9483284 0.9828340
#> mu[20] 0.9523858 0.04064267 0.8482031 0.9388682 0.9634449 0.9790414 0.9935240
#> mu[21] 0.9423036 0.05345404 0.8079626 0.9252097 0.9576179 0.9765132 0.9929560
#> mu[22] 0.8599000 0.09675708 0.6174361 0.8150042 0.8830981 0.9301439 0.9780386
#> mu[23] 0.8979243 0.08272481 0.6840770 0.8666935 0.9197052 0.9547404 0.9873338
#> mu[24] 0.7702797 0.12624023 0.4685095 0.6983110 0.7940664 0.8632374 0.9504542
#> mu[25] 0.9015438 0.06783134 0.7184096 0.8721805 0.9183943 0.9488643 0.9830123
#> mu[26] 0.8286458 0.10063687 0.5774598 0.7751328 0.8499142 0.9044448 0.9631409
#> mu[27] 0.7837436 0.15410418 0.4021578 0.7035826 0.8233122 0.9012793 0.9742931
#> mu[28] 0.9604387 0.03689736 0.8701376 0.9489823 0.9709026 0.9847133 0.9953791
#> mu[29] 0.7949657 0.12035117 0.5032860 0.7294183 0.8203461 0.8832584 0.9574039
#> mu[30] 0.8785039 0.10589779 0.5984718 0.8345568 0.9086267 0.9535693 0.9891460

Coefficient Estimation

example_mix$coefficient
#>           Mean        SD        2.5%        25%       50%       75%     97.5%
#> b[0] 0.8228225 0.3628065  0.11128702 0.57418580 0.8210162 1.0749201 1.5145888
#> b[1] 0.2422992 0.2818270 -0.30842846 0.05436379 0.2327066 0.4262353 0.8132562
#> b[2] 0.2996899 0.1141810  0.08317369 0.22068957 0.2997649 0.3770250 0.5255571
#> b[3] 0.9411094 0.5220184 -0.07962099 0.59345354 0.9432683 1.3014204 1.9184867

Random Effect Variance Estimation

example_mix$refvar
#> [1] 0.5476974

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.   :0.1110   Min.   :0.5458  
#>  1st Qu.:0.8404   1st Qu.:0.8250  
#>  Median :0.9941   Median :0.8958  
#>  Mean   :0.8696   Mean   :0.8673  
#>  3rd Qu.:1.0000   3rd Qu.:0.9362  
#>  Max.   :1.0000   Max.   :0.9600

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.   :0.001290   Min.   :0.001361  
#>  1st Qu.:0.007948   1st Qu.:0.002901  
#>  Median :0.018572   Median :0.006975  
#>  Mean   :0.021281   Mean   :0.009489  
#>  3rd Qu.:0.033778   3rd Qu.:0.011106  
#>  Max.   :0.053116   Max.   :0.043840

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.   :  3.592   Min.   : 3.842  
#>  1st Qu.:  9.021   1st Qu.: 5.737  
#>  Median : 13.658   Median : 9.416  
#>  Mean   : 24.129   Mean   :10.838  
#>  3rd Qu.: 22.290   3rd Qu.:12.180  
#>  Max.   :184.551   Max.   :37.072
note : you can use dataHBMEbetaNS for using dataset with non-sampled area