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.8966052 0.07107975 0.7082006 0.8639809 0.9139528 0.9475523 0.9804602
#> mu[2]  0.7794031 0.13656777 0.4467804 0.7037122 0.8079012 0.8837839 0.9603528
#> mu[3]  0.9275998 0.05790586 0.7721361 0.9056191 0.9431986 0.9675139 0.9901330
#> mu[4]  0.9120231 0.06759491 0.7376775 0.8862149 0.9302282 0.9590194 0.9869932
#> mu[5]  0.9105315 0.06732917 0.7354999 0.8825938 0.9290603 0.9580125 0.9873871
#> mu[6]  0.8526129 0.10197630 0.5973356 0.8014539 0.8765094 0.9280774 0.9801539
#> mu[7]  0.5793726 0.21485192 0.1624214 0.4123838 0.5937530 0.7604050 0.9267392
#> mu[8]  0.9412993 0.04716914 0.8111250 0.9242626 0.9541895 0.9729693 0.9916134
#> mu[9]  0.8593551 0.10454571 0.5839037 0.8124180 0.8880255 0.9333864 0.9774820
#> mu[10] 0.9314645 0.05400160 0.7877229 0.9114519 0.9456599 0.9682042 0.9894242
#> mu[11] 0.9438586 0.04936934 0.8116530 0.9283336 0.9583609 0.9766212 0.9927813
#> mu[12] 0.8158019 0.11375998 0.5366798 0.7566435 0.8376464 0.9006296 0.9662720
#> mu[13] 0.7814371 0.14014643 0.4342847 0.6989076 0.8101832 0.8874073 0.9695780
#> mu[14] 0.9008555 0.07109109 0.7252457 0.8694219 0.9179534 0.9508816 0.9844362
#> mu[15] 0.9374843 0.05095017 0.8037618 0.9198150 0.9510894 0.9713687 0.9906655
#> mu[16] 0.8715305 0.08717359 0.6392453 0.8324225 0.8938421 0.9340584 0.9772492
#> mu[17] 0.6765117 0.17409439 0.2992777 0.5528822 0.6976144 0.8207857 0.9381034
#> mu[18] 0.9498699 0.04670209 0.8266845 0.9368566 0.9630710 0.9792442 0.9938966
#> mu[19] 0.9171951 0.05845025 0.7577866 0.8926341 0.9311227 0.9587601 0.9854903
#> mu[20] 0.9523706 0.04083128 0.8431754 0.9379825 0.9645005 0.9799125 0.9937263
#> mu[21] 0.9329083 0.06020401 0.7770252 0.9127939 0.9500979 0.9725677 0.9917846
#> mu[22] 0.8465910 0.10377398 0.5892556 0.7950692 0.8701691 0.9213412 0.9771723
#> mu[23] 0.9006290 0.08680138 0.6689496 0.8720526 0.9239582 0.9587375 0.9887530
#> mu[24] 0.7744177 0.12480097 0.4713509 0.7052920 0.7936544 0.8701107 0.9505603
#> mu[25] 0.8828875 0.07779094 0.6874571 0.8495217 0.9014733 0.9387413 0.9784153
#> mu[26] 0.8494934 0.09318886 0.6164362 0.8059916 0.8696476 0.9165160 0.9686744
#> mu[27] 0.8257561 0.13912349 0.4492664 0.7595428 0.8667129 0.9291140 0.9830235
#> mu[28] 0.9537760 0.04446682 0.8392131 0.9402412 0.9666283 0.9822266 0.9942711
#> mu[29] 0.8019630 0.12534767 0.4914379 0.7362147 0.8284638 0.8978379 0.9662236
#> mu[30] 0.8781884 0.11760693 0.5402954 0.8407962 0.9149768 0.9570393 0.9910066

Coefficient Estimation

example$coefficient
#>           Mean        SD       2.5%       25%       50%       75%     97.5%
#> b[0] 1.1077401 0.3341415  0.4589386 0.8792082 1.1119604 1.3298543 1.7704675
#> b[1] 0.2887721 0.2781500 -0.2306025 0.1010361 0.2768392 0.4733573 0.8421794
#> b[2] 0.3562527 0.1110676  0.1436468 0.2805527 0.3565517 0.4313557 0.5776276

Random Effect Variance Estimation

example$refvar
#> [1] 0.5144564

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.5794  
#>  1st Qu.:0.8404   1st Qu.:0.8310  
#>  Median :0.9941   Median :0.8897  
#>  Mean   :0.8696   Mean   :0.8661  
#>  3rd Qu.:1.0000   3rd Qu.:0.9305  
#>  Max.   :1.0000   Max.   :0.9538

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.001667  
#>  1st Qu.:0.007948   1st Qu.:0.003369  
#>  Median :0.018572   Median :0.006793  
#>  Mean   :0.021281   Mean   :0.009992  
#>  3rd Qu.:0.033778   3rd Qu.:0.013609  
#>  Max.   :0.053116   Max.   :0.046161

Comparing Relative Standard Error (RSE)

RSE=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta))
summary(RSE)
#>    RSE_direct       RSE_HBMEbeta   
#>  Min.   :  3.592   Min.   : 4.287  
#>  1st Qu.:  9.021   1st Qu.: 6.275  
#>  Median : 13.658   Median : 9.224  
#>  Mean   : 24.129   Mean   :11.168  
#>  3rd Qu.: 22.290   3rd Qu.:13.806  
#>  Max.   :184.551   Max.   :37.084

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.8997566 0.06606045 0.7301016 0.8709959 0.9167253 0.9468090 0.9807084
#> mu[2]  0.8032813 0.12920350 0.4847188 0.7315212 0.8337645 0.9018922 0.9671470
#> mu[3]  0.9368740 0.04993152 0.8074889 0.9190930 0.9491480 0.9709575 0.9907664
#> mu[4]  0.9185495 0.06521573 0.7455522 0.8934535 0.9371569 0.9625052 0.9879488
#> mu[5]  0.9293911 0.06076299 0.7711550 0.9087811 0.9456590 0.9688253 0.9910176
#> mu[6]  0.8486172 0.10132864 0.5910571 0.7969081 0.8701693 0.9245230 0.9762852
#> mu[7]  0.5436563 0.20620037 0.1670513 0.3802063 0.5485583 0.7062395 0.9055225
#> mu[8]  0.9431247 0.04457284 0.8177684 0.9271026 0.9554502 0.9735018 0.9917768
#> mu[9]  0.8534181 0.10124415 0.6012222 0.8056377 0.8771144 0.9273468 0.9754516
#> mu[10] 0.9348562 0.05370679 0.7986648 0.9172699 0.9495636 0.9701450 0.9899353
#> mu[11] 0.9426899 0.04899812 0.8068171 0.9264925 0.9567556 0.9748596 0.9921103
#> mu[12] 0.8566249 0.09941307 0.6175252 0.8126695 0.8796139 0.9292816 0.9768085
#> mu[13] 0.7809432 0.13799539 0.4515600 0.7022171 0.8090009 0.8853057 0.9701598
#> mu[14] 0.8810927 0.08334851 0.6726975 0.8445695 0.9009756 0.9404615 0.9816574
#> mu[15] 0.9446682 0.04465520 0.8268789 0.9279743 0.9567121 0.9751424 0.9920616
#> mu[16] 0.8549774 0.09708462 0.6030030 0.8097453 0.8803099 0.9251541 0.9732220
#> mu[17] 0.7512595 0.16326669 0.3596250 0.6553480 0.7874817 0.8816817 0.9627492
#> mu[18] 0.9507843 0.04437714 0.8325047 0.9377263 0.9638197 0.9791568 0.9938962
#> mu[19] 0.9050254 0.06752979 0.7283703 0.8762084 0.9212121 0.9522722 0.9833808
#> mu[20] 0.9543826 0.03916730 0.8570894 0.9413783 0.9646139 0.9804419 0.9941139
#> mu[21] 0.9474053 0.04775056 0.8164958 0.9324227 0.9608633 0.9787435 0.9938933
#> mu[22] 0.8630789 0.09507920 0.6103769 0.8184184 0.8866483 0.9310506 0.9791651
#> mu[23] 0.9016685 0.08126665 0.6859748 0.8732054 0.9235123 0.9566789 0.9890723
#> mu[24] 0.7679872 0.12661163 0.4587715 0.6940441 0.7904281 0.8631487 0.9499277
#> mu[25] 0.9016905 0.06954442 0.7204682 0.8698178 0.9200230 0.9520031 0.9837891
#> mu[26] 0.8294923 0.10344686 0.5681045 0.7759099 0.8506228 0.9048317 0.9648607
#> mu[27] 0.7858635 0.15400721 0.4002341 0.7017670 0.8257890 0.9048925 0.9764993
#> mu[28] 0.9643389 0.03399663 0.8709409 0.9549444 0.9747567 0.9855256 0.9959476
#> mu[29] 0.7929322 0.12675425 0.4741817 0.7253322 0.8158263 0.8858126 0.9675000
#> mu[30] 0.8823266 0.10688255 0.5896449 0.8430603 0.9145082 0.9543854 0.9911797

Coefficient Estimation

example_mix$coefficient
#>           Mean        SD        2.5%        25%       50%       75%     97.5%
#> b[0] 0.7861455 0.3454170  0.09789761 0.55228198 0.7847992 1.0143183 1.4669211
#> b[1] 0.2310974 0.2860504 -0.32152759 0.03704843 0.2239674 0.4215811 0.7911352
#> b[2] 0.3272224 0.1123814  0.10917474 0.24655387 0.3289483 0.4031225 0.5457641
#> b[3] 0.9468793 0.5172092 -0.09253922 0.61786987 0.9347353 1.2910048 1.9456825

Random Effect Variance Estimation

example_mix$refvar
#> [1] 0.5673716

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.5794  
#>  1st Qu.:0.8404   1st Qu.:0.8310  
#>  Median :0.9941   Median :0.8897  
#>  Mean   :0.8696   Mean   :0.8661  
#>  3rd Qu.:1.0000   3rd Qu.:0.9305  
#>  Max.   :1.0000   Max.   :0.9538

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.001156  
#>  1st Qu.:0.007948   1st Qu.:0.002591  
#>  Median :0.018572   Median :0.006776  
#>  Mean   :0.021281   Mean   :0.009522  
#>  3rd Qu.:0.033778   3rd Qu.:0.011243  
#>  Max.   :0.053116   Max.   :0.042519

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.525  
#>  1st Qu.:  9.021   1st Qu.: 5.433  
#>  Median : 13.658   Median : 9.236  
#>  Mean   : 24.129   Mean   :10.851  
#>  3rd Qu.: 22.290   3rd Qu.:12.382  
#>  Max.   :184.551   Max.   :37.928
note : you can use dataHBMEbetaNS for using dataset with non-sampled area