# \donttest{
## Data Generation
set.seed(123)
m <- 30
x1 <- runif(m, 0, 1)
x2 <- runif(m, 0, 1)
b0 <- b1 <- b2 <- 0.5
u <- rnorm(m, 0, 1)
phi <- rgamma(m, 0.5, 0.5)
vardir <- 1 / phi
mu <- exp(b0 + b1 * x1 + b2 * x2 + u)
A <- mu^2 * phi
B <- mu * phi
y <- rgamma(m, A, B)
dataGamma <- as.data.frame(cbind(y, x1, x2, vardir))
dataGammaNs <- dataGamma
dataGammaNs$y[c(3, 14, 22, 29, 30)] <- NA
dataGammaNs$vardir[c(3, 14, 22, 29, 30)] <- NA
## Compute Fitted Model
## y ~ x1 +x2
## For data without any nonsampled area
model_formula = y ~ x1 +x2
v <- c(1, 1, 1)
c <- c(0, 0, 0)
## Using parameter coef and var.coef
saeHBGamma <- Gamma(model_formula, coef = c, var.coef = v, iter.update = 10, data = dataGamma)
saeHBGamma$Est # Small Area mean Estimates
saeHBGamma$refVar # Random effect variance
saeHBGamma$coefficient # coefficient
# Load Library 'coda' to execute the plot
# autocorr.plot(saeHBGamma$plot[[3]]) is used to generate ACF Plot
# plot(saeHBGamma$plot[[3]]) is used to generate Density and trace plot
## Do not using parameter coef and var.coef
saeHBGamma <- Gamma(model_formula, data = dataGamma) #'
## For data with nonsampled area use dataGammaNs
# }
Run the code above in your browser using DataLab