R/methods.R

In logitr: Logit Models w/Preference & WTP Space Utility Parameterizations

#' Methods for logitr objects
#'
#' Miscellaneous methods for `logitr` class objects.
#'
#' @name miscmethods.logitr
#' @aliases logLik.logitr terms.logitr coef.logitr coef.summary.logitr
#' summary.logitr print.logitr print.summary.logitr
#' @param x is an object of class `logitr`.
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param digits the number of digits for printing, defaults to `3`.
#' @param width the width of the printing.
#' @param ... further arguments.
#'
#' @rdname miscmethods.logitr
#' @export
logLik.logitr <- function(object, ...) {
  return(structure(
    object$logLik,
    df    = object$n$pars,
    null  = sum(object$freq * log(object$freq / object$n$obs)),
    class = "logLik"
  ))
}

#' @rdname miscmethods.logitr
#' @export
terms.logitr <- function(x, ...) {
  return(stats::terms(x$formula))
}

#' @rdname miscmethods.logitr
#' @export
coef.logitr <- function(object, ...) {
    return(object$coefficients)
}

#' @rdname miscmethods.logitr
#' @method coef summary.logitr
#' @export
coef.summary.logitr <- function(object, ...) {
    return(object$coefTable)
}

#' @rdname miscmethods.logitr
#' @export
summary.logitr <- function (object, ...) {
    object$modelInfoTable <- getModelInfoTable(object)
    coefs <- stats::coef(object)
    standErr <- se(object)
    object$coefTable <- getCoefTable(coefs, standErr)
    object$statTable <- getStatTable(object)
    if (object$modelType == "mxl") {
        object$randParSummary <- getRandParSummary(object)
    }
    class(object) <- c("summary.logitr", "logitr")
    return(object)
}

#' @rdname miscmethods.logitr
#' @export
print.logitr <- function (
  x,
  digits = max(3, getOption("digits") - 2),
  width = getOption("width"),
  ...
) {
  cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
  modelType <- getModelType(x)
  modelSpace <- getModelSpace(x)
  cat("A", modelType, "model estimated in the", modelSpace, "space\n\n")
  cat("Exit Status: ", x$status, ", ", getExitMessage(x), "\n\n", sep = "")
  # Print which run was "best" if a multistart was used
  if (x$n$multiStarts > 1) {
    modelRun <- getModelRun(x)
    cat(
      "Results below are from run", modelRun, "multistart runs\n",
      "as it had the largest log-likelihood value\n")
    cat("\n")
  }
  # Print coefficients & log-likelihood
  if (!any(is.na(stats::coef(x)))) {
      cat("Coefficients:\n")
      print.default(
        format(x$coefficients, digits = digits), print.gap = 2, quote = FALSE)
      print(x$logLik)
  } else {
    cat("No coefficients\n")
  }
  cat("\n")
  invisible(x)
}

#' @rdname miscmethods.logitr
#' @method print summary.logitr
#' @export
print.summary.logitr <- function(
  x,
  digits = max(3, getOption("digits") - 2),
  width = getOption("width"),
  ...
) {
  cat("=================================================", "\n", sep = "")
  if (is.null(x$date)) {x$date <- "date missing"}
  if (is.null(x$version)) {x$date <- "version missing"}
  cat("\nModel estimated on:", x$date, "\n")
  cat("\nUsing logitr version:", x$version, "\n\n")
  cat("Call:\n")
  print(x$call)
  cat("\n")
  cat("Frequencies of alternatives:\n")
  print(prop.table(x$freq), digits = digits)
  # Print multistart summary
  if (!is.null(x$multistartSummary)) {
    cat("\n")
    cat("Summary Of Multistart Runs:\n")
    print(x$multistartSummary)
    cat("\n")
    cat("Use statusCodes() to view the meaning of each status code\n")
  }
  cat("\nExit Status: ", x$status, ", ", getExitMessage(x), "\n", sep = "")
  print(x$modelInfoTable)
  cat("\n")
  cat("Model Coefficients:", "\n")
  stats::printCoefmat(x$coefTable, digits = digits)
  print(x$statTable)
  if (x$modelType == "mxl") {
      cat("\n")
      cat("Summary of 10k Draws for Random Coefficients:", "\n")
      print(x$randParSummary)
  }
  invisible(x)
}

getModelInfoTable <- function(object) {
  modelType <- getModelType(object)
  modelSpace <- getModelSpace(object)
  modelRun <- getModelRun(object)
  modelTime <- convertTime(object$time)
  algorithm <- object$options$algorithm
  modelInfoTable <- data.frame(c(
    modelType, modelSpace, modelRun, object$iterations,
    modelTime, algorithm, object$weightsUsed
  ))
  colnames(modelInfoTable) <- ""
  row.names(modelInfoTable) <- c(
    "Model Type:", "Model Space:", "Model Run:", "Iterations:",
    "Elapsed Time:", "Algorithm:", "Weights Used?:"
  )
  panelID <- object$inputs$panelID
  if (!is.null(panelID)) {
    modelInfoTable <- rbind(modelInfoTable, panelID)
    row.names(modelInfoTable)[nrow(modelInfoTable)] <- "Panel ID:"
  }
  if (!is.null(object$n$clusters)) {
    if (object$n$clusters > 0) {
      modelInfoTable <- rbind(modelInfoTable, object$inputs$clusterID)
      row.names(modelInfoTable)[nrow(modelInfoTable)] <- "Cluster ID:"
    }
  }
  robust <- object$inputs$robust
  if (!is.null(robust)) {
    modelInfoTable <- rbind(modelInfoTable, robust)
    row.names(modelInfoTable)[nrow(modelInfoTable)] <- "Robust?"
  }
  return(modelInfoTable)
}

getCoefTable <- function(coefs, standErr) {
    z <- coefs / standErr
    p <- 2 * (1 - stats::pnorm(abs(z)))
    coefTable <- cbind(coefs, standErr, z, p)
    colnames(coefTable) <- c("Estimate", "Std. Error", "z-value", "Pr(>|z|)")
    return(as.data.frame(coefTable))
}

getStatTable <- function(object) {
  aic <- stats::AIC(object)
  bic <- round(log(object$n$obs) * object$n$pars - 2 * object$logLik, 4)
  mcR2 <- 1 - (object$logLik / object$nullLogLik)
  adjMcR2 <- 1 - ((object$logLik - object$n$pars) / object$nullLogLik)
  statTable <- data.frame(c(
    object$logLik, object$nullLogLik, aic, bic, mcR2, adjMcR2, object$n$obs
  ))
  colnames(statTable) <- ""
  row.names(statTable) <- c(
    "Log-Likelihood:", "Null Log-Likelihood:", "AIC:", "BIC:", "McFadden R2:",
    "Adj McFadden R2:" , "Number of Observations:")
  if (!is.null(object$n$clusters)) { # Added for backwards compatibility
    if (object$n$clusters > 0) {
      statTable <- rbind(statTable, object$n$clusters)
      row.names(statTable)[nrow(statTable)] <- "Number of Clusters"
    }
  }
  return(statTable)
}

getRandParSummary <- function(object) {
  parSetup <- object$parSetup
  parIDs <- object$parIDs
  n <- object$n
  n$draws <- 10^4
  standardDraws <- getStandardDraws(parIDs, n$draws, 'halton')
  betaDraws <- makeBetaDraws(
      stats::coef(object), parIDs, n, standardDraws, object$inputs$correlation)
  randParSummary <- apply(betaDraws, 2, summary)
  # Add names to summary
  nIDs <- parIDs$n
  lnIDs <- parIDs$ln
  cnIDs <- parIDs$cn
  distName <- rep("", length(parSetup))
  distName[nIDs] <- "normal"
  distName[lnIDs] <- "log-normal"
  distName[cnIDs] <- "zero-censored normal"
  summaryNames <- paste(names(parSetup), " (", distName, ")", sep = "")
  colnames(randParSummary) <- summaryNames
  randParSummary <- as.data.frame(t(randParSummary))
  # Set min and max values for unbounded distributions
  if (length(nIDs) > 0) {
      randParSummary[nIDs,]$Min. <- -Inf
      randParSummary[nIDs,]$Max. <- Inf
  }
  if (length(lnIDs) > 0) {
      randParSummary[lnIDs,]$Min. <- 0
      randParSummary[lnIDs,]$Max. <- Inf
  }
  if (length(cnIDs) > 0) {
      randParSummary[cnIDs,]$Max. <- Inf
  }
  # Add names and drop fixed pars
  randParSummary <- randParSummary[parIDs$r,]
  row.names(randParSummary) <- names(parIDs$r)
  return(randParSummary)
}

getModelType <- function(x) {
  return(ifelse(x$modelType == "mnl", "Multinomial Logit", "Mixed Logit"))
}

getModelSpace <- function(x) {
  return(ifelse(
    x$modelSpace == "pref", "Preference", "Willingness-to-Pay"))
}

getModelRun <- function(x) {
  return(paste(x$multistartNumber, "of", x$n$multiStarts))
}

getExitMessage <- function(x) {
  codes <- getStatusCodes()
  return(codes$message[which(codes$code == x$status)])
}

#' Extract standard errors
#'
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param ... further arguments.
#' @export
se <- function(object, ...) {
  UseMethod("se")
}

#' Extract standard errors
#'
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param ... further arguments.
#' @export
se.logitr <- function(object, ...) {
  return(sqrt(diag(stats::vcov(object))))
}

#' Calculate the variance-covariance matrix
#'
#' Returns the variance-covariance matrix of the main parameters of a fitted
#' model object.
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param ... further arguments.
#' @export
vcov.logitr <- function(object, ...) {
  if (!is.null(object$vcov)) {
      # vcov was already computed during model estimation
      return(object$vcov)
  }
  if (is.null(object$data$clusterID) | object$inputs$robust == FALSE) {
    return(getCovarianceNonRobust(object$hessian))
  }
  return(getCovarianceRobust(object))
}

getCovarianceNonRobust <- function(hessian) {
  covariance <- hessian*NA
  tryCatch(
    {
      covariance <- solve(-1*hessian)
    },
    error = function(e) {}
  )
  return(covariance)
}

getCovarianceRobust <- function(object) {
  numClusters <- object$n$clusters
  inputs <- object$inputs
  parSetup <- object$parSetup
  modelInputs <- list(
    logitFuncs = setLogitFunctions(object$modelSpace),
    evalFuncs  = setEvalFunctions(object$modelType, inputs$useAnalyticGrad),
    inputs     = inputs,
    modelType  = object$modelType,
    modelSpace = object$modelSpace,
    n          = object$n,
    parSetup   = parSetup,
    parIDs     = object$parIDs,
    panel      = !is.null(inputs$panelID),
    standardDraws = object$standardDraws,
    data_diff = makeDiffData(object$data, object$modelType)
  )
  clusterID <- modelInputs$data_diff$clusterID
  clusters <- unique(clusterID)
  pars <- stats::coef(object)
  gradMat <- matrix(NA, nrow = numClusters, ncol = length(pars))
  for (i in seq_len(length(clusters))) {
    indices <- which(clusterID == i)
    tempMI <- getClusterModelInputs(indices, modelInputs, i)
    gradMat[i, ] <- -1 * modelInputs$evalFuncs$negGradLL(pars, tempMI)
  }
  gradMeanMat <- repmat(matrix(colMeans(gradMat), nrow = 1), numClusters, 1)
  diffMat <- gradMat - gradMeanMat
  M <- t(diffMat) %*% diffMat
  M <- M * (numClusters / (numClusters - 1)) # small sample correction
  D <- getCovarianceNonRobust(object$hessian)
  if (any(is.na(D))) { return(D) } # If there are NAs the next line will error
  return(D %*% M %*% D)
}

getClusterModelInputs <- function (indices, mi, i) {
  X <- mi$data_diff$X[indices,]
  # Cast to matrix in cases where there is 1 independent variable
  if (length(indices) == 1) {
      X <- matrix(X, nrow = 1)
  }
  mi$data_diff$X <- X
  mi$data_diff$scalePar <- mi$data_diff$scalePar[indices]
  obsID <- mi$data_diff$obsID[indices]
  unique_obsID <- unique(obsID)
  mi$data_diff$obsID <- obsID
  panelID <- mi$data_diff$panelID[unique_obsID]
  weights <- mi$data_diff$weights[unique_obsID]
  if (!is.null(panelID)) {
    weights <- mi$data_diff$weights[i]
  }
  mi$data_diff$panelID <- panelID
  mi$data_diff$weights <- weights
  mi$data_diff$clusterID <- NULL
  if (isMxlModel(mi$parSetup)) {
    mi$partials <- makePartials(mi, mi$data_diff)
  }
  mi$n$rowX <- nrow(X)
  return(mi)
}

#' @rdname miscmethods.logitr
#' @export
print.logitr_wtp <- function (
  x,
  digits = max(3, getOption("digits") - 2),
  width = getOption("width"),
  ...
) {
  stats::printCoefmat(x, digits = digits)
}

#' Extract Model Fitted Values
#'
#' Returns fitted values from an object of class `logitr`.
#' @keywords logitr fitted fitted.values
#'
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param probs Predicted probabilities for an object of class `logitr` to use
#' in computing fitted values Defaults to `NULL`.
#' @param ... further arguments.
#'
#' @return A data frame of the `obsID` and the fitted values extracted from
#' `object`.
#' @export
#' @examples
#' library(logitr)
#'
#' # Estimate a preference space model
#' mnl_pref <- logitr(
#'   data    = yogurt,
#'   outcome = "choice",
#'   obsID   = "obsID",
#'   pars    = c("price", "feat", "brand")
#' )
#'
#' # Extract the fitted values from the model
#' fitted(mnl_pref)
fitted.logitr <- function(object, probs = NULL, ...) {
  if (is.null(probs)) {
    probs <- stats::predict(object, type = "prob")
  }
  outcome <- object$data$outcome
  fitted <- probs[which(outcome == 1),]
  names(fitted)[which(names(fitted) == 'predicted_prob')] <- "fitted_value"
  return(fitted)
}

#' Extract Model Residuals
#'
#' Returns model residuals from an object of class `logitr`.
#' @keywords logitr residuals resid
#'
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param fitted Fitted values for an object of class `logitr` to use in
#' computing residuals. Defaults to `NULL`.
#' @param ... further arguments.
#'
#' @return A data frame of the `obsID` and the residuals (response minus fitted
#' values) extracted from `object`.
#' @export
#' @examples
#' library(logitr)
#'
#' # Estimate a preference space model
#' mnl_pref <- logitr(
#'   data    = yogurt,
#'   outcome = "choice",
#'   obsID   = "obsID",
#'   pars    = c("price", "feat", "brand")
#' )
#'
#' # Extract the residuals from the model
#' residuals(mnl_pref)
residuals.logitr <- function(object, fitted = NULL, ...) {
  if (is.null(fitted)) {
    fitted <- stats::fitted(object)
  }
  reps <- table(object$data$obsID)
  residuals <- fitted[rep(seq_along(reps), reps),]
  resids <- object$data$outcome - residuals$fitted_value
  residuals$fitted_value <- NULL
  residuals$residual <- as.vector(resids)
  return(residuals)
}

#' Extract Model Confidence Interval
#'
#' Returns confidence intervals from an object of class `logitr`.
#' @keywords logitr confint
#'
#' @param object is an object of class `logitr` (a model estimated using
#' the 'logitr()` function).
#' @param parm A specification of which parameters are to be given confidence
#' intervals, either a vector of numbers or a vector of names.
#' If missing, all parameters are considered.
#' @param level The confidence level required.
#' @param ... further arguments.
#'
#' @return A data frame of the confidence intervals of model coefficients.
#' @export
#' @examples
#' library(logitr)
#'
#' # Estimate a preference space model
#' mnl_pref <- logitr(
#'   data    = yogurt,
#'   outcome = "choice",
#'   obsID   = "obsID",
#'   pars    = c("price", "feat", "brand")
#' )
#'
#' # Compute a confidence interval
#' confint(mnl_pref)
confint.logitr <- function(object, parm, level = 0.95, ...) {
    draws <- getUncertaintyDraws(object, numDraws = 10^4)
    lower <- (1 - level)/2
    upper <- 1 - lower
    df <- data.frame(
        lower = apply(draws, 2, function(x) fquantile(x, lower, na.rm = TRUE)),
        upper = apply(draws, 2, function(x) fquantile(x, upper, na.rm = TRUE))
    )
    names(df) <- c(paste(lower*100, "%"), paste(upper*100, "%"))
    return(df)
}

#' Construct Design Matrices
#'
#' Creates a design (or model) matrix, e.g., by expanding factors to a set of
#' dummy variables (depending on the contrasts) and expanding interactions
#' similarly.
#' @keywords logitr model.matrix
#'
#' @param object an object of an appropriate class. For the default method,
#' a model `formula` or a `terms` object.
#' @param ... further arguments.
#'
#' @return A design matrix
#' @export
#' @examples
#' library(logitr)
#'
#' # Estimate a preference space model
#' mnl_pref <- logitr(
#'   data    = yogurt,
#'   outcome = "choice",
#'   obsID   = "obsID",
#'   pars    = c("price", "feat", "brand")
#' )
#'
#' # Get the model.matrix design matrix
#' model.matrix(mnl_pref)
model.matrix.logitr <- function(object, ...) {
    return(object$data$X)
}

#' Extracting the Model Frame from a Formula or Fit
#'
#' Returns a data.frame with the variables needed to use formula and
#' any `...` arguments.
#' @keywords logitr model.frame
#'
#' @param formula a model `formula` or `terms` object or an R object.
#' @param ... further arguments.
#'
#' @return A data.frame with the variables needed to use formula and
#' any `...` arguments.
#' @export
#' @examples
#' library(logitr)
#'
#' # Estimate a preference space model
#' mnl_pref <- logitr(
#'   data    = yogurt,
#'   outcome = "choice",
#'   obsID   = "obsID",
#'   pars    = c("price", "feat", "brand")
#' )
#'
#' # Get the model.frame data frame
#' model.frame(mnl_pref)
model.frame.logitr <- function(formula, ...) {
    stats::model.frame.default(
        formula$formula,
        data = eval(formula$call$data, envir = parent.frame())
    )
}