R/inv_ci.R

In baruuum/btoolbox: random toolbox

#' Calculate Confidence Intervals for Predictions from Regression Models
#' 
#' Calculates confidence intervals of predicted outcomes by applying inverse-link function of GLMs to endpoints of confidence intervals of linear predictor
#' 
#' @param m regression model
#' @param newdata dataset to make predictions
#' @param Sigma variance-covariance matrix of estimator; if NULL, vcov(m) is used
#' @param level confidence level; defaults to .95
#' @param m_class model class; either glm or lm
#' @return Confidence interval of predicted response calculated at values specified in \code{newdata}.
#' @export
.inv_ci_pred = function(
    m = NULL,                  
    newdata = NULL,
    Sigma = NULL,
    level = .95,
    m_class = c("glm", "lm")) 
{

    # get class
    m_class = match.arg(m_class)
    
    message(
        paste0("Calculating ", 100 * level, 
               "% confidence intervals by transforming endpoints of confidence intervals of linear predictor ...")
    )

    p_type = if (m_class == "lm") "response" else "link"
    z = abs(qnorm(0.5 * (1.0 - level)))

    if (is.null(Sigma)) {
        
        xb = predict(m, newdata, type = p_type, se = TRUE)
        
        link_ci  = cbind(
            xb$fit - z * xb$se.fit,
            xb$fit + z * xb$se.fit
        )
        
    } else {
        
        X = create_mm(m, newdata)
        fit = X %*% coef(m)
        se = apply(X, 1L, function(w) {
            sqrt(w %*% Sigma %*% w)
        })
        
        link_ci = cbind(
            fit - z * se,
            fit + z * se
        )
        
    }
    
    if (m_class == "lm") {
        
        colnames(link_ci) = c("lower", "upper")
        return(link_ci)
        
    } 
    
    # swap upper/lower if inverse-link is used
    if (m$family$link == "inverse")
        link_ci = link_ci[, c(2L, 1L)]

    colnames(link_ci) = c("lower", "upper")
    
    return(m$family$linkinv(link_ci))

}