R/simulCoef.R

In cirls: Constrained Iteratively Reweighted Least Squares

################################################################################
#
#    Simulate coefficients from truncated multivariate normal
#
################################################################################

#' Simulate coefficients, calculate Confidence Intervals and Variance-Covariance Matrix for a `cirls` object.
#'
#' @description
#' Simulates coefficients for a fitted `cirls` object. `confint` and `vcov` compute confidence intervals and the Variance-Covariance matrix for coefficients from a fitted `cirls` object. These methods supersede the default methods for `cirls` objects.
#'
#' @param object A fitted `cirls` object.
#' @param nsim The number of simulations to perform.
#' @param seed Either NULL or an integer that will be used in a call to [set.seed()] before simulating the coefficients.
#' @param complete If FALSE, doesn't return inference for undetermined coefficients in case of an over-determined model.
#' @param parm A specification of which parameters to compute the confidence intervals for. Either a vector of numbers or a vector of names. If missing, all parameters are considered.
#' @param level The confidence level required.
#' @param constrained A logical switch indicating Whether to simulate from the constrained (the default) or unconstrained coefficients distribution.
#' @param trunc If set to `FALSE` the unmodified GLM variance-covariance computed within [summary.glm()] is returned.
#' @param ... Further arguments passed to or from other methods. For `vcov` and `confint` can be used to provide a `seed` for the internal coefficient simulation.
#'
#' @details
#' `confint` and `vcov` are custom methods for [cirls][cirls.fit()] objects to supersede the default methods used for [glm][stats::glm()] objects. Internally, they both call `simulCoef` to generate coefficient vectors from a Truncated Multivariate Normal Distribution using the [TruncatedNormal::rtmvnorm()] function. This distribution accounts for truncation by constraints, ensuring all coefficients are feasible with respect to the constraint matrix. `simulCoef` typically doesn't need to be used directly for confidence intervals and variance-covariance matrices, but it can be used to check other summaries of the coefficients distribution.
#'
#' These methods only work when `Cmat` is of full row rank, i.e. if there are less constraints than variables in `object`.
#'
#' @note
#' By default, the Variance-Covariance matrix generated by `vcov` is different than the one returned by `summary(obj)$cov.scaled`. The former accounts for the reduction in degrees of freedom resulting from the constraints, while the latter is the unconstrained GLM Variance-Covariance. Note that the unconstrained one can be obtained from `vcov` by setting `constrained = FALSE`.
#'
#' @returns
#' For `simulCoef`, a matrix with `nsim` rows containing simulated coefficients.
#'
#' For `confint`, a two-column matrix with columns giving lower and upper confidence limits for each parameter.
#'
#' For `vcov`, a matrix of the estimated covariances between the parameter estimates of the model.
#'
#' @references
#' Geweke, J.F., 1996. Bayesian Inference for Linear Models Subject to Linear Inequality Constraints, in: Lee, J.C., Johnson, W.O., Zellner, A. (Eds.), Modelling and Prediction Honoring Seymour Geisser. *Springer, New York, NY*, pp. 248–263. \doi{10.1007/978-1-4612-2414-3_15}
#'
#' Botev, Z.I., 2017, The normal law under linear restrictions: simulation and estimation via minimax tilting, *Journal of the Royal Statistical Society, Series B*, **79** (**1**), pp. 1–24.
#'
#' @seealso [rtmvnorm][TruncatedNormal::tmvnorm()] for the underlying routine to simulate from a TMVN. [checkCmat()] to check if the contraint matrix can be reduced.
#'
#' @example man/examples/new_methods.R
#'
#' @order 1
#' @export
simulCoef <- function(object, nsim = 1, seed = NULL, complete = TRUE,
  constrained = TRUE)
{

  # Extract unconstrained coefficients
  ufit <- uncons(object)
  ubeta <- stats::coef(ufit, complete = FALSE)
  aliased <- stats::summary.glm(ufit)$aliased
  uvcov <- stats::.vcov.aliased(aliased, stats::summary.glm(ufit)$cov.scaled,
    complete = FALSE)

  # Check uvcov exists
  if (any(is.na(uvcov))){
    warning("Impossible to perform inference: unconstrained vcov matrix undefined. Returning NAs")
    return(matrix(NA, nsim, ifelse(complete, length(aliased), sum(aliased)),
      dimnames = list(NULL, names(aliased))))
  }

  # Set seed if necessary
  if (!is.null(seed)){
    R.seed <- get(".Random.seed", envir = .GlobalEnv)
    set.seed(seed)
    on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
  }

  # Extract contraints
  Cmat <- object$Cmat
  lb <- object$lb
  ub <- object$ub

  # If Cmat is empty, switch of the constrained simulation
  if (NROW(Cmat) == 0) constrained <- FALSE

  #----- Extract constraints and transform to "square" domain then simulate
  if (constrained){

    # Remove constraints affected by aliased coefficients
    Cmat <- Cmat[,!aliased, drop = F]
    keep <- rowSums(Cmat != 0)
    Cmat <- Cmat[as.logical(keep),, drop = F]
    lb <- lb[as.logical(keep)]
    ub <- ub[as.logical(keep)]

    # Check constraint matrix
    rowrk <- qr(t(Cmat))$rank
    if (nrow(Cmat) > rowrk){
      warning("Impossible to perform inference: constraint matrix not of full row rank. Returning NAs")
      return(matrix(NA, nsim, ifelse(complete, length(aliased), sum(aliased)),
        dimnames = list(NULL, names(aliased))))
    }

    # To allow back transformation, we "augment" the constraint matrix with
    #   its row null space (Tallis 1965)
    Hmat <- nullspace(t(Cmat))
    Bmat <- rbind(Cmat, t(Hmat))

    # Transform parameters
    simvcov <- Bmat %*% uvcov %*% t(Bmat)
    simbeta <- drop(Bmat %*% ubeta)

    # Expand bounds for simulation
    lowervec <- c(lb, rep(-Inf, ncol(Hmat)))
    uppervec <- c(ub, rep(Inf, ncol(Hmat)))

    # Initiate matrix with zeros for equality constraints
    truncres <- matrix(0, nrow = nsim, ncol = nrow(Bmat))
    eqind <- lowervec == uppervec

    # Simulate from truncated MVN
    truncres[,!eqind] <- suppressWarnings(TruncatedNormal::rtmvnorm(n = nsim,
      mu = simbeta, sigma = simvcov, lb = lowervec, ub = uppervec,
      check = FALSE))

    # Backtransform simulations
    simu <- t(solve(Bmat) %*% t(truncres))
  } else {

    #----- If not, simulate without bounds
    simu <- TruncatedNormal::rtmvnorm(n = nsim, mu = ubeta, sigma = uvcov,
      check = FALSE)
  }

  #----- Return, including NAs if complete == TRUE

  if (complete) {
    outsimu <- matrix(NA, nrow = nsim, ncol = length(aliased),
      dimnames = list(NULL, names(aliased)))
    outsimu[, which(!aliased)] <- simu
  } else {
    outsimu <- simu
    colnames(outsimu) <- names(aliased[!aliased])
  }

  # Export
  outsimu
}