saeHB-Twofold-Beta

STEP 1 Load package and load the data

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

STEP 2 Data Exploration

dataBeta$CV <- sqrt(dataBeta$vardir)/dataBeta$y
explore(y~X1+X2, CV = "CV", data = dataBeta, normality = TRUE)
#>                   y         X1         X2
#> Min.    0.002826007 0.04205953 0.02461368
#> 1st Qu. 0.677417203 0.34818477 0.20900053
#> Median  0.986658040 0.58338771 0.39409355
#> Mean    0.786269096 0.57240016 0.43864087
#> 3rd Qu. 0.999116512 0.85933934 0.73765722
#> Max.    1.000000000 0.99426978 0.96302423
#> NA      0.000000000 0.00000000 0.00000000
#> 
#> Normality test for y :
#> Decision: Data do NOT follow normal distribution, with p.value = 0 < 0.05

STEPS 3 Fitting Twofold HB Beta Model

model <- betaTF(y~X1+X2,area="codearea",weight="w",iter.mcmc = 10000, burn.in = 3000, iter.update = 5, thin = 10, data=dataBeta)

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 Result for Area

Area Mean Estimation

model$Est_area
#>         Mean         SD      2.5%       25%       50%       75%     97.5%
#> 1  1.2195835 0.08136864 1.0220924 1.1832158 1.2341439 1.2784699 1.3247949
#> 2  0.9688439 0.10752455 0.7361556 0.9078791 0.9827910 1.0476230 1.1379279
#> 3  1.2755892 0.10623761 1.0228377 1.2220985 1.2962507 1.3524717 1.4163335
#> 4  0.9691133 0.06863815 0.7900602 0.9348939 0.9816288 1.0188424 1.0578658
#> 5  0.9770920 0.16135194 0.6135998 0.8775442 0.9857523 1.0929231 1.2593776
#> 6  1.4818411 0.13679931 1.1593274 1.4106433 1.5062895 1.5822590 1.6692599
#> 7  0.8325981 0.17642277 0.4576188 0.7162542 0.8405524 0.9428438 1.1751190
#> 8  0.8398394 0.15075759 0.5359534 0.7373188 0.8467865 0.9525999 1.0907456
#> 9  0.7218695 0.18920000 0.3600545 0.5857106 0.7223206 0.8560675 1.0888686
#> 10 0.7987668 0.06318449 0.6587559 0.7668503 0.8087824 0.8448400 0.8851057

Area Random Effect

model$area_randeff
#>            f_mean
#> f[1]   0.29608768
#> f[2]   0.01303191
#> f[3]   0.46917624
#> f[4]   0.73444469
#> f[5]  -0.37974895
#> f[6]   0.54180458
#> f[7]  -0.69091732
#> f[8]  -0.40666414
#> f[9]  -0.94530268
#> f[10]  0.38830953

Calculate Area Relative Standard Error (RSE) or CV and Mean Squared Error (MSE)

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.   : 6.672   Min.   :0.003992  
#>  1st Qu.: 8.015   1st Qu.:0.007787  
#>  Median :10.165   Median :0.015138  
#>  Mean   :13.219   Mean   :0.017257  
#>  3rd Qu.:17.591   3rd Qu.:0.025208  
#>  Max.   :26.210   Max.   :0.035797

STEP 6 Extract Result for Subarea

Subarea Mean Estimation

model$Est_area
#>         Mean         SD      2.5%       25%       50%       75%     97.5%
#> 1  1.2195835 0.08136864 1.0220924 1.1832158 1.2341439 1.2784699 1.3247949
#> 2  0.9688439 0.10752455 0.7361556 0.9078791 0.9827910 1.0476230 1.1379279
#> 3  1.2755892 0.10623761 1.0228377 1.2220985 1.2962507 1.3524717 1.4163335
#> 4  0.9691133 0.06863815 0.7900602 0.9348939 0.9816288 1.0188424 1.0578658
#> 5  0.9770920 0.16135194 0.6135998 0.8775442 0.9857523 1.0929231 1.2593776
#> 6  1.4818411 0.13679931 1.1593274 1.4106433 1.5062895 1.5822590 1.6692599
#> 7  0.8325981 0.17642277 0.4576188 0.7162542 0.8405524 0.9428438 1.1751190
#> 8  0.8398394 0.15075759 0.5359534 0.7373188 0.8467865 0.9525999 1.0907456
#> 9  0.7218695 0.18920000 0.3600545 0.5857106 0.7223206 0.8560675 1.0888686
#> 10 0.7987668 0.06318449 0.6587559 0.7668503 0.8087824 0.8448400 0.8851057

Subarea Random Effect

model$sub_randeff
#>            u_mean
#> u[1]   0.11954322
#> u[2]  -0.06004139
#> u[3]   0.18445205
#> u[4]   0.16662170
#> u[5]  -0.54732210
#> u[6]   0.33113019
#> u[7]  -0.11287492
#> u[8]   0.15000192
#> u[9]   0.38059170
#> u[10]  0.53898413
#> u[11]  0.07701193
#> u[12]  0.29220387
#> u[13] -0.21388668
#> u[14]  0.63962542
#> u[15] -0.89899382
#> u[16]  0.04158116
#> u[17]  0.92457164
#> u[18] -0.38263658
#> u[19] -1.10470647
#> u[20]  0.25160464
#> u[21]  0.03855711
#> u[22]  0.30673295
#> u[23] -0.33886331
#> u[24] -0.27597231
#> u[25] -0.74352091
#> u[26] -0.62991979
#> u[27]  0.43736251
#> u[28]  0.51104113
#> u[29] -0.45462690
#> u[30]  0.36571720

Calculate Subarea Relative Standard Error (RSE) or CV and Mean Squared Error (MSE)

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.   : 7.01   Min.   :0.004223  
#>  1st Qu.:10.20   1st Qu.:0.008060  
#>  Median :14.92   Median :0.015218  
#>  Mean   :19.59   Mean   :0.019167  
#>  3rd Qu.:27.31   3rd Qu.:0.030544  
#>  Max.   :54.76   Max.   :0.046438

Extract Coefficient Estimation \(\beta\)

model$coefficient
#>           Mean        SD       2.5%        25%       50%       75%     97.5%
#> b[0] 0.8249460 0.2502618  0.3350546 0.66517225 0.8294142 0.9815195 1.3199930
#> b[1] 0.2971333 0.3084037 -0.2788817 0.06519255 0.2916986 0.5063736 0.9023682
#> b[2] 1.2900051 0.4020163  0.4508700 1.03118281 1.2834671 1.5552284 2.0997840

Extract Area Random Effect Variance \(\sigma_u^2\) and Subarea Random Effect Variance \(\sigma_v^2\)

model$refVar
#>      b_var     a_var
#> 1 1.008098 0.9872566

STEP 7 Visualize The Result

library(ggplot2)

Save the output of Subarea estimation and the Direct Estimation (y)

df <- data.frame(
  area = seq_along(model$Est_sub$Mean),             
  direct = dataBeta$y,              
  mean_estimate = model$Est_sub$Mean
)

Area Mean Estimation

ggplot(df, aes(x = area)) +
  geom_point(aes(y = direct), size = 2, colour = "#388894", alpha = 0.6) +   # scatter points
  geom_point(aes(y = mean_estimate), size = 2, colour = "#2b707a") +   # scatter points
  geom_line(aes(y = direct), linewidth = 1, colour = "#388894", alpha = 0.6) +  # line connecting points
  geom_line(aes(y = mean_estimate), linewidth = 1, colour = "#2b707a") +  # line connecting points
  labs(
    title = "Scatter + Line Plot of Estimated Means",
    x = "Area Index",
    y = "Value"
  )

ggplot(df, aes(x = , direct, y = mean_estimate)) +
  geom_point( size = 2, colour = "#2b707a") +
   geom_abline(intercept = 0, slope = 1, color = "gray40", linetype = "dashed") +
  geom_smooth(method = "lm", color = "#2b707a", se = FALSE) +
  ylim(0, 1) +
  labs(
    title = "Scatter Plot of Direct vs Model-Based",
    x = "Direct",
    y = "Model Based"
  )

Combine the CV of direct estimation and CV from output

df_cv <- data.frame(
  direct = sqrt(dataBeta$vardir)/dataBeta$y*100,              
  cv_estimate = CV_sub
)
df_cv <- df_cv[order(df_cv$direct), ]
df_cv$area <- seq_along(dataBeta$y)

Relative Standard Error of Subarea Mean Estimation

ggplot(df_cv, aes(x = area)) +
  geom_point(aes(y = direct), size = 2, colour = "#388894", alpha = 0.5) +
  geom_point(aes(y = cv_estimate), size = 2, colour = "#2b707a") +
  ylim(0, 100) +
  labs(
    title = "Scatter Plot of Direct vs Model-Based CV",
    x = "Area",
    y = "CV (%)"
  )