From 9b351e78609ff97f06cb67ffacead376ee845011 Mon Sep 17 00:00:00 2001 From: Mauricio 'Pacha' Vargas Sepulveda Date: Sun, 8 Sep 2024 03:31:18 -0400 Subject: [PATCH 1/2] consistent generics --- R/generics_augment.R | 2 +- man/broom.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/generics_augment.R b/R/generics_augment.R index c731f2a..d970367 100644 --- a/R/generics_augment.R +++ b/R/generics_augment.R @@ -5,7 +5,7 @@ generics::augment #' @title Broom Integration #' #' @description The provided `broom` methods do the following: -#' 1. `augment`: takes the input data and adds additional columns with the +#' 1. `augment`: Takes the input data and adds additional columns with the #' fitted values and residuals. #' 2. `glance`: Extracts the deviance, null deviance, and the number of #' observations.` diff --git a/man/broom.Rd b/man/broom.Rd index 0014c1b..8eb1e58 100644 --- a/man/broom.Rd +++ b/man/broom.Rd @@ -42,7 +42,7 @@ and \code{tidy} methods. \description{ The provided \code{broom} methods do the following: \enumerate{ -\item \code{augment}: takes the input data and adds additional columns with the +\item \code{augment}: Takes the input data and adds additional columns with the fitted values and residuals. \item \code{glance}: Extracts the deviance, null deviance, and the number of observations.` From c3317f215a5e785b497a37e9681ea6f3c6d8defa Mon Sep 17 00:00:00 2001 From: Mauricio 'Pacha' Vargas Sepulveda Date: Sun, 8 Sep 2024 03:33:50 -0400 Subject: [PATCH 2/2] fix conflicting generics vcov --- R/generics_vcov.R | 233 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 233 insertions(+) create mode 100644 R/generics_vcov.R diff --git a/R/generics_vcov.R b/R/generics_vcov.R new file mode 100644 index 0000000..b3194e7 --- /dev/null +++ b/R/generics_vcov.R @@ -0,0 +1,233 @@ +#' @title Covariance matrix for APEs +#' +#' @description Covariance matrix for the estimator of the +#' average partial effects from objects returned by \code{\link{apes}}. +#' +#' @param object an object of class \code{"apes"}. +#' @param ... additional arguments. +#' +#' @return A named matrix of covariance estimates. +#' +#' @seealso \code{\link{apes}} +#' +#' @export +#' +#' @noRd +vcov.apes <- function(object, ...) { + object[["vcov"]] +} + +#' @title Covariance matrix for GLMs +#' +#' @description Covariance matrix for the estimator of the structural parameters +#' from objects returned by \code{\link{feglm}}. The covariance is computed +#' from the hessian, the scores, or a combination of both after convergence. +#' +#' @param object an object of class \code{"feglm"}. +#' @param type the type of covariance estimate required. \code{"hessian"} refers +#' to the inverse of the negative expected hessian after convergence and is the +#' default option. \code{"outer.product"} is the outer-product-of-the-gradient +#' estimator. \code{"sandwich"} is the sandwich estimator (sometimes also +#' referred as robust estimator), and \code{"clustered"} computes a clustered +#' covariance matrix given some cluster variables. +#' +#' @param ... additional arguments. +#' +#' @return A named matrix of covariance estimates. +#' +#' @references Cameron, C., J. Gelbach, and D. Miller (2011). "Robust Inference +#' With Multiway Clustering". Journal of Business & Economic Statistics 29(2). +#' +#' @seealso \code{\link{feglm}} +#' +#' @examples +#' # same as the example in feglm but extracting the covariance matrix +#' +#' # subset trade flows to avoid fitting time warnings during check +#' set.seed(123) +#' trade_2006 <- trade_panel[trade_panel$year == 2006, ] +#' trade_2006 <- trade_2006[sample(nrow(trade_2006), 500), ] +#' +#' mod <- fepoisson( +#' trade ~ log_dist + lang + cntg + clny | exp_year + imp_year | pair, +#' trade_2006 +#' ) +#' +#' round(vcov(mod, type = "clustered"), 5) +#' +#' @return A named matrix of covariance estimates. +#' +#' @export +vcov.feglm <- function( + object, + type = c("hessian", "outer.product", "sandwich", "clustered"), + ...) { + # Check validity of input argument 'type' + type <- match.arg(type) + + # Extract cluster from formula + # it is totally fine not to have a cluster variable + cl_vars <- vcov_feglm_vars_(object) + k <- length(cl_vars) + + if (isTRUE(k >= 1L) && type != "clustered") { + type <- "clustered" + } + + # Compute requested type of covariance matrix + h <- object[["hessian"]] + p <- ncol(h) + + if (type == "hessian") { + # If the hessian is invertible, compute its inverse + v <- vcov_feglm_hessian_covariance_(h, p) + } else { + g <- get_score_matrix_(object) + if (type == "outer.product") { + # Check if the OP is invertible and compute its inverse + v <- vcov_feglm_outer_covariance_(g, p) + } else { + v <- vcov_feglm_covmat_( + object, type, h, g, + cl_vars, k, p + ) + } + } + + v +} + +vcov_feglm_vars_ <- function(object) { + suppressWarnings({ + attr(terms(object[["formula"]], rhs = 3L), "term.labels") + }) +} + +vcov_feglm_hessian_covariance_ <- function(h, p) { + v <- try(solve(h), silent = TRUE) + if (inherits(v, "try-error")) { + v <- matrix(Inf, p, p) + } + v +} + +vcov_feglm_outer_covariance_ <- function(g, p) { + v <- try(solve(g), silent = TRUE) + if (inherits(v, "try-error")) { + v <- matrix(Inf, p, p) + } + v +} + +vcov_feglm_covmat_ <- function( + object, type, h, g, + cl_vars, k, p) { + # Check if the hessian is invertible and compute its inverse + v <- try(solve(h), silent = TRUE) + if (inherits(v, "try-error")) { + v <- matrix(Inf, p, p) + } else { + # Compute inner part of the sandwich formula + if (type == "sandwich") { + b <- crossprod(g) + } else { + if (isFALSE(k >= 1L)) { + vcov_feglm_cluster_nocluster_() + } + d <- vcov_feglm_cluster_data_(object, cl_vars) + d <- mutate(d, across(all_of(cl_vars), check_factor_)) + sp_vars <- colnames(g) + g <- cbind(d, g) + rm(d) + b <- vcov_feglm_clustered_cov_(g, cl_vars, sp_vars, p) + } + # Sandwich formula + v <- v %*% b %*% v + } + + # Return covariance estimate + v +} + +vcov_feglm_cluster_nocluster_ <- function() { + stop( + paste( + "No cluster variable was found.", + "Please specify a cluster variable", + "in the model formula." + ), + call. = FALSE + ) +} + +vcov_feglm_cluster_data_ <- function(object, cl_vars) { + d <- try(object[["data"]][, get("cl_vars"), with = FALSE], silent = TRUE) + if (inherits(d, "try-error")) { + vcov_feglm_cluster_notfound_() + } + d +} + +vcov_feglm_cluster_notfound_ <- function() { + stop( + paste( + "At least one cluster variable was not found.", + "Ensure to pass variables that are not part of the model", + "itself, but are required to compute clustered standard errors", + "to 'feglm'. This can be done via 'formula'. See documentation", + "for details." + ), + call. = FALSE + ) +} + +# Ensure cluster variables are factors ---- + +vcov_feglm_clustered_cov_ <- function(g, cl_vars, sp_vars, p) { + # Multiway clustering by Cameron, Gelbach, and Miller (2011) + b <- matrix(0.0, p, p) + for (i in seq.int(length(cl_vars))) { + # Compute outer product for all possible combinations + cl_combn <- combn(cl_vars, i) + br <- matrix(0.0, p, p) + for (j in seq.int(ncol(cl_combn))) { + cl <- cl_combn[, j] + br <- br + crossprod( + as.matrix( + g %>% + group_by(!!sym(cl)) %>% + summarise(across(all_of(sp_vars), sum), .groups = "drop") %>% + select(-!!sym(cl)) + ) + ) + } + + # Update outer product + if (i %% 2L) { + b <- b + br + } else { + b <- b - br + } + } + b +} + +#' @title Covariance matrix for LMs +#' +#' @description Covariance matrix for the estimator of the structural parameters +#' from objects returned by \code{\link{felm}}. The covariance is computed +#' from the hessian, the scores, or a combination of both after convergence. +#' +#' @param object an object of class \code{"felm"}. +#' +#' @inherit vcov.feglm +#' +#' @seealso \code{\link{felm}} +#' +#' @export +vcov.felm <- function( + object, + type = c("hessian", "outer.product", "sandwich", "clustered"), + ...) { + vcov.feglm(object, type) +}