# \donttest{
# Data Generation
set.seed(123)
m <- 30
x1 <- runif(m, 0, 1)
x2 <- runif(m, 1, 5)
x3 <- runif(m, 10, 15)
x4 <- runif(m, 10, 20)
b0 <- b1 <- b2 <- b3 <- b4 <- 0.5
u <- rnorm(m, 0, 1)
vardir <- 1 / rgamma(m, 1, 1)
Mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 + u
y <- rnorm(m, Mu, sqrt(vardir))
dataNormal <- as.data.frame(cbind(y, x1, x2, x3, x4, vardir))
dataNormalNs <- dataNormal
dataNormalNs$y[c(3, 10, 15, 29, 30)] <- NA
dataNormalNs$vardir[c(3, 10, 15, 29, 30)] <- NA
## Compute Fitted Model
## y ~ x1 +x2 +x3 +x4
## For data without any nonsampled area
formula <- y~x1 + x2 + x3 + x4
var <- "vardir"
v <- c(1, 1, 1, 1, 1)
c <- c(0, 0, 0, 0, 0)
## Using parameter coef and var.coef
saeHBnormal <- Normal(formula, vardir = var, coef = c, var.coef = v, data = dataNormal)
saeHBnormal$Est # Small Area mean Estimates
saeHBnormal$refVar # Random effect variance
saeHBnormal$coefficient # coefficient
# Load Library 'coda' to execute the plot
# autocorr.plot(saeHBnormal$plot[[3]]) is used to generate ACF Plot
# plot(saeHBnormal$plot[[3]]) is used to generate Density and trace plot
## Do not using parameter coef and var.coef
saeHBnormal <- Normal(formula, vardir = "vardir", data = dataNormal)
## For data with nonsampled area use dataNormalNs
# }
Run the code above in your browser using DataLab