dedicated lm function

Author: pachadotdev

DESCRIPTION

@@ -2,7 +2,7 @@ Package: capybara
 Type: Package
 Title: Fast and Memory Efficient Fitting of Linear Models With High-Dimensional
     Fixed Effects
-Version: 0.7.0
+Version: 0.8.0
 Authors@R: c(
     person(
         given = "Mauricio",

NEWS.md

@@ -1,3 +1,8 @@
+# capybara 0.8.0
+
+* Dedicated functions for linear models to avoid the overhead of running
+  the GLM function with a Gaussian link.
+
 # capybara 0.7.0
 
 * The predict method now allows to pass new data to predict the outcome.

R/autoplot.r renamed to R/autoplot.R

@@ -56,8 +56,8 @@ autoplot.feglm <- function(object, ...) {
   }
 
   # stop if the object is not of class feglm or felm
-  if (!inherits(object, "feglm") && !inherits(object, "felm")) {
-    stop("The object must be of class 'feglm' or 'felm'")
+  if (!inherits(object, "feglm")) {
+    stop("The object must be of class 'feglm'")
   }
 
   # if conf_level is not provided, set it to 0.95
@@ -137,5 +137,64 @@ autoplot.feglm <- function(object, ...) {
 #'
 #' @export
 autoplot.felm <- function(object, ...) {
-  autoplot.feglm(object, ...)
+  # stop if ggplot2 is not installed
+  if (!requireNamespace("ggplot2", quietly = TRUE)) {
+    stop("The 'ggplot2' package is required to use this function")
+  }
+
+  # stop if the object is not of class feglm or felm
+  if (!inherits(object, "felm")) {
+    stop("The object must be of class 'felm'")
+  }
+
+  # if conf_level is not provided, set it to 0.95
+  if (!"conf_level" %in% names(list(...))) {
+    conf_level <- 0.95
+  } else {
+    conf_level <- list(...)$conf_level
+  }
+
+  # check that conf_level is between 0 and 1
+  if (conf_level <= 0 || conf_level >= 1) {
+    stop("The confidence level must be between 0 and 1")
+  }
+
+  # Extract the summary of the felm object
+  res <- summary(object)$cm
+  colnames(res) <- c("estimate", "std.error", "statistic", "p.value")
+
+  # Calculate the critical value for the specified confidence conf_level
+  alpha <- 1 - conf_level
+  z <- qnorm(1 - alpha / 2)
+
+  # Compute the confidence intervals
+  conf_data <- data.frame(
+    term = rownames(res),
+    estimate = res[, "estimate"],
+    conf_low = res[, "estimate"] - z * res[, "std.error"],
+    conf_high = res[, "estimate"] + z * res[, "std.error"]
+  )
+
+  p <- ggplot(conf_data, aes(x = !!sym("term"), y = !!sym("estimate"))) +
+    geom_errorbar(
+      aes(
+        ymin = !!sym("conf_low"),
+        ymax = !!sym("conf_high")
+      ),
+      width = 0.1,
+      color = "#165976"
+    ) +
+    geom_point() +
+    labs(
+      title = sprintf(
+        "Coefficient Estimates with Confidence Intervals at %s%%",
+        round(conf_level * 100, 0)
+      ),
+      x = "Term",
+      y = "Estimate"
+    ) +
+    theme_minimal() +
+    coord_flip()
+
+  p
 }

R/cpp11.R

@@ -31,3 +31,7 @@ feglm_fit_ <- function(beta_r, eta_r, y_r, x_r, wt_r, theta, family, control, k_
 feglm_offset_fit_ <- function(eta_r, y_r, offset_r, wt_r, family, control, k_list) {
   .Call(`_capybara_feglm_offset_fit_`, eta_r, y_r, offset_r, wt_r, family, control, k_list)
 }
+
+felm_fit_ <- function(y_r, x_r, wt_r, control, k_list) {
+  .Call(`_capybara_felm_fit_`, y_r, x_r, wt_r, control, k_list)
+}

R/feglm_helpers.R

@@ -511,7 +511,7 @@ get_index_list_ <- function(k_vars, data) {
 #' @description Computes the score matrix
 #' @param object Result list
 #' @noRd
-get_score_matrix_ <- function(object) {
+get_score_matrix_feglm_ <- function(object) {
   # Extract required quantities from result list
   control <- object[["control"]]
   data <- object[["data"]]

R/felm.R

@@ -1,5 +1,5 @@
 #' srr_stats
-#' @srrstats {G1.0} Implements a wrapper around `feglm` for linear models with high-dimensional fixed effects.
+#' @srrstats {G1.0} Implements linear models with high-dimensional fixed effects.
 #' @srrstats {G2.1a} Ensures the input `formula` is correctly specified and includes fixed effects.
 #' @srrstats {G2.1b} Validates that the input `data` is non-empty and of class `data.frame`.
 #' @srrstats {G2.3a} Uses structured checks for parameters like `weights` and starting values.
@@ -98,21 +98,87 @@ NULL
 #' summary(mod)
 #'
 #' @export
-felm <- function(formula = NULL, data = NULL, weights = NULL) {
-  # Use 'feglm' to estimate the model
-  # Using felm_fit_ directly leads to the incorrect yhat = Xb
-  # we need iteratively reweighted least squares
-  reslist <- feglm(
-    formula = formula, data = data, weights = weights, family = gaussian()
-  )
-
-  names(reslist)[which(names(reslist) == "eta")] <- "fitted.values"
-
-  reslist[["conv"]] <- NULL
-  reslist[["iter"]] <- NULL
-  reslist[["family"]] <- NULL
-  reslist[["deviance"]] <- NULL
-
-  # Return result list
-  structure(reslist, class = "felm")
+felm <- function(formula = NULL, data = NULL, weights = NULL, control = NULL) {
+  # Check validity of formula ----
+  check_formula_(formula)
+
+  # Check validity of data ----
+  check_data_(data)
+
+  # Check validity of control + Extract control list ----
+  control <- check_control_(control)
+
+  # Update formula and do further validity check ----
+  formula <- update_formula_(formula)
+
+  # Generate model.frame
+  lhs <- NA # just to avoid global variable warning
+  nobs_na <- NA
+  nobs_full <- NA
+  model_frame_(data, formula, weights)
+
+  # Get names of the fixed effects variables and sort ----
+  k_vars <- attr(terms(formula, rhs = 2L), "term.labels")
+
+  # Generate temporary variable ----
+  tmp_var <- temp_var_(data)
+
+  # Transform fixed effects and clusters to factors ----
+  data <- transform_fe_(data, formula, k_vars)
+
+  # Determine the number of dropped observations ----
+  nt <- nrow(data)
+  nobs <- nobs_(nobs_full, nobs_na, nt)
+
+  # Extract model response and regressor matrix ----
+  nms_sp <- NA
+  p <- NA
+  model_response_(data, formula)
+
+  # Check for linear dependence ----
+  # check_linear_dependence_(x, p)
+  check_linear_dependence_(cbind(y,x), p + 1L)
+
+  # Extract weights if required ----
+  if (is.null(weights)) {
+    wt <- rep(1.0, nt)
+  } else {
+    wt <- data[[weights]]
+  }
+
+  # Check validity of weights ----
+  check_weights_(wt)
+
+  # Get names and number of levels in each fixed effects category ----
+  nms_fe <- lapply(select(data, all_of(k_vars)), levels)
+  lvls_k <- vapply(nms_fe, length, integer(1))
+
+  # Generate auxiliary list of indexes for different sub panels ----
+  k_list <- get_index_list_(k_vars, data)
+
+  # Fit linear model ----
+  if (is.integer(y)) {
+    y <- as.numeric(y)
+  }
+  fit <- felm_fit_(y, x, wt, control, k_list)
+
+  y <- NULL
+  x <- NULL
+
+  # Add names to beta, hessian, and mx (if provided) ----
+  names(fit[["coefficients"]]) <- nms_sp
+  if (control[["keep_mx"]]) {
+    colnames(fit[["mx"]]) <- nms_sp
+  }
+
+  # Add to fit list ----
+  fit[["nobs"]] <- nobs
+  fit[["lvls_k"]] <- lvls_k
+  fit[["nms_fe"]] <- nms_fe
+  fit[["formula"]] <- formula
+  fit[["data"]] <- data
+  fit[["control"]] <- control
+
+  # Return result list ----
+  structure(fit, class = "felm")
 }

R/felm_helpers.R

@@ -0,0 +1,42 @@
+#' srr_stats
+#' @srrstats {G1.0} Provides modular helper functions for computations in linear models with fixed effects.
+#' @noRd
+NULL
+
+#' @title Get score matrix
+#' @description Computes the score matrix
+#' @param object Result list
+#' @noRd
+get_score_matrix_felm_ <- function(object) {
+  # Extract required quantities from result list
+  control <- object[["control"]]
+  data <- object[["data"]]
+  w <- object[["weights"]]
+
+  # Update weights and dependent variable
+  y <- data[[1L]]
+
+  # Center regressor matrix (if required)
+  if (control[["keep_mx"]]) {
+    mx <- object[["mx"]]
+  } else {
+    # Extract additional required quantities from result list
+    formula <- object[["formula"]]
+    k_vars <- names(object[["lvls_k"]])
+
+    # Generate auxiliary list of indexes to project out the fixed effects
+    k_list <- get_index_list_(k_vars, data)
+
+    # Extract regressor matrix
+    x <- model.matrix(formula, data, rhs = 1L)[, -1L, drop = FALSE]
+    nms_sp <- attr(x, "dimnames")[[2L]]
+    attr(x, "dimnames") <- NULL
+
+    # Center variables
+    mx <- center_variables_r_(x, w, k_list, control[["center_tol"]], 10000L)
+    colnames(mx) <- nms_sp
+  }
+
+  # Return score matrix
+  mx * (y * w)
+}