R/nter.R

In BayerSe/trexreg: Truncated Expectation Regression

#' Non-Tuncated Expectation Regression
#'
#' TBA
#'
#' @export
nter <- function(...) UseMethod('nter')

nter.default <- function(...) {
  stop('This method should not be called directly.
       Use the formula interface instead.')
}

nter.formula <- function(formula, data = parent.frame(), ...) {
  chkDots(...)
  if(is.matrix(data)) {
    data <- as.data.frame(data)
  }

  # Preparation steps
  formula <- Formula::Formula(formula)
  formula <- nter.check_formula(formula)
  mf <- stats::model.frame(formula = formula, data = data, na.action = stats::na.pass)
  terms <- attr(mf, 'terms')

  # Extract data
  xm <- stats::model.matrix(formula, mf, rhs=1)
  y <- stats::model.response(mf)

  # Fit model
  fit <- nter.fit(xm = xm, y = y)

  # Store multiple objects
  fit$formula <- formula
  fit$call <- match.call()
  fit$terms <- terms
  fit$xlevels <- stats::.getXlevels(terms, mf)
  fit$data <- data
  fit$y <- y
  fit$x <- xm
  fit$lq <- 0
  fit$uq <- 1

  fit
}

#' @export
nter.fit <- function(xm, y, ...) {
  chkDots(...)
  xm <- check.x(xm)
  y <- check.y(xm, y)
  fit <- do.nter.fit(xm = xm, y = y)
  class(fit) <- c('ter', 'nter')

  fit
}

#' @export
do.nter.fit <- function(xm, y) {
  # QR-decomposition of x
  qx <- qr(xm)
  if(qx$rank < ncol(xm)) {
    stop("'x' is singular (it has ", ncol(xm),
         " columns but its rank is ", qx$rank, ")")
  }

  # Estimate the parameters: compute (x'x)^(-1) x'y
  coef <- solve.qr(qx, y)
  stopifnot(!anyNA(coef))
  colnames(coef) <- colnames(y)

  # Store the fitted values
  fitted.values <- qr.fitted(qx, y)
  names(fitted.values) <- rownames(xm)

  # Store the loss
  loss <- mean((fitted.values - y)^2)

  # Retun results
  list(
    coefficients  = list(
      coef_mean=coef,
      coef_lq=NA,
      coef_uq=NA
    ),
    fitted.values = fitted.values,
    loss          = loss
  )
}

vcov.nter <- function(object, ...) {
  # Original White Covariance Estimator
  chkDots(...)
  y <- object$y
  x <- object$x
  e <- object$y - object$fitted.values
  coef <- object$coefficients$coef_mean
  xxi <- solve(crossprod(x))
  xd <- x * replicate(ncol(x), e)
  vcov <- xxi %*% crossprod(xd) %*% xxi

  colnames(vcov) <- rownames(vcov) <- colnames(x)
  vcov
}