# \donttest{
# Data Generation
set.seed(123)
m <- 30
x1 <- runif(m, 0, 1)
x2 <- runif(m, 0, 1)
x3 <- runif(m, 0, 1)
x4 <- runif(m, 0, 1)
b0 <- b1 <- b2 <- b3 <- b4 <- 0.5
u <- rnorm(m, 0, 1)
pi <- rgamma(1, 1, 0.5)
m <- exp(b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 + u)
Mu <- m / (1 + m)
A <- Mu * pi
B <- (1 - Mu) * pi
y <- rbeta(m, A, B)
y <- ifelse(y < 1, y, 0.99999999)
y <- ifelse(y > 0, y, 0.00000001)
MU <- A / (A + B)
vardir <- A * B / ((A + B)^2 * (A + B + 1))
dataBeta <- as.data.frame(cbind(y, x1, x2, x3, x4, vardir))
dataBetaNs <- dataBeta
dataBetaNs$y[c(3, 14, 22, 29, 30)] <- NA
dataBetaNs$vardir[c(3, 14, 22, 29, 30)] <- NA
## Compute Fitted Model
## y ~ x1 +x2
## For data without any nonsampled area
vc <- c(1, 1, 1)
c <- c(0, 0, 0)
formula <- y~x1 + x2
dat <- dataBeta[1:10, ]
## Using parameter coef and var.coef
saeHBbeta <- Beta(formula, var.coef = vc, iter.update = 10, coef = c, data = dat)
saeHBbeta$Est # Small Area mean Estimates
saeHBbeta$refVar # Random effect variance
saeHBbeta$coefficient # coefficient
# Load Library 'coda' to execute the plot
# autocorr.plot(saeHBbeta$plot[[3]]) # is used to generate ACF Plot
# plot(saeHBbeta$plot[[3]]) # is used to generate Density and trace plot
## Do not use parameter coef and var.coef
saeHBbeta <- Beta(formula, data = dataBeta)
# }
Run the code above in your browser using DataLab