# R/pclm_CI.R In ungroup: Penalized Composite Link Model for Efficient Estimation of Smooth Distributions from Coarsely Binned Data

#### Documented in AIC.pclmAIC.pclm2DBIC.pclmBIC.pclm2Dcompute_standard_errorspclm.confidencepclm.confidence.dxpclm.confidence.mx

```# --------------------------------------------------- #
# Author: Marius D. Pascariu
# 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)
```

## Try the ungroup package in your browser

Any scripts or data that you put into this service are public.

ungroup documentation built on June 28, 2021, 5:08 p.m.