robustsae {robustsae}

Robust Small Area Estimation Modeling Both Means and Variances

Description

This function provides full Bayesian Analysis for specific area-level small area models when data are provided for modeling both the mean and the variance.

Usage

robustsae(formula, S2, ni, nsim = 1000, burnin = 500, data, truemean)

Arguments

formula

a symbolic description of the model to be fitted. The details of model specification are given under Details.

S2

a vector contain the sampling variances which are given for estimating the true variances.

ni

a vector containing the sample sizes for each area.

nsim

user-specified number of MCMC draws. See German (2006).

burnin

the number of burning iterations for the sampler. See German (2006).

data

an optional data frame containing the variables named in formula, S2 and ni.

truemean

true mean values for each area.

Details

Let \theta_i denotes interest parameter for each area i, x_i the available area-specific auxiliary data, \beta the regression coefficients and m the number of small areas. A typical area level model is given by

y_i= x_i \beta + u_i + e_i, (i = 1, \ldots, m).

Assume that the random effects u_i and the sampling errors e_i are to be independently distributed with the u_i ~ N(0, \sigma^2) and the e_i ~ N(0, v_i). To foster robustness in small area estimation procedures, student t distribution is used for the random effects. Also, due to the availability of additional data purported to estimate the error variances, this considers modeling of both the means and the variances.

The robust Bayesian small area estimation model is

y_i | \theta_i ~ N(\theta_i, v_i)

S_i^2 | v_i ~ Gamma((n_i-1)/2, 1/(2v_i))

\theta_i | \beta, \sigma^2, df ~ t(x_i\beta, \sigma^2, df)

, where df is degrees of freedom parameter. For a full Bayesian analysis, this function uses the modified Jeffrey' prior which is the product of the general Jeffrey' prior and e^(-a/(2*\sigma^2)) where a is chosen as 1:

\pi(\beta) ~ 1

\pi(v_i)~ 1/v_i

\sigma^2 ~ Inv-Gamma(p/2, a/2), for a>0

\pi(df) ~ df^{-1/2} (df+1)^{p/2 -1} (df+3)^{- p/2 - 1/2}

The estimates of interest parameters are obtained by Rao-Balackwellization with Gibbs sampling with Metropolis-Hastings algorithm.

Value

The function returns a object of class "robustsae" containing the following components:

mean

Rao-Balackwellization estimates of theta's

variance

Rao-Balackwellization estimates of v's

Criteria

a list containing the following comparison criteria : Returns NA if truemean is not provided.

• ASD: average squared deviation, defined as 1/m \sum_{i=1}^{m} (\hat{\theta}_i - \theta_i)^2

• AAB: average absolute bias, defined as 1/m \sum_{i=1}^{m} |\hat{\theta}_i - \theta_i|

• ASRB: average squared relative bias, defined as 1/m \sum_{i=1}^{m} ((\hat{\theta}_i - \theta_i)/{\theta_i})^2

• ARB: average relative bias, defined as 1/m \sum_{i=1}^{m} |(\hat{\theta}_i - \theta_i)/\theta_i|

Author(s)

Malay Ghosh, Jiyoun Myung, Fernando Moura

References

Rao, J. N. K. (2003) Small Area Estimation. John Wiley and Sons.

Chip, S., and Green berg, E. (1995). Understanding the Metropolis-Hastings Algorithm. The American Statistician, 49, 327-335.

Examples

# If there is truemean data,
# load data set
data(BZdata)
attach(BZdata)

result <- robustsae(y ~ X1 + X2, S2, ni = BZdata$ni, nsim = 1000, burnin = 500, 
                      data = BZdata, truemean = truemean)
result

detach(BZdata)

# If there is no truemean data,
#load data set
data(corndata)
attach(corndata)

result2 <- robustsae(Xi ~ Z1i, Si^2, ni=corndata$ni, data = corndata) # no truemean
result2$mean
result2$variance

detach(corndata)