#' Sensitivity Analysis for Correlated Effects RVE
#'
#' \code{sensitivity} is used to assess the impact of differing rho values on
#' the correlated effects meta-regression model.
#'
#'
#' @aliases sensitivity
#' @param x A dataframe containing values of rho, tau squared, coefficient
#' estimates, and standard errors.
#' @references
#'
#' Hedges, L.V., Tipton, E., Johnson, M.C. (2010) Robust variance estimation in
#' meta-regression with dependent effect size estimates. \emph{Research
#' Synthesis Methods}. \bold{1}(1): 39--65. Erratum in \bold{1}(2): 164--165.
#' DOI: 10.1002/jrsm.5
#'
#' Tipton, E. (in press) Small sample adjustments for robust variance
#' estimation with meta-regression. \emph{Psychological Methods}.
#' @keywords robu
#' @examples
#'
#'
#' # Correlated Effects Model
#' CorrMod <- robu(formula = effectsize ~ followup + males + binge + college,
#' data = corrdat, studynum = studyid, var.eff.size = var,
#' rho = .8, modelweights = "CORR", small = FALSE)
#'
#' sensitivity(CorrMod) # Output sensitivity
#'
#' @export
sensitivity <- function(x){
modelweights <- x$modelweights
user_weighting <- x$user_weighting
if(modelweights == "HIER")
stop("Sensitivity analysis is not available for hierarchical effects.")
if(user_weighting == TRUE)
stop("Sensitivity analysis is not available for user specified weights.")
mod_info <- x$mod_info
p <- x$p
N <- x$N
Xreg <- x$Xreg
y <- x$y
X <- x$X
data.full <- x$data.full
X.full <- x$X.full
k <- data.full$k
k_list <- x$k_list
ml <- x$ml
term1 <- mod_info$term1
term2 <- mod_info$term2
small <- x$small
labels <- x$labels
mod_label <- x$mod_label
rho.test <- seq(0, 1, .2)
rho_labels <- c(paste("Rho = ", seq(0, 1, .2), sep=""))
var_labels <- rep("", 2 * (p + 1))
var_labels[seq(1, length(var_labels), by = 2)] <- labels
var_labels <- c(var_labels, "Tau.sq")
col2_labels <- rep("", 2 * (p + 1))
col2_labels[seq(1, length(col2_labels), by = 2)] <- "Coefficient"
col2_labels[seq(2, length(col2_labels), by = 2)] <- "Std. Error"
col2_labels <- c(col2_labels, "Estimate")
sen <- data.frame(cbind(var_labels, col2_labels))
for (i in (1: length(rho.test))){
tau.sq1 <- term1 + rho.test[i] * term2
tau.sq <- ifelse(tau.sq1 < 0, 0, tau.sq1)
data.full$r.weights <- 1 / (as.vector(data.full$k) *
(as.vector(data.full$avg.var.eff.size)
+ as.vector(tau.sq)))
W.r.big <- diag(data.full$r.weights) # W
W.r <- by(data.full$r.weights, data.full$study,
function(x) diag(x, nrow = length(x)))
sumXWX.r <- Reduce("+", Map(function(X, W)
t(X) %*% W %*% X,
X, W.r))
sumXWy.r <- Reduce("+", Map(function(X, W, y)
t(X) %*% W %*% y,
X, W.r, y))
b.r <- solve(sumXWX.r) %*% sumXWy.r
data.full$pred.r <- Xreg %*% b.r
data.full$e.r <- cbind(data.full$effect.size) -
data.full$pred.r
data.full$e.r <- as.numeric(data.full$e.r)
sigma.hat.r <- by(data.full$e.r, data.full$study,
function(x) tcrossprod(x))
if (!small) { # Begin small = FALSE
sumXWeeWX.r <- Reduce("+", Map(function(X, W, V)
t(X) %*% W %*% V %*% W %*% X,
X, W.r, sigma.hat.r))
VR.r <- solve(sumXWX.r) %*% sumXWeeWX.r %*%
solve(sumXWX.r)
SE <- sqrt(diag(VR.r)) * sqrt(N / (N - (p + 1)))
} else {
Q <- solve(sumXWX.r) #
Q.list <- rep(list(Q), N)
H <- Xreg %*% Q %*% t(Xreg) %*% W.r.big
ImH <- diag(c(1), dim(Xreg)[1], dim(Xreg)[1]) - H
data.full$ImH <- cbind(ImH)
ImHj <- by(data.full$ImH, data.full$study,
function(x) as.matrix(x))
dfS <- c(rep(0, p + 1))
diag_one <- by(rep(1, nrow(X.full)), X.full$study,
function(x) diag(x, nrow = length(x)))
ImHii <- Map(function(X, Q, W, D)
D - X %*% Q %*% t(X) %*% W,
X, Q.list, W.r, diag_one)
eigenvec <- lapply(ImHii, function(x) eigen(x)$vectors)
eigenval <- lapply(ImHii, function(x) eigen(x)$values)
I <- ImHii
A.MBB <- Map(function (eigenvec, eigenval, k_list)
eigenvec %*% diag(1/sqrt(eigenval),
k_list, k_list)
%*% t(eigenvec),
eigenvec, eigenval, k_list)
A.MBB1 <- Map(function(K, A, I)
if (K > 1) A else matrix(sqrt(solve(I))),
k_list, A.MBB, I)
A.MBB2 <- A.MBB
sumXWA.MBBeeA.MBBWX.r <- Reduce("+", Map(function(X,W,A,S)
t(X) %*% W %*% A %*%
S %*% A %*% W %*% X,
X, W.r, A.MBB2,
sigma.hat.r))
data.full$ImH <- ImH
ImH <- lapply(split(data.full$ImH,
data.full$study),
matrix, ncol=nrow(data.full))
giTemp <- Map(function(I, A, W, X, Q)
t(I) %*% A %*% W %*% X %*% Q,
ImHj, A.MBB2, W.r, X, Q.list)
dfs <- c(rep(0, p + 1))
for (i in 1:(p + 1)) {
L <- c(rep(0,p + 1))
L[i] <- 1
Ll <- rep(list(L), N)
gi <- Map(function(G, L) G %*% cbind(L), giTemp, Ll)
G <- Reduce("+", lapply(gi, function(x) tcrossprod(x)))
B <- solve(sqrt(W.r.big) )%*% G %*% solve(sqrt(W.r.big))
e.val2 <- eigen(B)
dfs[i] <- sum(e.val2$values)^2 / sum(e.val2$values^2)
}
VR.MBB1 <- solve(sumXWX.r) %*% sumXWA.MBBeeA.MBBWX.r %*%
solve(sumXWX.r)
VR.r <- VR.MBB1
SE <- sqrt(diag(VR.r))
}
vals <- c()
temp_vals <- c()
for (i in 1:(p + 1)) {
temp_vals <- c(b.r[i], SE[i])
vals <- c(vals, temp_vals)
}
vals <- c(vals, tau.sq)
vals <- format(vals, digits=3, justify="centre")
sen <- cbind(sen, vals)
}
colnames(sen) <- c(" ", " ", rho_labels)
format(sen[,1], justify = "left")
cat(mod_label, "\n")
cat("Model:",paste(x$ml[2]), paste(x$ml[[1]]), paste(x$ml[3]),"\n\n")
cat(paste("Sensitivity Analysis"), "\n\n")
print.data.frame(sen, quote = FALSE, row.names = FALSE, right = FALSE)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.