# R/mse_msaeRBns.R In msaeRB: Ratio Benchmarking for Multivariate Small Area Estimation

#' @title Parametric Bootstrap Mean Squared Error Estimators of Ratio Benchmarking for Multivariate Non Sampled Area in Small Area Estimation
#'
#' @description Calculates the parametric bootstrap mean squared error estimates of ratio benchmarking for multivariate non sampled area in small area estimation
#'
#' @param formula an object of class list of formula describe the fitted models
#' @param vardir matrix containing sampling variances of direct estimators. The order is: \code{var1, cov12, ..., cov1r, var2, cov23, ..., cov2r, ..., cov(r-1)(r), var(r)}
#' @param weight matrix containing proportion of units in small areas. The order is: \code{w1, w2, ..., w(r)}
#' @param cluster matrix containing cluster of auxiliary variables. The order is: \code{c1, c2, ..., c(r)}
#' @param samevar logical. If \code{TRUE}, the varians is same. Default is \code{FALSE}
#' @param B number of bootstrap. Default is 1000
#' @param MAXITER maximum number of iterations for Fisher-scoring. Default is 100
#' @param PRECISION coverage tolerance limit for the Fisher Scoring algorithm. Default value is \code{1e-4}
#' @param data dataframe containing the variables named in formula, vardir, and weight
#'
#' @return
#' \item{mse.eblup}{estimated mean squared errors of the EBLUPs for the small domains based on Prasad Rao}
#' \item{pbmse.eblupRB}{parametric bootstrap mean squared error estimates of the ratio benchmark}
#' \item{running.time}{time for running function}
#'
#' @export mse_msaeRBns
#'
#' @import abind
#' @importFrom magic adiag
#' @importFrom Matrix forceSymmetric
#' @importFrom stats model.frame na.omit model.matrix median pnorm rnorm
#' @importFrom MASS mvrnorm
#'
#' @examples
#' \donttest{
#' ## load dataset
#' data(datamsaeRBns)
#'
#' # Compute MSE EBLUP and Ratio Benchmark
#' # This is the long running example
#' ## Using parameter 'data'
#' Fo = list(f1 = Y1 ~ X1 + X2,
#'           f2 = Y2 ~ X1 + X2,
#'           f3 = Y3 ~ X1 + X2)
#' vardir = c("v1", "v12", "v13", "v2", "v23", "v3")
#' weight = c("w1", "w2", "w3")
#' cluster = c("c1", "c2", "c3")
#'
#' mse_msae = mse_msaeRBns(Fo, vardir, weight, cluster, data = datamsaeRBns)
#'
#' ## Without parameter 'data'
#' Fo = list(f1 = datamsaeRBns$Y1 ~ datamsaeRBns$X1 + datamsaeRBns$X2, #' f2 = datamsaeRBns$Y2 ~ datamsaeRBns$X1 + datamsaeRBns$X2,
#'           f3 = datamsaeRBns$Y3 ~ datamsaeRBns$X1 + datamsaeRBns$X2) #' vardir = datamsaeRBns[, c("v1", "v12", "v13", "v2", "v23", "v3")] #' weight = datamsaeRBns[, c("w1", "w2", "w3")] #' cluster = datamsaeRBns[, c("c1", "c2", "c3")] #' #' mse_msae = mse_msaeRBns(Fo, vardir, weight, cluster) #' #' ## Return #' mse_msae$pbmse.eblupRB # to see the MSE of Ratio Benchmark
#' }
mse_msaeRBns = function (formula, vardir, weight, cluster, samevar = FALSE, B = 1000, MAXITER = 100, PRECISION = 1e-04, data) {
start_time = Sys.time()
r = length(formula)
if (r <= 1)
stop("You should use mse_saeRBns() for univariate")
R_function = function(vardir, n, r) {
if (r == 1) {
R = diag(vardir)
}
else {
R = matrix(rep(0, times = n * r * n * r), nrow = n *
r, ncol = n * r)
k = 1
for (i in 1:r) {
for (j in 1:r) {
if (i <= j) {
mat0 = matrix(rep(0, times = r * r), nrow = r,
ncol = r)
mat0[i, j] = 1
matVr = diag(vardir[, k], length(vardir[,
k]))
R_hasil = kronecker(mat0, matVr)
R = R + R_hasil
k = k + 1
}
}
}
R = forceSymmetric(R)
}
return(as.matrix(R))
}
temporary_vardir = function(vardir, n, r) {
dt = list()
for (h in 1:n) {
var_mat = matrix(0, nrow = r, ncol = r)
k = 1
for (i in 1:r) {
for (j in 1:r) {
if (i <= j) {
var_mat[i, j] = vardir[h, k]
k = k + 1
}
}
}
var_mat = forceSymmetric(var_mat)
dt[[h]] = as.matrix(var_mat)
}
return(dt)
}
eblup_inside = function(r, n, samevar, y, X, cluster, vardir, MAXITER = 100, PRECISION = 1e-04){
y.matrix = matrix(as.vector(y), n, r)
y.matrix[y.matrix==0] = NA
indexns = unique(which(is.na(y.matrix), arr.ind=TRUE)[,1])
indexs = c(1:n)[-indexns]
n.s = length(indexs)
n.ns = n - n.s
y.s = y.matrix[-indexns, ]
y.s.vec = as.vector(y.s)
corr = cor(y.s)
hasil_cor = c()
for (i in 1:r){
for (j in 1:r){
if (i < j && i != j){
hasil_cor = c(hasil_cor, corr[i, j])
}
}
}
vardir[indexns, ] = NA
index.vardir = sum(1:r)
for (i in 2:r){
index.vardir[i] = index.vardir[1] - sum(2:i)
}
index.vardir = sort(index.vardir)
index.covardir = c(1:sum(1:r))[-index.vardir]
varians.direct = vardir[, index.vardir]
covarians.direct = vardir[, -index.vardir]
for (i in 1:r) {
df.vardir = data.frame(cluster[, i], varians.direct[, i])
for (j in 1:n) {
if (df.vardir[j, 2] %in% NA) {
df.vardir[j, 2] = mean(df.vardir[df.vardir[, 1] == df.vardir[j,1], 2], na.rm = TRUE)
}
}
varians.direct[ ,i] = df.vardir[, 2]
}
index.cor = 0
if (r > 2){
for (p in 1:r) {
for (q in 1:r) {
if (p < q) {
index.cor = index.cor + 1
for (col in 1:n) {
if (covarians.direct[col, index.cor] %in% NA) {
covarians.direct[col, index.cor] = hasil_cor[index.cor] * sqrt(varians.direct[col, p] * varians.direct[col, q])
}
}
}
}
}
} else {
index.cor = index.cor + 1
for (col in 1:n) {
if (covarians.direct[col] %in% NA) {
covarians.direct[col] = hasil_cor[index.cor] * sqrt(varians.direct[col, 1] * varians.direct[col, 2])
}
}
}
varians.direct = rbind(index.vardir, varians.direct)
if (r > 2) {
covarians.direct = rbind(index.covardir, covarians.direct)
} else {
covarians.direct = c(index.covardir, covarians.direct)
}
vardir.full = cbind(varians.direct, covarians.direct)
vardir.full = vardir.full[, order(vardir.full[1, ])]
vardir.full = vardir.full[-1, ]
R.ns = R_function(vardir.full[indexns, ], n.ns, r)
vardir.s = vardir[-indexns, ]
R = R_function(vardir.s, n.s, r)
W.s = W[-indexns, ]
W.s = prop.table(W.s, 2)
Xindexns = c()
for (i in 1:r) {
Xindexns = c(Xindexns, indexns + rep(n, times = n.ns)*(i - 1))
}
X.s = X[-Xindexns, ]
X.ns = X[Xindexns, ]
y_names = sapply(formula, "[[", 2)
Ir = diag(r)
In = diag(n.s)
dV = list()
dV1 = list()
for (i in 1:r){
dV[[i]] = matrix(0, nrow = r, ncol = r)
dV[[i]][i, i] = 1
dV1[[i]] = kronecker(dV[[i]], In)
}
convergence = TRUE
if (samevar) {
Vu = median(diag(R))
k = 0
diff = rep(PRECISION + 1, r)
while (any(diff > PRECISION) & (k < MAXITER)) {
k = k + 1
Vu1 = Vu
Gr = kronecker(Vu1, Ir)
Gn = kronecker(Gr, In)
V = as.matrix(Gn + R)
Vinv = solve(V)
XstVinv = t(Vinv %*% X.s)
Q = solve(XstVinv %*% X.s)
P = Vinv - t(XstVinv) %*% Q %*% XstVinv
Py = P %*% y.s.vec
s = (-0.5) %*% sum(diag(P)) + 0.5 %*% (t(Py) %*% Py)
iF =  0.5 %*% sum(diag(P %*% P))
Vu = Vu1 + solve(iF) %*% s
diff = abs((Vu - Vu1)/Vu1)
}
Vu = as.vector((rep(max(Vu, 0), r)))
names(Vu) = y_names
if (k >= MAXITER && diff >= PRECISION) {
convergence = FALSE
}
Gn = kronecker(diag(Vu), In)
Gn.ns = kronecker(diag(Vu), diag(n.ns))
V = as.matrix(Gn + R)
Vinv = solve(V)
XstVinv = t(Vinv %*% X.s)
Q = solve(XstVinv %*% X.s)
P = Vinv - t(XstVinv) %*% Q %*% XstVinv
Py = P %*% y.s.vec
beta = Q %*% XstVinv %*% y.s.vec
res = y.s.vec - X.s %*% beta
random.effect = data.frame(matrix(Gn %*% Vinv %*% res, n.s, r))
names(random.effect) = y_names
se.b = sqrt(diag(Q))
t.value = beta/se.b
p.value = 2 * pnorm(abs(as.numeric(t.value)), lower.tail = FALSE)
coef = as.matrix(cbind(beta, se.b, t.value, p.value))
colnames(coef) = c("beta", "std. error", "t value", "p-value")
rownames(coef) = colnames(X)
Gn.ns = kronecker(diag(Vu), diag(n.ns))
Vinv.ns = solve(Gn.ns + R.ns)
Q.ns = ginv(t(Vinv.ns %*% X.ns) %*% X.ns)
g1.s = diag(Gn %*% Vinv %*% R)
g2.s = diag(R %*% Vinv %*% X.s %*% Q %*% t(X.s) %*% t(R %*% Vinv))
g1.ns = diag(Gn.ns %*% Vinv.ns %*% R.ns)
g2.ns = diag(R.ns %*% Vinv.ns %*% X.ns %*% Q.ns %*% t(X.ns) %*% t(R.ns %*% Vinv.ns))
g12.s = matrix(g1.s + g2.s, n.s, r)
names(g12.s) = y_names
g12.s = data.frame(index = indexs, g12.s)
g12.ns = matrix(g1.ns + g2.ns, n.ns, r)
names(g12.ns) = y_names
g12.ns = data.frame(index = indexns, g12.ns)
g12 = rbind(g12.s, g12.ns)
g12 = g12[order(g12$index), ] rownames(g12) = g12$index
g12 = g12[, -1]
dg = Vinv - Gn %*% Vinv %*% Vinv
gg3 = (dg %*% V %*% t(dg))/iF
g3 = diag(gg3)
g3.matrix = matrix(g3, n.s, r)
g3.s = data.frame(index = indexs, g3.matrix)
g3.ns = data.frame(index = indexns, matrix(NA, n.ns, r))
names(g3.ns) = names(g3.s)
g3.all = rbind(g3.s, g3.ns)
g3.all = g3.all[order(g3.all$index), ] rownames(g3.all) = g3.all$index
g3.all = g3.all[ ,-1]
g3.full = matrix(NA, nrow = n, ncol = r)
colnames(g3.full) = y_names
for (i in 1:r) {
df = data.frame(cluster[, i], g3.all[, i])
for (j in 1:n) {
if (df[j, 2] %in% NA) {
df[j, 2] = mean(df[df[, 1] == df[j,1], 2], na.rm = TRUE)
}
}
g3.full[ ,i] = df[, 2]
}
names(g3.full) = y_names
mse.df = as.data.frame(g12 + 2 * g3.full)
} else {
Vu = apply(matrix(diag(R), nrow = n.s, ncol = r), 2, median)
k = 0
diff = rep(PRECISION + 1, r)
while (any(diff > rep(PRECISION, r)) & (k < MAXITER)) {
k = k + 1
Vu1 = Vu
if (r == 1) {
Gr = Vu1
} else {
Gr = diag(as.vector(Vu1))
}
Gn = kronecker(Gr, In)
V = as.matrix(Gn + R)
Vinv = solve(V)
XstVinv = t(Vinv %*% X.s)
Q = solve(XstVinv %*% X.s)
P = Vinv - t(XstVinv) %*% Q %*% XstVinv
Py = P %*% y.s.vec
s = sapply(dV1, function(x) (-0.5) * sum(diag(P %*% x)) + 0.5 * (t(Py) %*% x %*% Py))
iF = matrix(unlist(lapply(dV1, function(x) lapply(dV1, function(y) 0.5 * sum(diag(P %*% x %*% P %*% y))))), r)
Vu = Vu1 + solve(iF) %*% s
diff = abs((Vu - Vu1)/Vu1)
}
Vu = as.vector(sapply(Vu, max, 0))
if (k >= MAXITER && diff >= PRECISION){
convergence = FALSE
}
if (r == 1) {
Gr = Vu
} else {
Gr = diag(as.vector(Vu))
}
Gn = kronecker(Gr, In)
V = as.matrix(Gn + R)
Vinv = solve(V)
XtVinv = t(Vinv %*% X.s)
Q = solve(XtVinv %*% X.s)
P = Vinv - t(XtVinv) %*% Q %*% XtVinv
Py = P %*% y.s.vec
beta = Q %*% XtVinv %*% y.s.vec
res = y.s.vec - X.s %*% beta
random.effect = data.frame(matrix(Gn %*% Vinv %*% res, n.s, r))
names(random.effect) = y_names
se.b = sqrt(diag(Q))
t.value = beta/se.b
p.value = 2 * pnorm(abs(as.numeric(t.value)), lower.tail = FALSE)
coef = as.matrix(cbind(beta, se.b, t.value, p.value))
colnames(coef) = c("beta", "std. error", "t value", "p-value")
rownames(coef) = colnames(X)
FI = solve(iF)
Gn.ns = kronecker(Gr, diag(n.ns))
Vinv.ns = solve(Gn.ns + R.ns)
XtVinv.ns = t(Vinv.ns %*% X.ns)
Q.ns = ginv(XtVinv.ns %*% X.ns)
g1.s = diag(Gn %*% Vinv %*% R)
g2.s = diag(R %*% Vinv %*% X.s %*% Q %*% t(X.s) %*% t(R %*% Vinv))
g1.ns = diag(Gn.ns %*% Vinv.ns %*% R.ns)
g2.ns = diag(R.ns %*% Vinv.ns %*% X.ns %*% Q.ns %*% t(X.ns) %*% t(R.ns %*% Vinv.ns))
g12.s = matrix(g1.s + g2.s, n.s, r)
names(g12.s) = y_names
g12.s = data.frame(index = indexs, g12.s)
g12.ns = matrix(g1.ns + g2.ns, n.ns, r)
names(g12.ns) = y_names
g12.ns = data.frame(index = indexns, g12.ns)
g12 = rbind(g12.s, g12.ns)
g12 = g12[order(g12$index), ] rownames(g12) = g12$index
g12 = g12[, -1]
dg = lapply(dV1, function(x) x %*% Vinv - Gn %*% Vinv %*% x %*% Vinv)
gg3 = list()
for (i in 1:r) {
for (j in 1:r) {
gg3[[(i - 1) * r + j]] = FI[i, j] * (dg[[i]] %*% V %*% t(dg[[j]]))
}
}
g3 = diag(Reduce("+", gg3))
g3.matrix = matrix(g3, n.s, r)
g3.s = data.frame(index = indexs, g3.matrix)
g3.ns = data.frame(index = indexns, matrix(NA, n.ns, r))
names(g3.ns) = names(g3.s)
g3.all = rbind(g3.s, g3.ns)
g3.all = g3.all[order(g3.all$index), ] rownames(g3.all) = g3.all$index
g3.all = g3.all[ ,-1]
g3.full = matrix(NA, nrow = n, ncol = r)
colnames(g3.full) = y_names
for (i in 1:r) {
df = data.frame(cluster[, i], g3.all[, i])
for (j in 1:n) {
if (df[j, 2] %in% NA) {
df[j, 2] = mean(df[df[, 1] == df[j,1], 2], na.rm = TRUE)
}
}
g3.full[ ,i] = df[, 2]
}
names(g3.full) = y_names
mse.df = as.data.frame(g12 + 2 * g3.full)
}
random.effect = data.frame(index = indexs, random.effect)
random.effect.ns = data.frame(index = indexns, matrix(NA, n.ns, r))
names(random.effect.ns) = names(random.effect)
random.effect = rbind(random.effect, random.effect.ns)
random.effect = random.effect[order(random.effect$index), ] rownames(random.effect) = random.effect$index
random.effect = random.effect[ ,-1]
random.effect.full = matrix(NA, nrow = n, ncol = r)
colnames(random.effect.full) = y_names
for (i in 1:r) {
df = data.frame(cluster[, i], random.effect[, i])
for (j in 1:n) {
if (df[j, 2] %in% NA) {
df[j, 2] = mean(df[df[, 1] == df[j,1], 2], na.rm = TRUE)
}
}
random.effect.full[ ,i] = df[, 2]
}
Xbeta = X %*% beta
eblup = matrix(Xbeta, nrow = n, ncol = r) + random.effect.full
eblup = as.data.frame(eblup)
names(eblup) = y_names
random.effect.full = as.data.frame(random.effect.full)
eblup.mat = as.matrix(eblup)
eblup.ratio = matrix(0, nrow = n, ncol = r)
alfa = c()
for (i in 1:r) {
eblup.ratio[, i] = eblup.mat[, i] * (sum(W.s[, i] * y.s[, i])/sum(W[, i] * eblup.mat[, i]))
}
eblup.ratio = as.data.frame(eblup.ratio)
names(eblup.ratio) = y_names
agregation.direct = diag(t(W.s) %*% y.s)
agregation.eblup = diag(t(W) %*% eblup.mat)
agregation.eblup.ratio = diag(t(W) %*% as.matrix(eblup.ratio))
agregation = as.matrix(rbind(agregation.direct, agregation.eblup, agregation.eblup.ratio))
colnames(agregation) = y_names
result = list(eblup = list(est.eblup = NA, est.eblupRB = NA), fit = list(method = NA, convergence = NA, iteration = NA, estcoef = NA, refvar = NA), random.effect = NA, agregation = NA, mse = NA, g12 = NA, R.s = NA, R.ns = NA)
result$eblup$est.eblup = eblup
result$eblup$est.eblupRB = eblup.ratio
result$fit$method = "REML"
result$fit$convergence = convergence
result$fit$iteration = k
result$fit$estcoef = coef
result$fit$refvar = t(Vu)
result$random.effect = random.effect.full result$mse = mse.df
result$g12 = g12 result$R.s = R
result$R.ns = R.ns return(result) } if (!missing(data)) { formuladata = lapply(formula, function(x) model.frame(x, na.action = na.omit, data)) y = unlist(lapply(formula, function(x) model.frame(x, na.action = na.omit, data)[[1]])) X = Reduce(adiag, lapply(formula, function(x) model.matrix(x, data))) W = as.matrix(data[, weight]) n = length(y)/r if (!all(vardir %in% names(data))) stop("Object vardir is not appropriate with data.") if (length(vardir) != sum(1:r)) stop("Length of vardir is not appropriate with data. The length must be ", sum(1:r)) if (any(is.na(data[, weight]))) stop("Object weight contains NA values.") if (!all(weight %in% names(data))) stop("Object weight is not appropriate with data.") if (length(weight) != r) stop("Length of weight is not appropriate with data. The length must be ", r) if (!all(cluster %in% names(data))) stop("Object cluster is not appropriate with data.") if (length(cluster) != r) stop("Length of cluster is not appropriate with data. The length must be ", r) vardir = data[, vardir] cluster = as.matrix(data[, cluster]) samevar = samevar } else { formuladata = lapply(formula, function(x) model.frame(x, na.action = na.omit)) y = unlist(lapply(formula, function(x) model.frame(x, na.action = na.omit)[[1]])) X = Reduce(adiag, lapply(formula, function(x) model.matrix(x))) W = as.matrix(weight) n = length(y)/r if ((dim(vardir)[1] != n) || (dim(vardir)[2] != sum(1:r))) stop("Object vardir is not appropriate with data. It must be ", n ," x ", sum(1:r) ," matrix.") if (any(is.na(weight))) stop("Object weight contains NA values.") if ((dim(weight)[1] != n) || (dim(weight)[2] != r)) stop("Object weight is not appropriate with data. It must be ", n ," x ", r ," matrix.") if (any(is.na(cluster))) stop("Object cluster contains NA values.") if ((dim(cluster)[1] != n) || (dim(cluster)[2] != r)) stop("Object cluster is not appropriate with data. It must be ", n ," x ", r ," matrix.") cluster = as.matrix(cluster) samevar = samevar } y.matrix = matrix(as.vector(y), n, r) y.matrix[y.matrix==0] = NA indexns = unique(which(is.na(y.matrix), arr.ind=TRUE)[,1]) indexs = c(1:n)[-indexns] n.s = length(indexs) n.ns = n - n.s y.s = y.matrix[-indexns, ] y.s.vec = as.vector(y.s) W.s = W[-indexns, ] W.s = prop.table(W.s, 2) Xindexns = c() for (i in 1:r) { Xindexns = c(Xindexns, indexns + rep(n, times = n.ns)*(i - 1)) } X.s = X[-Xindexns, ] X.ns = X[Xindexns, ] y_names = sapply(formula, "[[", 2) temp_vardir = temporary_vardir(vardir[-indexns, ], n = n.s, r = r) eblup_first = eblup_inside(r = r, n = n, samevar = samevar, y = y, X = X, cluster = cluster, vardir = vardir) R.s = eblup_first$R.s
R.ns = eblup_first$R.ns beta = eblup_first$fit$estcoef[, 1] A = eblup_first$fit$refvar A_mat = kronecker(diag(as.vector(A)), diag(n.s)) Vinv = solve(A_mat + R.s) XtVinv = t(Vinv %*% X.s) Q = solve(XtVinv %*% X.s) mse_prasad = eblup_first$mse
g12 = eblup_first$g12 sumg12.pb = rep(0, n * r) sumg3.pb = rep(0, n * r) boot = 1 while (boot <= B) { u.boot = mvrnorm(n = n.s, mu = rep(0, r), Sigma = (diag(as.vector(A), nrow = r, ncol = r))) theta.boot = X.s %*% beta + as.vector(u.boot) e.boot = matrix(0, nrow = n.s, ncol = r) for (i in 1:n.s) { e.boot[i, ] = mvrnorm(1, mu = rep(0, r), Sigma = (temp_vardir[[i]])) } direct.boot = theta.boot + as.vector(e.boot) direct.boot.mat = matrix(direct.boot, nrow = n.s, ncol = r) direct.boot.mat.s = data.frame(index = indexs, direct.boot.mat) direct.boot.mat.ns = data.frame(index = indexns, matrix(0, n.ns, r)) direct.boot.mat.full = rbind(direct.boot.mat.s, direct.boot.mat.ns) direct.boot.mat.full = direct.boot.mat.full[order(direct.boot.mat.full$index), ]
rownames(direct.boot.mat.full) = direct.boot.mat.full$index direct.boot.mat.full = direct.boot.mat.full[, -1] resultEBLUP = eblup_inside(r = r, n = n, samevar = samevar, y = unlist(direct.boot.mat.full), X = X, cluster = cluster, vardir = vardir) sigma2.simula = resultEBLUP$fit$refvar beta.simula = resultEBLUP$fit$estcoef[, 1] Gn.simula = kronecker(diag(as.vector(sigma2.simula)), diag(n.s)) Vinv.simula = as.matrix(solve(Gn.simula + R.s)) Xbeta.simula = X %*% beta.simula Xsbeta.simula = X.s %*% beta.simula XtVi.simula = t(Vinv.simula %*% X.s) Q.simula = solve(XtVi.simula %*% X.s) #randomEffect random.effect.simula = Gn.simula %*% Vinv.simula %*% (direct.boot - Xsbeta.simula) random.effect.simula = matrix(random.effect.simula, nrow = n.s, ncol = r) random.effect.simula = data.frame(index = indexs, random.effect.simula) random.effect.simula.ns = data.frame(index = indexns, matrix(NA, n.ns, r)) names(random.effect.simula.ns) = names(random.effect.simula) random.effect.simula = rbind(random.effect.simula, random.effect.simula.ns) random.effect.simula = random.effect.simula[order(random.effect.simula$index), ]
rownames(random.effect.simula) = random.effect.simula$index random.effect.simula = random.effect.simula[ ,-1] random.effect.simula.full = matrix(NA, nrow = n, ncol = r) colnames(random.effect.simula.full) = y_names for (i in 1:r) { df = data.frame(cluster[, i], random.effect.simula[, i]) for (j in 1:n) { if (df[j, 2] %in% NA) { df[j, 2] = mean(df[df[, 1] == df[j,1], 2], na.rm = TRUE) } } random.effect.simula.full[ ,i] = df[, 2] } thetaEBLUP.boot1 = Xbeta.simula + as.vector(random.effect.simula.full) thetaEBLUP.boot1.mat = matrix(thetaEBLUP.boot1, nrow = n, ncol = r) thetaRATIO.boot1.mat = matrix(0, nrow = n, ncol = r) for (i in 1:r) { thetaRATIO.boot1.mat[, i] = thetaEBLUP.boot1.mat[, i] * (sum(W.s[, i] * direct.boot.mat[, i]) / sum(W[, i] * thetaEBLUP.boot1.mat[, i])) } g1boot.s = diag(Gn.simula %*% Vinv.simula %*% R.s) g2boot.s = diag(R.s %*% Vinv.simula %*% X.s %*% Q.simula %*% t(X.s) %*% t(R.s %*% Vinv.simula)) Gn.simula.ns = kronecker(diag(as.vector(sigma2.simula)), diag(n.ns)) Vinv.simula.ns = as.matrix(solve(Gn.simula.ns + R.ns)) XtVinv.simula.ns = t(Vinv.simula.ns %*% X.ns) Q.ns = ginv(XtVinv.simula.ns %*% X.ns) g1boot.ns = diag(Gn.simula.ns %*% Vinv.simula.ns %*% R.ns) g2boot.ns = diag(R.ns %*% Vinv.simula.ns %*% X.ns %*% Q.ns %*% t(X.ns) %*% t(R.ns %*% Vinv.simula.ns)) g12boot.s = matrix(g1boot.s + g2boot.s, n.s, r) names(g12boot.s) = y_names g12boot.s = data.frame(index = indexs, g12boot.s) g12boot.ns = matrix(g1boot.ns + g2boot.ns, n.ns, r) names(g12boot.ns) = y_names g12boot.ns = data.frame(index = indexns, g12boot.ns) g12boot = rbind(g12boot.s, g12boot.ns) g12boot = g12boot[order(g12boot$index), ]
rownames(g12boot) = g12boot$index g12boot = g12boot[, -1] Bstim.eblup = solve(XtVinv %*% X.s) %*% XtVinv %*% direct.boot Xbeta.eblup = X %*% Bstim.eblup Xsbeta.eblup = X.s %*% Bstim.eblup random.effect2.simula = A_mat %*% Vinv %*% (direct.boot - Xsbeta.eblup) random.effect2.simula = matrix(random.effect2.simula, nrow = n.s, ncol = r) random.effect2.simula = data.frame(index = indexs, random.effect2.simula) random.effect2.simula.ns = data.frame(index = indexns, matrix(NA, n.ns, r)) names(random.effect2.simula.ns) = names(random.effect2.simula) random.effect2.simula = rbind(random.effect2.simula, random.effect2.simula.ns) random.effect2.simula = random.effect2.simula[order(random.effect2.simula$index), ]
rownames(random.effect2.simula) = random.effect2.simula$index random.effect2.simula = random.effect2.simula[ ,-1] random.effect2.simula.full = matrix(NA, nrow = n, ncol = r) colnames(random.effect2.simula.full) = y_names for (i in 1:r) { df = data.frame(cluster[, i], random.effect2.simula[, i]) for (j in 1:n) { if (df[j, 2] %in% NA) { df[j, 2] = mean(df[df[, 1] == df[j,1], 2], na.rm = TRUE) } } random.effect2.simula.full[ ,i] = df[, 2] } thetaEBLUP.boot2 = Xbeta.eblup + as.vector(random.effect2.simula.full) thetaEBLUP.boot2.mat = matrix(thetaEBLUP.boot2, nrow = n, ncol = r) thetaRATIO.boot2.mat = matrix(0, nrow = n, ncol = r) for (i in 1:r) { thetaRATIO.boot2.mat[, i] = thetaEBLUP.boot2.mat[, i] * (sum(W.s[, i] * direct.boot.mat[, i]) / sum(W[, i] * thetaEBLUP.boot2.mat[, i])) } g3boot = (thetaRATIO.boot1.mat - thetaRATIO.boot2.mat)^2 sumg12.pb = sumg12.pb + as.vector(g12boot) sumg3.pb = sumg3.pb + as.vector(g3boot) boot = boot + 1 } g12.pb = sumg12.pb/B g3.pb = sumg3.pb/B msebootratio.df = 2 * (g12) - g12.pb + g3.pb names(mse_prasad) = y_names names(msebootratio.df) = y_names end_time = Sys.time() running_time = end_time - start_time result1 = list(mse.eblup = NA, pbmse.eblupRB = NA, running.time = NA) result1$mse.eblup = mse_prasad
result1$pbmse.eblupRB = msebootratio.df result1$running.time = running_time
return(result1)
}


## Try the msaeRB package in your browser

Any scripts or data that you put into this service are public.

msaeRB documentation built on June 13, 2021, 1:06 a.m.