R/pclm_CI.R

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

# --------------------------------------------------- #
# Author: Marius D. Pascariu
# Last update: Fri Oct 02 17:06:29 2026
# --------------------------------------------------- #


#' Compute Standard Errors and Confidence Intervals
#' @param B B-spline basis matrix
#' @param QmQ Object generated by pclm.fit
#' @param H0 Inverse of the penalized system, taken from pclm.fit so that the
#'   system is factorized once per fit rather than twice.
#' @keywords internal
compute_standard_errors <- function(B, QmQ, H0) {
  H1 <- H0 %*% QmQ %*% H0 # vcov matrix sandwich estimator
  # Only the diagonal of B H1 B' is wanted. Building the whole n x n product to
  # read it costs O(n^2) memory where O(n kB^2) of work would do.
  SE <- sqrt(rowSums((B %*% H1) * 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, qn = qn)
    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, qn = qn)
      
    } 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, qn = qn)
    }
    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 ("deep" %in% names(object)) 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
#' Sakamoto Y, Ishiguro M, Kitagawa G (1986).
#' \emph{Akaike Information Criterion Statistics}.
#' KTK Scientific Publishers.
#' ISBN 9027722536.
#' @keywords internal
#' @export
AIC.pclm2D <- function(object, ..., k = 2) AIC.pclm(object, k = k)


#' PCLM Bayes Information Criterion
#' 
#' @inherit AIC.pclm description params details return
#' @references
#' Sakamoto Y, Ishiguro M, Kitagawa G (1986).
#' \emph{Akaike Information Criterion Statistics}.
#' KTK Scientific Publishers.
#' ISBN 9027722536.
#' @keywords internal
#' @export
BIC.pclm <- function(object, ...) {
  # Neither stats::BIC nor this method has a k. Anything arriving here, such as
  # the k that BIC.pclm2D used to declare and then drop, would be swallowed by
  # ... without a word.
  if (...length() != 0L) {
    stop("BIC() for a PCLM fit takes no extra arguments", call. = FALSE)
  }
  X <- if ("deep" %in% names(object)) object$deep else object
  with(X, {
    bic <- dev + log(ny_) * trace
    return(bic)
  })
}


#' PCLM-2D Bayes Information Criterion
#'
#' @inherit BIC.pclm description params details return
#' @references
#' Sakamoto Y, Ishiguro M, Kitagawa G (1986).
#' \emph{Akaike Information Criterion Statistics}.
#' KTK Scientific Publishers.
#' ISBN 9027722536.
#' @keywords internal
#' @export
BIC.pclm2D <- function(object, ...) BIC.pclm(object, ...)