Nothing
#-------------------------------------
# vcovCR with defaults
#-------------------------------------
#' Cluster-robust variance-covariance matrix for a gls object.
#'
#' \code{vcovCR} returns a sandwich estimate of the variance-covariance matrix
#' of a set of regression coefficient estimates from a \code{\link[nlme]{gls}} object.
#'
#' @param cluster Optional expression or vector indicating which observations
#' belong to the same cluster. If not specified, will be set to
#' \code{getGroups(obj)}.
#' @param target Optional matrix or vector describing the working
#' variance-covariance model used to calculate the \code{CR2} and \code{CR4}
#' adjustment matrices. If not specified, the target is taken to be the
#' estimated variance-covariance structure of the \code{gls} object.
#' @inheritParams vcovCR
#'
#' @return An object of class \code{c("vcovCR","clubSandwich")}, which consists
#' of a matrix of the estimated variance of and covariances between the
#' regression coefficient estimates.
#'
#' @seealso \code{\link{vcovCR}}
#'
#' @examples
#'
#' if (requireNamespace("nlme", quietly = TRUE)) {
#'
#' library(nlme)
#' data(Ovary, package = "nlme")
#' Ovary$time_int <- 1:nrow(Ovary)
#' lm_AR1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data = Ovary,
#' correlation = corAR1(form = ~ time_int | Mare))
#' vcovCR(lm_AR1, type = "CR2")
#'
#' }
#'
#' @export
vcovCR.gls <- function(obj, cluster, type, target, inverse_var, form = "sandwich", ...) {
if (missing(cluster)) cluster <- nlme::getGroups(obj)
if (missing(target)) target <- NULL
if (missing(inverse_var) ) inverse_var <- is.null(target)
vcov_CR(obj, cluster = cluster, type = type,
target = target, inverse_var = inverse_var, form = form)
}
# residuals_CS()
# coef()
# nobs()
#-------------------------------------
# model_matrix()
#-------------------------------------
get_data <- function (object) {
if ("data" %in% names(object)) {
data <- object$data
} else {
dat_call <- object$call$data
envir_names <- sys.frames()
data <- simpleError("start")
i <- 1L
while (inherits(data, "simpleError") & i <= length(envir_names)) {
data <- tryCatch(eval(dat_call, envir = envir_names[[i]]), error = function(e) e)
i <- i + 1L
}
}
if (inherits(data, "simpleError")) return(NULL)
naAct <- object[["na.action"]]
if (!is.null(naAct)) {
data <- if (inherits(naAct, "omit")) {
data[-naAct, ]
} else if (inherits(naAct, "exclude")) {
data
} else eval(object$call$na.action)(data)
}
subset <- object$call$subset
if (!is.null(subset)) {
subset <- eval(asOneSidedFormula(subset)[[2]], data)
data <- data[subset, ]
}
data
}
#' @export
model_matrix.gls <- function(obj) {
dat <- get_data(obj)
model.matrix(formula(obj), data = dat)
}
#-------------------------------------
# Get (model-based) working variance matrix
#-------------------------------------
#' @export
targetVariance.gls <- function(obj, cluster = nlme::getGroups(obj)) {
groups <- nlme::getGroups(obj)
if (is.null(groups)) groups <- cluster
if (is.null(obj$modelStruct$corStruct)) {
if (is.null(obj$modelStruct$varStruct)) {
V_list <- matrix_list(rep(1, length(cluster)), cluster, "both")
} else {
wts <- nlme::varWeights(obj$modelStruct$varStruct)
V_list <- matrix_list(1 / wts^2, cluster, "both")
}
} else {
R_list <- nlme::corMatrix(obj$modelStruct$corStruct)
if (is.null(obj$modelStruct$varStruct)) {
V_list <- R_list
} else {
sd_vec <- 1 / nlme::varWeights(obj$modelStruct$varStruct)[order(order(groups))]
sd_list <- split(sd_vec, groups)
V_list <- Map(function(R, s) tcrossprod(s) * R, R = R_list, s = sd_list)
}
}
# check if clustering level is higher than highest level of random effects
tb_groups <- table(groups)
tb_cluster <- table(cluster)
if (length(tb_groups) < length(tb_cluster) |
any(as.vector(tb_groups) != rep(as.vector(tb_cluster), length.out = length(tb_groups))) |
any(names(tb_groups) != rep(names(tb_cluster), length.out = length(tb_groups)))) {
# check that random effects are nested within clusters
tb_cross <- table(groups, cluster)
nested <- apply(tb_cross, 1, function(x) sum(x > 0) == 1)
if (!all(nested)) stop("Random effects are not nested within clustering variable.")
# expand target_list to level of clustering
crosswalk <- data.frame(groups, cluster)
V_list <- add_bdiag(small_mats = V_list,
big_mats = matrix_list(rep(0, length(cluster)), cluster, dim = "both"),
crosswalk = crosswalk)
}
V_list
}
#-------------------------------------
# Get weighting matrix
#-------------------------------------
#' @export
weightMatrix.gls <- function(obj, cluster = nlme::getGroups(obj)) {
V_list <- targetVariance(obj, cluster)
lapply(V_list, function(v) chol2inv(chol(v)))
}
#---------------------------------------
# Get bread matrix and scaling constant
#---------------------------------------
#' @export
bread.gls <- function(x, ...) {
vcov(x) * nobs(x) / x$sigma^2
}
# v_scale() is default
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.