R/get_APEs.R

In bife: Binary Choice Models with Fixed Effects

#' @title
#' Compute average partial effects for binary choice models with fixed effects
#' @description
#' \code{\link{get_APEs}} is a post-estimation routine that can be used to estimate average partial 
#' effects with respect to all covariates in the model and the corresponding covariance matrix. The 
#' estimation of the covariance is based on a linear approximation (delta method). Note that 
#' the command automatically determines which of the regressors are continuous or binary.
#' 
#' \strong{Remark:} The routine currently does not allow to compute average partial effects based 
#' on functional forms like interactions and polynomials.
#' 
#' \strong{Note:} \code{\link{apeff_bife}} is deprecated and will be removed soon.
#' @param
#' object an object of class \code{"bife"}.
#' @param
#' n_pop unsigned integer indicating a finite population correction for the estimation of the 
#' covariance matrix of the average partial effects proposed by 
#' Cruz-Gonzalez, Fernández-Val, and Weidner (2017). The correction factor is computed as follows: 
#' \eqn{(n^{\ast} - n) / (n^{\ast} - 1)}{(n.pop - n) / (n.pop - 1)}, 
#' where \eqn{n^{\ast}}{n.pop} and \eqn{n}{n} are the size of the entire 
#' population and the full sample size. Default is \code{NULL}, which refers to a factor of zero and 
#' a covariance obtained by the delta method.
#' @param
#' sampling_fe a string equal to \code{"independence"} or \code{"unrestricted"} which imposes 
#' sampling assumptions about the unobserved effects. \code{"independence"} imposes that all 
#' unobserved effects are mutually independent sequences. \code{"unrestricted"} does not impose any
#' sampling assumptions. Note that this option only affects the optional finite population correction. 
#' Default is \code{"independence"}.
#' @param
#' weak_exo logical indicating if some of the regressors are assumed to be weakly exogenous (e.g. 
#' predetermined). If object is returned by \code{\link{bias_corr}}, the option will be 
#' automatically set to \code{TRUE} if the chosen bandwidth parameter is larger than zero. 
#' Note that this option only affects the estimation of the covariance matrix. 
#' Default is \code{FALSE}, which assumes that all regressors are strictly exogenous.
#' @param
#' ... arguments passed to the deprecated function \code{\link{apeff_bife}}.
#' @return
#' The function \code{\link{get_APEs}} returns a named list of class \code{"bifeAPEs"}.
#' @references
#' Cruz-Gonzalez, M., I. Fernández-Val, and M. Weidner. (2017). "Bias corrections for probit and 
#' logit models with two-way fixed effects". The Stata Journal, 17(3), 517-545.
#' @references
#' Fernández-Val, I. (2009). "Fixed effects estimation of structural parameters and marginal 
#' effects in panel probit models". Journal of Econometrics 150(1), 71-85.
#' @references
#' Fernández-Val, I. and M. Weidner (2018). "Fixed effects estimation of large-t panel data models". 
#' Annual Review of Economics, 10, 109-138.
#' @references
#' Neyman, J. and E. L. Scott (1948). "Consistent estimates based on partially consistent 
#' observations". Econometrica, 16(1), 1-32.
#' @references
#' Stammann, A., F. Heiss, and D. McFadden (2016). "Estimating Fixed Effects Logit Models with 
#' Large Panel Data". Working paper.
#' @seealso
#' \code{\link{bias_corr}}, \code{\link{bife}}
#' @examples 
#' \donttest{
#' # Load 'psid' dataset
#' library(bife)
#' dataset <- psid
#' 
#' # Fit a static logit model
#' mod <- bife(LFP ~ I(AGE^2) + log(INCH) + KID1 + KID2 + KID3 + factor(TIME) | ID, dataset)
#' summary(mod)
#' 
#' # Compute average partial effects
#' mod_ape <- get_APEs(mod)
#' summary(mod_ape)
#' 
#' # Apply analytical bias correction
#' mod_bc <- bias_corr(mod)
#' summary(mod_bc)
#' 
#' # Compute bias-corrected average partial effects
#' mod_ape_bc <- get_APEs(mod_bc)
#' summary(mod_ape_bc)
#' }
#' @export
get_APEs <- function(
    object,
    n_pop       = NULL,
    sampling_fe = c("independence", "unrestricted"),
    weak_exo    = FALSE
    ) {
  # Validity of input argument (object)
  if (!inherits(object, "bife")) {
    stop("'get_APEs' called on a non-'bife' object.", call. = FALSE)
  }
  
  # Check validity of 'sampling_fe'
  sampling_fe <- match.arg(sampling_fe)
  
  # Extract prior information if available
  biascorr <- !is.null(object[["bandwidth"]])
  if (biascorr) {
    L <- object[["bandwidth"]]
    if (L > 0L) {
      weak_exo <- TRUE
    } else {
      weak_exo <- FALSE
    }
  }
  
  # Get names of the fixed effects variables and extract individual number of time periods
  idvar <- attr(terms(object[["formula"]], rhs = 2L), "term.labels")
  Ti <- object[["data"]][, .N, by = eval(idvar)][[2L]]
  
  # Extract model information
  beta <- object[["coefficients"]]
  family <- object[["family"]]
  nt_full <- object[["nobs"]][["nobs_full"]]
  nt <- object[["nobs"]][["nobs"]]
  p <- length(beta)
  
  # Check validity of 'n_pop'
  if (!is.null(n_pop)) {
    n_pop <- as.integer(n_pop)
    if (n_pop < nt_full) {
      warning(paste("Size of the entire population is lower than the full sample size.",
                    "Correction factor set to zero."), call. = FALSE)
      adj <- 0.0
    } else {
      adj <- (n_pop - nt_full) / (n_pop - 1L)
    }
  } else {
    adj <- 0.0
  }
  
  # Extract model response and regressor matrix
  y <- object[["data"]][[1L]]
  X <- model.matrix(object[["formula"]], object[["data"]], rhs = 1L)[, - 1L, drop = FALSE]
  id <- as.integer(object[["data"]][[idvar]])
  attr(X, "dimnames") <- NULL
  
  # Determine which of the regressors are binary
  binary <- apply(X, 2L, function(x) all(x %in% c(0L, 1L)))
  
  # Compute required derivatives
  eta <- as.vector(X %*% beta + object[["alpha"]][id])
  mu <- family[["linkinv"]](eta)
  mu_eta <- family[["mu.eta"]](eta)
  if (family[["link"]] == "logit") {
    v <- y - mu
    w <- mu_eta
    z <- w * (1.0 - 2.0 * mu)
  } else {
    w <- mu_eta / family[["variance"]](mu)
    v <- w * (y - mu)
    w <- w * mu_eta
    z <- - eta * w
  }
  rm(y)
  
  # Compute average partial effects and Jacobian
  MX <- center_variables(X * sqrt(w), sqrt(w), Ti) / sqrt(w)
  PX <- X - MX
  Delta <- matrix(NA_real_, nt, p)
  Delta1 <- matrix(NA_real_, nt, p)
  J <- matrix(NA_real_, p, p)
  Delta[, !binary] <- mu_eta
  Delta1[, !binary] <- partial_mu_eta(eta, family, 2L)
  for (j in seq.int(p)) {
    if (binary[[j]]) {
      eta0 <- eta - X[, j] * beta[[j]]
      f1 <- family[["mu.eta"]](eta0 + beta[[j]])
      Delta[, j] <- family[["linkinv"]](eta0 + beta[[j]]) - family[["linkinv"]](eta0)
      Delta1[, j] <- f1 - family[["mu.eta"]](eta0)
      J[, j] <- - colSums(PX * Delta1[, j]) / nt_full
      J[j, j] <- sum(f1) / nt_full + J[j, j]
      J[- j, j] <- colSums(X[, - j, drop = FALSE] * Delta1[, j]) / nt_full + J[- j, j]
      rm(eta0, f1)
    } else {
      Delta[, j] <- beta[[j]] * Delta[, j]
      Delta1[, j] <- beta[[j]] * Delta1[, j]
      J[, j] <- colSums(MX * Delta1[, j]) / nt_full
      J[j, j] <- sum(mu_eta) / nt_full + J[j, j]
    }
  }
  delta <- colSums(Delta) / nt_full
  Delta <- t(t(Delta) - delta) / nt_full
  rm(mu, mu_eta, PX)
  
  # Compute projection of \Psi
  Psi <- - Delta1 / w
  MPsi <- center_variables(Psi * sqrt(w), sqrt(w), Ti) / sqrt(w)
  PPsi <- Psi - MPsi
  rm(Delta1, Psi)
  
  # Compute analytical bias correction of average partial effects
  if (biascorr) {
    # Compute second-order partial derivatives
    Delta2 <- matrix(NA_real_, nt, p)
    Delta2[, !binary] <- partial_mu_eta(eta, family, 3L)
    for (j in seq.int(p)) {
      if (binary[[j]]) {
        eta0 <- eta - X[, j] * beta[[j]]
        Delta2[, j] <- partial_mu_eta(eta0 + beta[[j]], family, 2L) - 
          partial_mu_eta(eta0, family, 2L)
        rm(eta0)
      } else {
        Delta2[, j] <- beta[[j]] * Delta2[, j]
      }
    }
    
    # Compute corrected average partial effect using the bias correction of Fernández-Val (2009)
    B <- as.vector(group_sums_bias(Delta2 + PPsi * z, w, Ti)) / 2.0
    rm(Delta2)
    if (L > 0L) {
      B <- B - as.vector(group_sums_spectral(MPsi * w, v, w, L, Ti))
    }
    delta <- delta - B / nt
  }
  rm(eta, w, z, MPsi)
  
  # Compute standard errors
  WinvJ <- solve(object[["Hessian"]] / nt_full, J)
  Gamma <- (MX %*% WinvJ - PPsi) * v / nt_full
  V <- crossprod(Gamma)
  if (adj > 0.0) {
    # Simplify covariance if sampling assumptions are imposed
    if (sampling_fe == "independence") {
      V <- V + adj * group_sums_var(Delta, Ti)
    }
    
    # Add covariance in case of weak exogeneity
    if (weak_exo) {
      C <- group_sums_cov(Delta, Gamma, Ti)
      V <- V + adj * (C + t(C))
    }
  }
  
  # Add names
  names(delta) <- names(beta)
  dimnames(V) <- list(names(beta), names(beta))
  
  # Return result list
  structure(list(
    delta       = delta,
    vcov        = V,
    sampling_fe = sampling_fe,
    weak_exo    = weak_exo
    ), class = "bifeAPEs")
}


### Deprecated functions

#' @rdname get_APEs
#' @aliases get_APEs
#' @export
apeff_bife <- function(...) {
  .Deprecated("get_APEs")
}