R/NegBinomGLM.R

In gofreg: Bootstrap-Based Goodness-of-Fit Tests for Parametric Regression

##' @title Generalized linear model with negative binomial distribution
##' @description This class represents a generalized linear model with negative
##'   binomial distribution. It inherits from [GLM] and implements its functions
##'   that, for example, evaluate the conditional density and distribution
##'   functions.
##' @param x matrix of covariates, each row representing one sample
##' @param params model parameters to use (`list()` with tags beta and shape),
##'   defaults to the fitted parameter values
##'
##' @importFrom R6 R6Class
##'
##' @export
##'
##' @examples
##' # Use the built-in cars dataset
##' x <- datasets::cars$speed
##' y <- datasets::cars$dist
##' data <- dplyr::tibble(x=x, y=y)
##'
##' # Create an instance of a NegBinomGLM
##' model <- NegBinomGLM$new()
##'
##' # Fit a Negative Binomial GLM to the cars dataset
##' model$fit(data, params_init = list(beta=3, shape=2), inplace=TRUE)
##' params_opt <- model$get_params()
##'
##' # Plot the resulting regression function
##' plot(datasets::cars)
##' abline(a = 0, b = params_opt$beta)
##'
##' # Generate a sample for y for given x following the same distribution
##' x.new <- seq(min(x), max(x), by=2)
##' y.smpl <- model$sample_yx(x.new)
##' points(x.new, y.smpl, col="red")
##'
##' # Evaluate the conditional density, distribution, quantile and regression
##' # function at given values
##' model$f_yx(y.smpl, x.new)
##' model$F_yx(y.smpl, x.new)
##' model$F1_yx(y.smpl, x.new)
##' y.pred <- model$mean_yx(x.new)
##' points(x.new, y.pred, col="blue")
NegBinomGLM <- R6::R6Class(
  classname = "NegBinomGLM",
  inherit = GLM,
  public = list(
    #' @description Calculates the maximum likelihood estimator for the model
    #'   parameters based on given data.
    #'
    #' @param data tibble containing the data to fit the model to
    #' @param params_init initial value of the model parameters to use for the
    #'   optimization (defaults to the fitted parameter values)
    #' @param loglik `function(data, model, params)` defaults to [loglik_xy()]
    #' @param inplace `logical`; if `TRUE`, default model parameters are set
    #'   accordingly and parameter estimator is not returned
    #'
    #' @return MLE of the model parameters for the given data, same shape as
    #'   `params_init`
    #' @export
    fit = function(data, params_init = private$params, loglik = loglik_xy, inplace = FALSE) {
      checkmate::assert_names(names(data), must.include = c("x"))
      private$check_params(params_init, data$x)
      params_opt <- super$fit(data, params_init = unlist(params_init, use.names = FALSE), loglik = loglik)
      params_opt <- list(beta = params_opt[-length(params_opt)], shape = params_opt[length(params_opt)])
      if (inplace) {
        private$params <- params_opt
        invisible(self)
      } else {
        params_opt
      }
    },

    #' @description Evaluates the conditional density function.
    #'
    #' @param t value(s) at which the conditional density shall be evaluated
    #'
    #' @return value(s) of the conditional density function, same shape as `t`
    #' @export
    f_yx = function(t, x, params = private$params) {
      super$check_params(params)
      # when computing the MLE, params is a plain vector and needs to be reshaped
      if (checkmate::test_atomic_vector(params)) {
        checkmate::assert_atomic_vector(params, len = 1 + ncol(x))
        params <- list(beta = params[-length(params)], shape = params[length(params)])
      } else {
        private$check_params(params, x)
      }
      mean <- self$mean_yx(x, params)
      shape <- params$shape
      dnbinom(t, mu = mean, size = shape)
    },

    #' @description Evaluates the conditional distribution function.
    #'
    #' @param t value(s) at which the conditional distribution shall be
    #'   evaluated
    #'
    #' @return value(s) of the conditional distribution function,  same shape as
    #'   `t`
    #' @export
    F_yx = function(t, x, params = private$params) {
      super$check_params(params)
      # when computing the MLE, params is a plain vector and needs to be reshaped
      if (checkmate::test_atomic_vector(params)) {
        checkmate::assert_atomic_vector(params, len = 1 + ncol(x))
        params <- list(beta = params[-length(params)], shape = params[length(params)])
      } else {
        private$check_params(params, x)
      }
      mean <- self$mean_yx(x, params)
      shape <- params$shape
      pnbinom(t, mu = mean, size = shape)
    },

    #' @description Evaluates the conditional quantile function.
    #'
    #' @param t value(s) at which the conditional quantile function shall be
    #'   evaluated
    #'
    #' @return value(s) of the conditional quantile function, same shape as
    #'   `t`
    #' @export
    F1_yx = function(t, x, params = private$params) {
      private$check_params(params, x)
      mean <- self$mean_yx(x, params)
      shape <- params$shape
      qnbinom(t, mu = mean, size = shape)
    },

    #' @description Generates a new sample of response variables with the same
    #'   conditional distribution.
    #'
    #' @return vector of sampled response variables, same length as `nrow(x)`
    #' @export
    sample_yx = function(x, params = private$params) {
      private$check_params(params, x)
      mean <- self$mean_yx(x, params)
      shape <- params$shape
      rnbinom(length(mean), mu = mean, size = shape)
    }
  ),
  private = list(
    # @description Check that `params` have the correct form
    check_params = function(params, x) {
      super$check_params(params)
      checkmate::assert_list(params, len = 2)
      checkmate::assert_names(names(params), identical.to = c("beta", "shape"))
      checkmate::assert_vector(params$beta, len = ncol(as.matrix(x)))
    }
  )
)