R/PLSstandarderrors.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Extract Empirical Estimating Functions
#'
#' @description Generic function for extracting the empirical estimating functions of a fitted model.
#' @param x a fitted model object.
#' @param ... arguments passed to methods.
#' @method estfun gpls
#' @importFrom sandwich estfun
#' @return A matrix containing the empirical estimating functions. Typically, this should be an n x k matrix corresponding to n observations and k parameters. The columns should be named as in coef or terms, respectively.
#' The estimating function (or score function) for a model is the derivative of the objective function with respect to the parameter vector. The empirical estimating functions is the evaluation of the estimating function at the observed data (n observations) and the estimated parameters (of dimension k).
#' @export
#' @keywords internal
estfun.gpls <- function (x, ...) {
  xmat <- x$model.matrix
  wres <- residuals.gpls(x, type = "working") * diag(x$irwls.wts)
  dispersion <- if (x$family == "binomial" || x$family == "poisson" || x$family == "negative.binomial")
    1
  else sum(wres^2, na.rm = TRUE)/sum(diag(x$irwls.wts))
  rval <- wres * xmat/dispersion
  attr(rval, "assign") <- NULL
  attr(rval, "contrasts") <- NULL
  return(rval)
}

#' Heteroskedasticity-Consistent Covariance Matrix Estimation
#'
#' @param x a fitted model object.
#' @param type a character string specifying the estimation type. For details see below.
#' @param omega a vector or a function depending on the arguments residuals (the working residuals of the model), diaghat (the diagonal of the corresponding hat matrix) and df (the residual degrees of freedom). For details see below.
#' @param ... arguments passed to sandwich (in vcovHC) and estfun (in meatHC), respectively.
#' @return A matrix containing the covariance matrix estimate.
#' @export meatHC.gpls
#' @keywords internal
meatHC.gpls <- function (x, type = c("HC4m", "const", "HC", "HC0", "HC1", "HC2", "HC3", "HC4", "HC5"), omega = NULL, ...){
  if (is.list(x) && !is.null(x$na.action))
    class(x$na.action) <- "omit"
  X <- x$model.matrix
  if (any(alias <- is.na(coef(x))))
    X <- X[, !alias, drop = FALSE]
  attr(X, "assign") <- NULL
  n <- NROW(X)
  diaghat <- x$hat
  df <- n - NCOL(X)
  ef <- estfun.gpls(x, ...)
  res <- rowMeans(ef/X, na.rm = TRUE)
  res[apply(abs(ef) < .Machine$double.eps, 1L, all)] <- 0
  if (is.null(omega)) {
    type <- match.arg(type)
    if (type == "HC")
      type <- "HC0"
    switch(type, const = {
      stop("Ordinary standard errors are not available for gpls models.")
    }, HC0 = {
      omega <- function(residuals, diaghat, df) residuals^2
    }, HC1 = {
      omega <- function(residuals, diaghat, df) residuals^2 * length(residuals)/df
    }, HC2 = {
      omega <- function(residuals, diaghat, df) residuals^2/(1 -diaghat)
    }, HC3 = {
      omega <- function(residuals, diaghat, df) residuals^2/(1 -diaghat)^2
    }, HC4 = {
      omega <- function(residuals, diaghat, df) {
        n <- length(residuals)
        p <- as.integer(round(sum(diaghat), digits = 0))
        delta <- pmin(4, n * diaghat/p)
        residuals^2/(1 - diaghat)^delta
      }
    }, HC4m = {
      omega <- function(residuals, diaghat, df) {
        gamma <- c(1, 1.5)
        n <- length(residuals)
        p <- as.integer(round(sum(diaghat), digits = 0))
        delta <- pmin(gamma[1], n * diaghat/p) + pmin(gamma[2], n * diaghat/p)
        residuals^2/(1 - diaghat)^delta
      }
    }, HC5 = {
      omega <- function(residuals, diaghat, df) {
        k <- 0.7
        n <- length(residuals)
        p <- as.integer(round(sum(diaghat), digits = 0))
        delta <- pmin(n * diaghat/p, pmax(4, n * k * max(diaghat)/p))
        residuals^2/sqrt((1 - diaghat)^delta)
      }
    })
  }
  if (is.function(omega))
    omega <- omega(res, diaghat, df)
  rval <- sqrt(omega) * X
  rval <- crossprod(rval)/n
  return(rval)
}

#' Bread for Sandwiches
#'
#' @description Generic function for extracting an estimator for the bread of sandwiches.
#' @param x a fitted model object.
#' @param ... arguments passed to methods.
#' @method bread gpls
#' @return A matrix containing an estimator for the expectation of the negative derivative of the estimating functions, usually the Hessian. Typically, this should be an k x k matrix corresponding to k parameters. The rows and columns should be named as in coef or terms, respectively.
#' The default method tries to extract vcov and nobs and simply computes their product.
#' @importFrom sandwich bread
#' @export
#' @keywords internal
bread.gpls <- function (x, ...) {
  if (!is.null(x$na.action))
    class(x$na.action) <- "omit"
  xmat <- x$model.matrix
  wres <- residuals.gpls(x, type = "working") * diag(x$irwls.wts)
  dispersion <- if (x$family == "binomial" || x$family == "poisson" || x$family == "negative.binomial")
    1
  else sum(wres^2, na.rm = TRUE)/sum(diag(x$irwls.wts))
  return(XtXinv(xmat, W = x$irwls.wts) * nrow(xmat) * dispersion)
}


#' Heteroskedasticity-Consistent Covariance Matrix Estimation
#'
#' @param x a fitted model object.
#' @param type a character string specifying the estimation type. For details see below.
#' @param omega a vector or a function depending on the arguments residuals (the working residuals of the model), diaghat (the diagonal of the corresponding hat matrix) and df (the residual degrees of freedom). For details see below.
#' @param sandwich 	logical. Should the sandwich estimator be computed? If set to FALSE only the meat matrix is returned.
#' @param ... arguments passed to sandwich (in vcovHC) and estfun (in meatHC), respectively.
#' @method vcovHC gpls
#' @return A matrix containing the covariance matrix estimate.
#' @importFrom sandwich vcovHC
#' @export
#' @keywords internal
vcovHC.gpls <- function (x, type = c("HC4m", "const", "HC", "HC0", "HC1", "HC2", "HC3", "HC4", "HC5"), omega = NULL, sandwich = TRUE, ...){
  type <- match.arg(type)
  if (type == "const"){
    stop("Ordinary standard errors are not available for gpls models.")
  }
  rval <- meatHC.gpls(x, type = type, omega = omega)
  bread <- bread.gpls(x)
  if (sandwich)

    if (x$family == "binomial" || x$family == "poisson" || x$family == "negative.binomial")
      rval <- (bread %*% rval %*% bread) / nrow(x$model.matrix)
    else
      rval <- (bread %*% rval %*% bread) / nrow(x$model.matrix)
  return(rval)
}