R/pclm_CI.R

In ungroup: Penalized Composite Link Model for Efficient Estimation of Smooth Distributions from Coarsely Binned Data

# --------------------------------------------------- #
# Author: Marius D. Pascariu
# License: MIT
# Last update: Thu Nov 07 11:46:20 2019
# --------------------------------------------------- #


#' Compute Standard Errors and Confidence Intervals
#' @param B B-spline basis matrix
#' @param QmQ Object generated by pclm.fit
#' @param QmQP Object generated by pclm.fit
#' @keywords internal
compute_standard_errors <- function(B, QmQ, QmQP) {
  H0 <- solve(QmQP)       # vcov matrix Bayesian approach
  H1 <- H0 %*% QmQ %*% H0 # vcov matrix sandwich estimator
  SE <- sqrt(diag(B %*% H1 %*% t(B)))
  return(SE)
}


#' Compute Confidence Intervals for PCLM output
#' 
#' @param fit Fitted values from pclm.fit
#' @param SE Standard Errors
#' @inheritParams pclm
#' @inheritParams pclm.fit
#' @keywords internal
pclm.confidence <- function(fit, out.step, y, SE, ci.level, type, offset) {
  I    <- c(as.list(environment()))
  I$ny <- length(y) 
  cil  <- 1 - ci.level/100
  I$qn <- qnorm(1 - cil/2)
  
  if (is.null(offset)) {
    out <- pclm.confidence.dx(I)
  } else {
    out <- pclm.confidence.mx(I)
  }
  return(out)
}


#' Compute Confidence Intervals for estimated hazard 
#' @param X List of inputs from pclm.confidence function
#' @keywords internal
pclm.confidence.mx <- function(X) {
  with(X, {
    lower <- as.numeric(exp(log(fit) - qn * SE))
    upper <- as.numeric(exp(log(fit) + qn * SE))
    
    if (type == "2D") {
      fit   <- matrix(fit, ncol = ny)
      upper <- matrix(upper, ncol = ny)
      lower <- matrix(lower, ncol = ny)
      SE    <- matrix(SE, ncol = ny)
    }
    out <- list(fit = fit, lower = lower, upper = upper, SE = SE)
    return(out)
  })
}


#' Compute Confidence Intervals for estimated distribution
#' @inheritParams pclm.confidence.mx
#' @keywords internal
pclm.confidence.dx <- function(X) {
  with(X, {
    if (X$type == "2D") {
      SE  <- matrix(SE, ncol = ny)
      fit <- matrix(fit, ncol = ny)
      Xi  <- X
      L   <- U <- fit * 0
      for (i in 1:ny) {
        Xi$fit <- fit[, i]
        Xi$SE  <- SE[, i]
        Xi$type <- "1D"
        ci <- pclm.confidence.dx(Xi)
        L[, i] <- ci$lower
        U[, i] <- ci$upper
      }
      out <- list(fit = fit, lower = L, upper = U, SE = SE)
      
    } else {
      Qn  <- t(as.matrix(c(-1, 1) * qn))
      fit <- as.numeric(fit)
      lx  <- rev(cumsum(rev(fit))) # dx to lx
      Lx  <- lx*out.step - fit*out.step/2
      mxu <- exp(log(fit) + SE %*% Qn)/Lx # compute CI for hazard rates
      qxu <- 1 - exp(-out.step*mxu)       # hazard to qx
      
      qx_lx_dx <- function(qx) {
        N <- length(qx)
        lx[1] * c(1, cumprod(1 - qx)[1:(N - 1)]) * qx 
      }
      
      dxU <- apply(qxu, 2, qx_lx_dx)       # qx to lx to dx
      # # make sure that the size of the population is the same
      dxU <- apply(dxU, 2, FUN = function(dx.hat) dx.hat * sum(fit)/sum(dx.hat))
      out <- list(fit = fit, lower = dxU[, 1], upper = dxU[, 2], SE = SE)
    }
    return(out)
  })
}


#' PCLM Akaike Information Criterion
#' 
#' @inherit stats::AIC description params details return references
#' @keywords internal
#' @export
AIC.pclm <- function(object, ..., k = 2) {
  X <- if (any(ls(object) == "deep"))  object$deep else object
  with(X, {
    aic <- dev + k * trace
    return(aic)  
  })
}


#' PCLM-2D Akaike Information Criterion
#' 
#' @inherit AIC.pclm description params details return references
#' @keywords internal
#' @export
AIC.pclm2D <- function(object, ..., k = 2) AIC.pclm(object, ..., k)


#' PCLM Bayes Information Criterion
#' 
#' @inherit AIC.pclm description params details return references
#' @keywords internal
#' @export
BIC.pclm <- function(object, ...) {
  X <- if (any(ls(object) == "deep"))  object$deep else object
  with(X, {
    bic <- dev + log(ny_) * trace
    return(bic)
  })
}


#' PCLM-2D Akaike Information Criterion
#' 
#' @inherit AIC.pclm description params details return references
#' @keywords internal
#' @export
BIC.pclm2D <- function(object, ..., k = 2) BIC.pclm(object, ..., k)