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.2198615 0.07921056 1.0210870 1.1815151 1.2384559 1.2773143 1.320807
#> 2  0.9695509 0.10830161 0.7167013 0.9101108 0.9815290 1.0460518 1.145187
#> 3  1.2737726 0.11356345 0.9563260 1.2179782 1.3006111 1.3527973 1.422594
#> 4  0.9689708 0.07263555 0.7881512 0.9373074 0.9822225 1.0201324 1.062958
#> 5  0.9831896 0.16376557 0.6143263 0.8811928 1.0045183 1.1010920 1.261260
#> 6  1.4824635 0.12682200 1.1811936 1.4089117 1.5023257 1.5768754 1.663475
#> 7  0.8193407 0.17565815 0.4702812 0.6999436 0.8209256 0.9430062 1.152576
#> 8  0.8538585 0.15597183 0.5332365 0.7538820 0.8680240 0.9643364 1.104871
#> 9  0.7216715 0.17151776 0.3798080 0.6097727 0.7160974 0.8420327 1.043955
#> 10 0.7996785 0.05956151 0.6568060 0.7693301 0.8098963 0.8446335 0.881005

Area Random Effect

model$area_randeff
#>            f_mean
#> f[1]   0.26392757
#> f[2]   0.01207837
#> f[3]   0.44271438
#> f[4]   0.72800668
#> f[5]  -0.39604354
#> f[6]   0.53313490
#> f[7]  -0.77557184
#> f[8]  -0.34589159
#> f[9]  -0.91406595
#> f[10]  0.35901381

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.493   Min.   :0.003548  
#>  1st Qu.: 7.761   1st Qu.:0.007638  
#>  Median :10.043   Median :0.014490  
#>  Mean   :13.021   Mean   :0.016723  
#>  3rd Qu.:17.864   3rd Qu.:0.026196  
#>  Max.   :23.767   Max.   :0.030856

STEP 6 Extract Result for Subarea

Subarea Mean Estimation

model$Est_area
#>         Mean         SD      2.5%       25%       50%       75%    97.5%
#> 1  1.2198615 0.07921056 1.0210870 1.1815151 1.2384559 1.2773143 1.320807
#> 2  0.9695509 0.10830161 0.7167013 0.9101108 0.9815290 1.0460518 1.145187
#> 3  1.2737726 0.11356345 0.9563260 1.2179782 1.3006111 1.3527973 1.422594
#> 4  0.9689708 0.07263555 0.7881512 0.9373074 0.9822225 1.0201324 1.062958
#> 5  0.9831896 0.16376557 0.6143263 0.8811928 1.0045183 1.1010920 1.261260
#> 6  1.4824635 0.12682200 1.1811936 1.4089117 1.5023257 1.5768754 1.663475
#> 7  0.8193407 0.17565815 0.4702812 0.6999436 0.8209256 0.9430062 1.152576
#> 8  0.8538585 0.15597183 0.5332365 0.7538820 0.8680240 0.9643364 1.104871
#> 9  0.7216715 0.17151776 0.3798080 0.6097727 0.7160974 0.8420327 1.043955
#> 10 0.7996785 0.05956151 0.6568060 0.7693301 0.8098963 0.8446335 0.881005

Subarea Random Effect

model$sub_randeff
#>            u_mean
#> u[1]   0.11886871
#> u[2]  -0.10718848
#> u[3]   0.25015214
#> u[4]   0.14210581
#> u[5]  -0.51331581
#> u[6]   0.30298135
#> u[7]  -0.10243186
#> u[8]   0.19048195
#> u[9]   0.39559507
#> u[10]  0.61554565
#> u[11]  0.09843486
#> u[12]  0.18308245
#> u[13] -0.15302587
#> u[14]  0.68070419
#> u[15] -0.92907007
#> u[16]  0.04992264
#> u[17]  0.93627138
#> u[18] -0.40715528
#> u[19] -1.02010304
#> u[20]  0.22223790
#> u[21]  0.05457168
#> u[22]  0.32847023
#> u[23] -0.40355500
#> u[24] -0.26265698
#> u[25] -0.79676652
#> u[26] -0.69647058
#> u[27]  0.39545778
#> u[28]  0.44096537
#> u[29] -0.50873886
#> u[30]  0.42589613

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.115   Min.   :0.004315  
#>  1st Qu.:10.697   1st Qu.:0.009007  
#>  Median :14.579   Median :0.014487  
#>  Mean   :19.383   Mean   :0.018670  
#>  3rd Qu.:25.879   3rd Qu.:0.028655  
#>  Max.   :52.988   Max.   :0.043302

Extract Coefficient Estimation \(\beta\)

model$coefficient
#>           Mean        SD       2.5%        25%       50%      75%     97.5%
#> b[0] 0.8213243 0.2642550  0.3243395 0.64883516 0.8135892 1.003231 1.3311506
#> b[1] 0.2938539 0.3303445 -0.3786486 0.07039325 0.2931766 0.508798 0.9429934
#> b[2] 1.3318983 0.3942418  0.5515258 1.06956619 1.3214754 1.596781 2.1284576

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

model$refVar
#>       b_var    a_var
#> 1 0.9453136 1.000931

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 (%)"
  )