pclm: Fit a composite link model

In JOPS: Practical Smoothing with P-Splines

Fit a composite link model

Description

Fit a smooth latent distribution using the penalized composite link model (PCLM).

Usage

pclm(y, C, B, lambda = 1, pord = 2, itmax = 50, show = FALSE)

Arguments

y	a vector of counts, length m.
C	a composition matrix, m by q.
B	a B-spline basis matrix, q by n.
lambda	the penalty parameter.
pord	the the order of the difference penalty (default = 2).
itmax	the maximum number of iterations (default = 50).
show	Set to TRUE or FALSE to display iteration history (default = FALSE).

Details

The composite link model assumes that E(y) = \mu = C\exp(B \alpha), where \exp(B\alpha) is a latent discrete distribution, usually on a finer grid than that for y.

Note that sum(gamma) == sum(mu).

Value

A list with the following items:

alpha	the estimated B-spline coefficients, length n.
gamma	the estimated latent distribution, length q.
mu	estimated values of y, length m.
dev	the deviance of the model.
ed	the effective model dimension.
aic	Akaike's Information Criterion.

Author(s)

Paul Eilers and Jutta Gampe

References

Eilers, P. H. C. (2007). III-posed problems with counts, the composite link model and penalized likelihood. Statistical Modelling, 7(3), 239–254.

Eilers, P.H.C. and Marx, B.D. (2021). Practical Smoothing, The Joys of P-splines. Cambridge University Press.

Examples

# Left and right boundaries, and counts, of wide intervals of the data
cb <- c( 0, 20, 30, 40, 50, 60)
ce <- c(20, 30, 40, 50, 60, 70)
y <- c(79, 54, 19, 1, 1, 0)

# Construct the composition matrix
m <- length(y)
n <- max(ce)
C <- matrix(0, m, n)
for (i in 1:m) C[i, cb[i]:ce[i]] <- 1

mids = (cb + ce) / 2 - 0.5
widths = ce - cb + 1
dens = y / widths / sum(y)
x = (1:n) - 0.5
B = bbase(x)
fit = pclm(y, C, B, lambda = 2, pord = 2, show = TRUE)
gamma = fit$gamma / sum(fit$gamma)
# Plot density estimate and data
plot(x, gamma, type = 'l', lwd = 2, xlab = "Lead Concentration", ylab = "Density")
rect(cb, 0, ce, dens, density = rep(10, 6), angle = rep(45, 6))

JOPS documentation built on Sept. 8, 2023, 5:42 p.m.