R/gamlss_functions.R

In abnormally-distributed/abdisttools: Random Functions for Personal Use

#' Calculate model weights for AIC, BIC, DIC, etc.
#'
#' Takes a vector of information criterion and returns the model probabilities
#'
#' @param AIC the information criterion scores
#' @export
#' @examples
#' akaike_weights()
#'
#'
akaike_weights = function(AIC) {
  delta = AIC - min(AIC)
  round(exp(-.5*delta)/sum(exp(-.5*delta)), digits=4)
}

#' Generalized Likelihood Ratio Test for two nested GAMLSS models
#'
#' @param m1 model 1
#' @param m2 model 2
#' @export
#' @examples
#' GBIC()
GLR.test = function(m1, m2){
  cat("Generalized Likelihood Ratio Test for nested GAMLSS models.\n")
  cat(paste0("GLRT: (log) ", signif(m1$G.deviance - m2$G.deviance, 3 ) , "\n"))
  cat(paste0("GLRT: ", signif(exp(m1$G.deviance - m2$G.deviance),5) , "\n"))
  .5 * pchisq(m1$G.deviance - m2$G.deviance,1,lower.tail=FALSE)
}

#' BIC (SBC) function for gamlss objects analagous to gamlss's GAIC function.
#'
#' BIC (SBC) function for gamlss objects analagous to gamlss's GAIC function.
#'
#' @param object the gamlss model or models
#' @param c should small sample correction be applied?
#' @export
#' @examples
#' GBIC()

GBIC = function (object, ... , c = FALSE)
{


  if (length(list(...))) {
    object <- list(object, ...)
    isgamlss <- unlist(lapply(object, is.gamlss))
    if (!any(isgamlss))
      stop("some of the objects are not gamlss")
    df <- as.numeric(lapply(object, function(x) x$df.fit))
    N <- as.numeric(lapply(object, function(x) x$N))
    Cor <- if (c == TRUE)
      (2 * df * (df + 1))/(N - df - 1)
    else rep(0, length(object))
    BIC <- as.numeric(lapply(object, function(x) x$G.dev +
                               x$df.fit * log(length(fitted(x))))) + Cor
    val <- as.data.frame(cbind(df, BIC))
    Call <- match.call()
    Call$k <- NULL
    Call$c <- NULL
    row.names(val) <- as.character(Call[-1])
    val <- val[order(BIC), ]
    val
  }
  else {
    k = log(length(fitted(object)))
    val <- if (is.gamlss(object))
      object$G.dev + object$df.fit * k
    else stop(paste("this is not a gamlss object"))
    if ((k == 2) && (c == TRUE))
      val <- val + (2 * object$df.fit * (object$df.fit +
                                           1))/(object$N - object$df.fit - 1)
    val
  }
}


#' Load essential gam and gamlss packages at once
#'
#' Loads gamlss, mgcv, mgcViz, gamlss.mx, gamlss.cens, gamlss.add, gamlss.dist,
#' gamlss.tr, and gamlss.nl
#' @param load load=TRUE
#' @export
#' @examples
#' load.gamlss()

load.gamlss = function(load=TRUE){
  require("gamlss")
  require("mgcv")
  require("mgcViz")
  require("gamlss.mx")
  require("gamlss.util")
  require("gamlss.cens")
  require("gamlss.add")
  require("gamlss.dist")
  require("gamlss.tr")
  require("gamlss.nl")
}

#' this function provides robust covariance CI for GAMLSS models
#' 
#' this function provides robust sandwich estimator CIs for a GAMLSS model along
#' with exact MCMC p-values and parameter estimates and the standard errors. 
#' 
#' @param object the fitted gamlss model
#' @param parm parm
#' @param level the confidence level
#' @param what one of "mu", "sigma", "nu", "tau". Defaults to "mu".
#' @param parameter same as what
#' @param robust whether robust standard errors should be used. Defaults to TRUE. 
#' @export
#' @examples
#' gamlss.confint()
#' 
gamlss.confint = function (object, level = 0.95) 
{

  Sim_Post = function (n = 10000, estimates, vcov, tol = 1e-06) 
  {
    set.seed(666)
    mu = estimates
    Sigma = vcov
    p <- length(mu)
    if (!all(dim(Sigma) == c(p, p))) 
      stop("incompatible arguments")
    Sigma = fBasics::makePositiveDefinite(Sigma)
    eS <- eigen(Sigma, symmetric = TRUE)
    ev <- eS$values
    X <- matrix(abdisttools:::rstudent_t(p * n , df = 20), n)
    X <- base::scale(X, TRUE, FALSE)
    X <- X %*% svd(X, nu = 0)$v
    X <- base::scale(X, FALSE, TRUE)
    X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)
    nm <- names(mu)
    if (is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
      nm <- dn[[1L]]
    dimnames(X) <- list(nm, NULL)
    if (n == 1) 
      drop(X)
    else t(X)
  }
  
  names = paste0(broom::tidy(object)$parameter,"_",broom::tidy(object)$term)
  VCOV = rvcov_ab(object, type="vcov")
  post = Sim_Post(500000, as.vector(abdisttools:::vcov_ab_gamlss(object, type = "coef")), vcov=VCOV)
  colnames(post) = names
  post = post %>% as.data.frame() %>% mutate_if(is.numeric, round, 4) %>% as.matrix()

  summary = post_summary(post, estimate.method = "median", cred.method = "ETI", cred.level = level, pvals = TRUE) %>%
    mutate_if(is.numeric, round, 6)
  colnames(summary) = c("term", "estimate", "std.error","conf.low", "conf.high", "p.value")
  
  summary = summary %>% 
    transmute(term=term, estimate=estimate, std.error=std.error, conf.low = conf.low, conf.high = conf.high
              , t.statistic = estimate/std.error, p.value = p.value)
  
  return(summary)
}


#' Robust Variance-Covariance matrix of the parameters from a fitted GAMLSS model
#'
#' The function rvcov() is design for providing robust standard errors for the 
#' parameters estimates of a GAMLSS fitted model. The same result can be achieved
#' by using vcov(fitted_model,robust=TRUE). The function get.() gets the 
#' K matrix (see details below). This is a slightly modified version of the
#' original gamlss function that performs better when the variance-covariance
#' matrix is not positive definite by using the Moore-Penrose Pseudoinverse when
#' the matrix is problematic. This allows robust standard errors (the sandwich estimator
#' of the standard errors) to be obtained when they would otherwise be undefined without
#' resorting to the less accurate default standard errors. 
#'
#' @inheritParams gamlss::rvcov
#' @examples
#' rvcov(g)
#' @export
#' 

rvcov_ab = 
  function (object, type = c("vcov", "cor", "se", "coef", "all"), 
            hessian.fun = "PB") 
  {
    type <- match.arg(type)
    hessian.fun <- match.arg(hessian.fun)
    if (!("gamlss" %in% class(object))) 
      stop("the null model is not a gamlss model")
    V <- abdisttools:::vcov_ab_gamlss(object, hessian.fun = hessian.fun)
    K <- get.K(object)
    varCov <- V %*% K %*% V
    se <- sqrt(diag(varCov))
    corr <- cov2cor(varCov)
    switch(type, vcov = varCov, cor = corr, se = se, all = list(se = se, 
                                                                vcov = varCov, cor = corr))
  }


#' @export
vcov_ab_gamlss = function (object, type = c("vcov", "cor", "se", "coef", "all"), 
                        robust = TRUE, hessian.fun = c("R", "PB"), ...) 
{
  HessianPB <- function(pars, fun, ..., .relStep = (.Machine$double.eps)^(1/3), 
                        minAbsPar = 0) {
    pars <- as.numeric(pars)
    npar <- length(pars)
    incr <- ifelse(abs(pars) <= minAbsPar, minAbsPar * .relStep, 
                   abs(pars) * .relStep)
    baseInd <- diag(npar)
    frac <- c(1, incr, incr^2)
    cols <- list(0, baseInd, -baseInd)
    for (i in seq_along(pars)[-npar]) {
      cols <- c(cols, list(baseInd[, i] + baseInd[, -(1:i)]))
      frac <- c(frac, incr[i] * incr[-(1:i)])
    }
    indMat <- do.call("cbind", cols)
    shifted <- pars + incr * indMat
    indMat <- t(indMat)
    Xcols <- list(1, indMat, indMat^2)
    for (i in seq_along(pars)[-npar]) {
      Xcols <- c(Xcols, list(indMat[, i] * indMat[, -(1:i)]))
    }
    coefs <-solve(do.call("cbind", Xcols), apply(shifted, 
                                                                  2, fun, ...))/frac
    Hess <- diag(coefs[1 + npar + seq_along(pars)], ncol = npar)
    Hess[row(Hess) > col(Hess)] <- coefs[-(1:(1 + 2 * npar))]
    list(mean = coefs[1], gradient = coefs[1 + seq_along(pars)], 
         Hessian = (Hess + t(Hess)))
  }
  type <- match.arg(type)
  hessian.fun <- match.arg(hessian.fun)
  if (!is.gamlss(object)) 
    stop(paste("This is not an gamlss object", "\n", ""))
  coefBeta <- list()
  for (i in object$par) {
    if (length(eval(parse(text = paste(paste("object$", i, 
                                             sep = ""), ".fix==TRUE", sep = "")))) != 0) {
      ff <- eval(parse(text = paste(paste(paste(object$family[1], 
                                                "()$", sep = ""), i, sep = ""), ".linkfun", sep = "")))
      fixvalue <- ff(fitted(object, i)[1])
      names(fixvalue) <- paste("fixed", i, sep = " ")
      coefBeta <- c(coefBeta, fixvalue)
    }
    else {
      if (!is.null(unlist(attr(terms(formula(object, i), 
                                     specials = .gamlss.sm.list), "specials")))) 
        paste("")
      nonNAcoef <- !is.na(coef(object, i))
      coefBeta <- c(coefBeta, coef(object, i)[nonNAcoef])
    }
  }
  betaCoef <- unlist(coefBeta)
  like.fun <- gen.likelihood(object)
  hess <- if (hessian.fun == "R") 
    optimHess(betaCoef, like.fun)
  else HessianPB(betaCoef, like.fun)$Hessian
  varCov <- try(solve(hess), silent = TRUE)
  if (any(class(varCov) %in% "try-error") || any(diag(varCov) < 
                                                 0)) {
    varCov <- try(solve(HessianPB(betaCoef, like.fun)$Hessian), 
                  silent = TRUE)
    if (any(class(varCov) %in% "try-error")) 
      varCov <- try(corpcor::pseudoinverse(HessianPB(betaCoef, like.fun)$Hessian), 
                    silent = TRUE)
  }
  rownames(varCov) <- colnames(varCov) <- rownames(hess)
  se <- sqrt(diag(varCov))
  corr <- cov2cor(varCov)
  coefBeta <- unlist(coefBeta)
  if (robust) {
    K <- get.K(object)
    varCov <- varCov %*% K %*% varCov
    se <- sqrt(diag(varCov))
    corr <- cov2cor(varCov)
  }
  switch(type, vcov = varCov, cor = corr, se = se, coef = coefBeta, 
         all = list(coef = coefBeta, se = se, vcov = varCov, cor = corr))
}


#' Simulate multiple posterior distributions with gaussian approximation
#'
#' Simulate the posterior distributions of a model fit with gaussian approximation 
#' by supplying a vector of coefficients and their associated variance-covariance
#' matrix. 
#'
#' @param n number of samples. Defaults to 10000.
#' @param estimates the vector of coefficients/estimates
#' @param vcov the variance-covariance matrix
#' @param tol the tolerance
#' @export
#' @examples
#' sim_post(10000, coefs, rvcov)
#'
#'
sim_post = function (n = 10000, estimates, vcov, tol = 1e-06) 
{
  mu = estimates
  Sigma = vcov
  p <- length(mu)
  if (!all(dim(Sigma) == c(p, p))) 
    stop("incompatible arguments")
  Sigma = fBasics::makePositiveDefinite(Sigma)
  eS <- eigen(Sigma, symmetric = TRUE)
  ev <- eS$values
  X <- matrix(rnorm(p * n), n)
  X <- base::scale(X, TRUE, FALSE)
  X <- X %*% svd(X, nu = 0)$v
  X <- base::scale(X, FALSE, TRUE)
  X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% 
    t(X)
  nm <- names(mu)
  if (is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
    nm <- dn[[1L]]
  dimnames(X) <- list(nm, NULL)
  if (n == 1) 
    drop(X)
  else t(X)
}