vcov.rma <- function(object, type="fixed", ...) {
mstyle <- .get.mstyle()
.chkclass(class(object), must="rma")
na.act <- getOption("na.action")
on.exit(options(na.action=na.act), add=TRUE)
if (!is.element(na.act, c("na.omit", "na.exclude", "na.fail", "na.pass")))
stop(mstyle$stop("Unknown 'na.action' specified under options()."))
type <- match.arg(type, c("fixed", "beta", "alpha", "delta", "obs", "fitted", "resid"))
#########################################################################
if (type=="fixed") {
out <- object$vb
if (inherits(object, "rma.ls"))
out <- list(beta = object$vb, alpha = object$va)
if (inherits(object, "rma.uni.selmodel"))
out <- list(beta = object$vb, delta = object$vd)
return(out)
}
if (type=="beta") {
out <- object$vb
return(out)
}
if (type=="alpha") {
if (!inherits(object, "rma.ls"))
stop(mstyle$stop("Can only extract var-cov matrix of alpha coefficients for location-scale models."))
out <- object$va
return(out)
}
if (type=="delta") {
if (!inherits(object, "rma.uni.selmodel"))
stop(mstyle$stop("Can only extract var-cov matrix of delta coefficients for selection models."))
out <- object$vd
return(out)
}
#########################################################################
if (type=="obs") {
if (inherits(object, c("rma.uni","rma.mv"))) {
out <- matrix(NA_real_, nrow=object$k.f, ncol=object$k.f)
out[object$not.na, object$not.na] <- as.matrix(object$M) # as.matrix() needed when sparse=TRUE
rownames(out) <- colnames(out) <- object$slab
if (na.act == "na.omit")
out <- out[object$not.na, object$not.na]
if (na.act == "na.fail" && any(!object$not.na))
stop(mstyle$stop("Missing values in data."))
return(out)
} else {
stop(mstyle$stop("Extraction of marginal var-cov matrix not available for objects of this class."))
}
}
#########################################################################
if (type=="fitted") {
out <- object$X.f %*% object$vb %*% t(object$X.f)
rownames(out) <- colnames(out) <- object$slab
if (na.act == "na.omit")
out <- out[object$not.na, object$not.na]
if (na.act == "na.exclude" || na.act == "na.pass") {
out[!object$not.na,] <- NA_real_
out[,!object$not.na] <- NA_real_
}
return(out)
}
#########################################################################
if (type=="resid") {
### the SEs of the residuals cannot be estimated consistently for "robust.rma" objects
if (inherits(object, c("robust.rma", "rma.gen")))
stop(mstyle$stop("Extraction of var-cov matrix of the residuals not available for objects of this type."))
options(na.action="na.omit")
H <- hatvalues(object, type="matrix")
options(na.action = na.act)
ImH <- diag(object$k) - H
if (inherits(object, "robust.rma")) {
ve <- ImH %*% tcrossprod(object$meat,ImH)
} else {
ve <- ImH %*% tcrossprod(as.matrix(object$M),ImH) # as.matrix() needed when sparse=TRUE
}
if (na.act == "na.omit") {
out <- ve
rownames(out) <- colnames(out) <- object$slab[object$not.na]
}
if (na.act == "na.exclude" || na.act == "na.pass") {
out <- matrix(NA_real_, nrow=object$k.f, ncol=object$k.f)
out[object$not.na, object$not.na] <- ve
rownames(out) <- colnames(out) <- object$slab
}
return(out)
}
#########################################################################
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.