Nothing
#' @title EBLUPs based on a Univariate Fay Herriot model with Additive Logistic Transformation for Non-Sampled Data
#' @description This function gives the transformed EBLUP based on a univariate Fay-Herriot model. Random effects for sampled domains are from the fitted model and random effects for non-sampled domains are from cluster information.
#' @param formula an object of class \code{\link[stats]{formula}} that describe the fitted model.
#' @param vardir vector containing the sampling variances of direct estimators for each domain. The values must be sorted as the variables in \code{formula}.
#' @param MAXITER maximum number of iterations allowed in the Fisher-scoring algorithm, Default: \code{100}.
#' @param PRECISION convergence tolerance limit for the Fisher-scoring algorithm, Default: \code{1e-4}.
#' @param cluster Default: \code{"auto"}. If \code{cluster = "auto"}, then the clustering will be performed by the function by finding optimal number of cluster. If cluster is a number, then clustering will be performed based on the chosen number of cluster. If cluster is a vector containing cluster information, then the vector will be used directly to find average of random effects. Clustering is performed with k-medoids algorithms using the function \code{\link[fpc]{pamk}}. If \code{"auto"} is chosen, \code{krange} are set to \code{2:(nrow(data)-1)}.
#' @param data optional data frame containing the variables named in \code{formula} and \code{vardir}.
#' @return The function returns a list with the following objects:
#' \item{est}{a data frame containing values of the estimators for each domains.}
#' \itemize{
#' \item \code{PC} : transformed EBLUP estimators using inverse alr.
#' \item \code{status} : status of corresponding domain, whether sampled or non-sampled.
#' \item \code{cluster} : cluster of corresponding domain.
#' }
#' \item{fit}{a list containing the following objects (model is fitted using REML):}
#' \itemize{
#' \item \code{convergence} : a logical value equal to \code{TRUE} if Fisher-scoring algorithm converges in less than \code{MAXITER} iterations.
#' \item \code{iterations} : number of iterations performed by the Fisher-scoring algorithm.
#' \item \code{estcoef} : a data frame that contains the estimated model coefficients, standard errors, t-statistics, and p-values of each coefficient.
#' \item \code{refvar} : estimated random effects variance.
#' \item \code{cluster.information} : a data frame containing average random effects of sampled domain in each cluster.
#' }
#' \item{components}{a data frame containing the following columns:}
#' \itemize{
#' \item \code{random.effects} : estimated random effect values of the fitted model.
#' \item \code{residuals} : residuals of the fitted model.
#' \item \code{status} : status of corresponding domain, whether sampled or non-sampled.
#' }
#'
#' @examples
#' \dontrun{
#' ## Load dataset
#' data(datasaeu.ns)
#'
#' ## If data is defined
#' Fo = y ~ x1 + x2
#' vardir = "vardir"
#' model.ns <- saeFH.ns.uprop(Fo, vardir, data = datasaeu.ns)
#'
#' ## If data is undefined (and option for cluster arguments)
#' Fo = datasaeu.ns$y ~ datasaeu.ns$x1 + datasaeu.ns$x2
#' vardir = datasaeu.ns$vardir
#'
#' ### "auto"
#' model.ns1 <- saeFH.ns.uprop(Fo, vardir, cluster = "auto")
#'
#' ### number of clusters
#' model.ns2 <- saeFH.ns.uprop(Fo, vardir, cluster = 2)
#'
#' ### vector containing cluster for each domain
#' model.ns3 <- saeFH.ns.uprop(Fo, vardir, cluster = datasaeu.ns$cluster)
#'
#' ## See the estimators
#' model.ns$est
#' }
#'
#' @export saeFH.ns.uprop
# SAE Univariate for Non-Sampled Area Function
saeFH.ns.uprop = function(formula, vardir,
MAXITER = 100,
PRECISION = 1e-4,
cluster = "auto",
data) {
# require(fpc)
# Setting List for Results
result = list(est = NA,
fit = list(convergence = TRUE,
iterations = 0,
estcoef = NA,
refvar = NA,
cluster.information = NA),
components = data.frame(random.effects = NA,
residuals = NA)
)
# Getting Data
if (!missing(data)) {
formuladata = model.frame(formula, na.action = na.pass, data)
X = model.matrix(formula, formuladata)
} else{
formuladata = model.frame(formula, na.action = na.pass)
X = model.matrix(formula, formuladata)
}
Z = formuladata[,1]
if (any(na.omit(Z) < 0 | na.omit(Z) > 1)) {
stop("Proportion in a domain must fall between 0 and 1")
}
D = length(Z)
non.sampled = which(Z == 0 | Z == 1 | is.na(Z))
# Getting Vardir
namevar = deparse(substitute(vardir))
if (is.numeric(vardir)) {
vardir = vardir
} else if(is.character(vardir)) {
if (missing(data)) {
stop("If vardir is character, data need to be defined")
} else {
vardir = data[, vardir]
}
}
if (length(non.sampled) > 0) {
if (any(is.na(vardir[-non.sampled]))) {
stop("If value of a domain is not [0, 1, or NA], vardir for corresponding domain must be defined")
}
} else {
if (any(is.na(vardir))) {
stop("All domain are sampled, all vardir must be defined")
}
}
# Data Transformation (alr)
y = log(Z / (1 - Z))
if (length(non.sampled) > 0) {
y.sm = y[-non.sampled]
} else {
y.sm = y
}
# Vardir Transformation
q = 2
H0 = q * (diag(1, q - 1) + matrix(1, nrow = q - 1) %*% t(matrix(1, nrow = q - 1)))
vardir = as.numeric(H0^2) * vardir
if (length(non.sampled) > 0) {
vardir.sm = vardir[-non.sampled]
} else {
vardir.sm = vardir
}
# Cluster information
if (length(non.sampled) > 0) {
## Clustering
if (length(cluster) == 1) {
if (cluster == "auto") {
klas = pamk(X[,-1], scaling = T, krange = 2:(D - 1))
} else if(is.numeric(cluster) & (cluster > 1)) {
klas = pamk(X[,-1], scaling = T, krange = cluster)
} else {
stop("Invalid choice of cluster clusters")
}
clust.df = klas$pamobject$clustering
} else {
if (length(cluster) != nrow(X)) {
stop("Cluster information length is not appropriate with the data")
}
clust.df = cluster
}
clust.data = model.matrix(~., data.frame(class = as.factor(clust.df)))
if (any(apply(clust.data[-non.sampled,], 2, function(x){all(x == 0)}))) {
stop("A cluster may not contain all non-sampled area, please select other number of cluster or give other cluster information")
}
## New X Matrix
nameX = c(colnames(X), colnames(clust.data)[-1])
X = cbind(X, clust.data[,-1])
colnames(X) = nameX
}
if (length(non.sampled) > 0) {
X.sm = X[-non.sampled,]
X.ns = X[non.sampled,]
} else {
X.sm = X
}
# Estimating Variance
## Fisher-scoring algorithm for REML estimator for variance
### Initial value of variance using median of sampling variance vardir
Vu.est = 0
Vu.est[1] = median(vardir.sm)
iter = 0
diff = PRECISION + 1
while ((diff > PRECISION) & (iter < MAXITER)) {
iter = iter + 1
V.Inv = 1 / (Vu.est[iter] + vardir.sm)
XtV.Inv = t(V.Inv * X.sm)
Q = solve(XtV.Inv %*% X.sm)
P = diag(V.Inv) - t(XtV.Inv) %*% Q %*% XtV.Inv
Py = P %*% y.sm
### Score function
s = -0.5 * sum(diag(P)) + 0.5 * (t(Py) %*% Py)
### Fisher information
Fi = 0.5 * sum(diag(P %*% P))
### Updating equation
Vu.est[iter + 1] = Vu.est[iter] + s/Fi
### Relative difference
diff = abs((Vu.est[iter + 1] - Vu.est[iter]) / Vu.est[iter])
}
result$fit$iterations = iter
if ((iter >= MAXITER)) {
result$fit$convergence = FALSE
return(result)
}
# Final Estimation of Variance of Random Effects
Vu = max(Vu.est[iter + 1], 0)
# Coefficient Estimator Beta
V.Inv = 1 / (Vu + vardir.sm)
XtV.Inv = t(V.Inv * X.sm)
Q = solve(XtV.Inv %*% X.sm)
beta.REML = Q %*% XtV.Inv %*% y.sm
# Std error & p-value
std.error = sqrt(diag(Q))
t.value = beta.REML / std.error
p.value = 2 * pnorm(abs(t.value), lower.tail = F)
Xbeta.REML = X.sm %*% beta.REML
resid = y.sm - Xbeta.REML
# Random Effects & EBLUP Predictor
u.hat = Vu * V.Inv * resid
EBLUP = Xbeta.REML + u.hat
# Cluster information in action
if (length(non.sampled) > 0) {
## Mean of Random Effects in each cluster
u.mean = setNames(aggregate(x = u.hat,
by = list(clust.df[-non.sampled]),
FUN = mean), c("cluster", "mean.random.effect"))
result$fit$cluster.information = u.mean
## Using Random Effects Means to Non-sampled Area
EBLUP.sm = EBLUP
u.hat.sm = u.hat
u.hat.ns = apply(matrix(clust.df[non.sampled]), 1, function(x){
u.mean[u.mean$cluster == x, 2]
})
EBLUP.ns = X.ns %*% beta.REML + u.hat.ns
EBLUP = matrix(nrow = D, ncol = 1)
EBLUP[-non.sampled] = EBLUP.sm
EBLUP[non.sampled] = EBLUP.ns
u.hat = matrix(nrow = D, ncol = 1)
u.hat[-non.sampled] = u.hat.sm
u.hat[non.sampled] = u.hat.ns
Xbeta.REML = X %*% beta.REML
## Missing values in y will be replaced by EBLUP estimator
y[non.sampled] = EBLUP.ns
## Mean of Vardir in each cluster
var.mean = setNames(aggregate(x = vardir.sm,
by = list(clust.df[-non.sampled]),
FUN = mean), c("cluster", "mean.var"))
var.ns = apply(matrix(clust.df[non.sampled]), 1, function(x){
var.mean[var.mean$cluster == x, 2]
})
vardir[non.sampled] = var.ns
}
# Compositional Plug-in Predictors
## Transformation to Proportion (alr)
PC = exp(EBLUP) / (1 + exp(EBLUP))
status = matrix("Sampled", nrow = D)
status[non.sampled,] = "Non-Sampled"
# Results
result$est = data.frame(PC = PC, status = status)
if (length(non.sampled) > 0) {
result$est = data.frame(result$est,
cluster = clust.df)
}
result$fit$estcoef = data.frame(beta = beta.REML,
std.error = std.error,
t.value = t.value,
p.value = p.value)
result$fit$refvar = Vu
result$components = data.frame(random.effects = u.hat,
residuals = (y - EBLUP),
status = status)
result$components$residuals[non.sampled] = NA
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.