# pclm: Fit a composite link model In JOPS: Practical Smoothing with P-Splines

 pclm R Documentation

## 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.