An Aplication to SAE HB ME under Beta Distribution On Sample Data

In saeHB.ME.beta: SAE with Measurement Error using HB under Beta Distribution

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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)

Ekstract Mean Estimation

Small Area Mean Estimation

example$Est

Coefficient Estimation

example$coefficient

Random Effect Variance Estimation

example$refvar

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)

Comparing Mean Squared Error (MSE)

MSE=as.data.frame(cbind(MSE_direct,MSE_HBMEbeta))
summary(MSE)

Comparing Relative Standard Error (RSE)

RSE=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta))
summary(RSE)

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)

Ekstract Mean Estimation

Small Area Mean Estimation

example_mix$Est

Coefficient Estimation

example_mix$coefficient

Random Effect Variance Estimation

example_mix$refvar

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)

Comparing Mean Squared Error (MSE)

MSE_mix=as.data.frame(cbind(MSE_direct,MSE_HBMEbeta_mix))
summary(MSE_mix)

Comparing Relative Standard Error (RSE)

RSE_mix=as.data.frame(cbind(RSE_direct,RSE_HBMEbeta_mix))
summary(RSE_mix)

note : you can use dataHBMEbetaNS for using dataset with non-sampled area

saeHB.ME.beta documentation built on July 9, 2023, 5:41 p.m.