Regression analysis of panel count data

In pcdreg: Semiparametric Regression for Panel Count Data

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6,
                      fig.height = 4)

library(pcdreg)

The data and the model

Panel count data arise when a recurrent event process is observed at periodic examinations. Each visit gives the number of events since the previous visit, but not their event times. Examples include hospitalisations counted at quarterly clinic visits, tumours counted at scheduled sacrifices, and insurance claims tallied at renewal.

Write $N_i(t)$ for the number of events subject $i$ has experienced by time $t$, observed at times $T_{i1} < \dots < T_{iJ_i}$, and $\Delta N_{ij}$ for the count in $(T_{i,j-1}, T_{ij}]$. pcdreg() fits the proportional rate model

$$E[dN(t) \mid X(t)] = \exp(\beta' X(t)) \, d\Lambda(t),$$

where $\Lambda$ is an unspecified non-decreasing baseline cumulative rate and $X(t)$ may vary over time.

The proportional means model is $E[N(t) \mid X(t)] = \mu(t)\exp(\beta'X(t))$. It requires $\mu$ to be non-decreasing. With covariates that fluctuate over time, this condition is difficult to arrange and interpret: a fitted means model can predict a declining mean count. The rate model requires only $\exp(\beta'X(t))\,d\Lambda(t)$ to be positive. Predicted means therefore increase automatically, and $\exp(\beta'X(t))$ has the instantaneous relative-risk interpretation of a hazard ratio in survival analysis.

The means model is the established alternative. The package also fits it with pcdreg(model = "mean"). The last section compares the models.

Getting data into shape

pcd() constructs the response. It plays the same role as Surv() in survival analysis and has two forms.

When covariates do not change over time, one row per examination is enough:

visits <- data.frame(
  id    = c(1, 1, 1, 2, 2),
  time  = c(0.5, 1.1, 1.9, 0.7, 1.4),
  count = c(2, 0, 3, 1, 4),
  dose  = c(10, 10, 10, 25, 25)
)
pcd(visits$id, visits$time, visits$count)

Interval starts are filled in from the previous examination of the same subject, with follow-up beginning at zero.

When covariates do change over time, use the counting process form: one row per interval with constant covariates. This is the layout coxph() uses for time-dependent covariates. A numeric count records an examination at tstop. An NA count records a covariate change only.

cp <- data.frame(
  id     = c(1, 1, 1, 2),
  tstart = c(0.0, 0.5, 0.9, 0.0),
  tstop  = c(0.5, 0.9, 1.2, 0.8),
  count  = c(2, NA, 0, 3),
  dose   = c(10, 10, 25, 25)
)
pcd(cp$id, cp$tstart, cp$tstop, cp$count)

Subject 1 was examined at 0.5 and 1.2, with two events in the first interval and none in the second. The dose changed at 0.9, within the second interval. The rate model uses this information; the means model cannot.

Within a subject the intervals must be contiguous, start at zero, and end at an examination. pcdreg() checks all of this and says which subject is at fault when it fails.

Converting data held in the older layout

An older layout repeats each panel count on every row of its interval and marks examinations with a separate indicator. Convert it in one line. A repeated count on a non-examination row carries no information beyond the examination row:

new <- old
new$count[old$event == 0] <- NA
pcdreg(pcd(id, tstart, tstop, count) ~ x1 + x2, data = new)

Fitting

set.seed(1)
d <- sim_pcd(200, beta = c(1, -1), lambda = function(t) 8 / (1 + t))
fit <- pcdreg(pcd(id, tstart, tstop, count) ~ x1 + x2, data = d)
summary(fit)

There is no intercept: a constant in $X$ would be indistinguishable from a rescaling of $\Lambda$, so pcdreg() drops it and reports the baseline separately.

Fitting maximises the likelihood implied by a nonhomogeneous Poisson process. The EM algorithm treats the unobserved per-time counts as missing data. This makes the integral $\int \exp(\beta'X(t)) d\Lambda(t)$ tractable and requires covariates only at observed examination times, not as complete trajectories.

The Poisson likelihood is a working device, not an assumption about the data. The estimator stays consistent and asymptotically normal when the counts are not Poisson at all.

Which standard errors

This distinction matters for the covariance. Three estimators are available.

se <- function(type) sqrt(diag(vcov(fit, type)))
rbind(robust = se("robust"), information = se("information"))

"robust", the default, is the sandwich $\Omega^{-1} S \Omega^{-1}/n$ built from the two matrices

fit$Omega
fit$S

Both matrices are by-products of the last EM iteration, so the robust standard errors are essentially free. The estimator is consistent whether or not the counts are Poisson.

"information" is $S^{-1}/n$. Under the Poisson assumption, $S$ is the efficient information and this is the right answer. It is also a cheap stand-in for the profile likelihood method, with which it agrees closely.

"profile" is the numerical profile likelihood estimator of Murphy and van der Vaart. It is included for comparison with the literature and is computed only on request because it is expensive:

pfit <- pcdreg(pcd(id, tstart, tstop, count) ~ x1 + x2, data = d,
                  profile = TRUE)
rbind(information = sqrt(diag(vcov(pfit, "information"))),
      profile     = sqrt(diag(vcov(pfit, "profile"))))

The two estimates track each other closely, as theory predicts, because both estimate $S$. They differ because the profile version uses a numerical derivative with a step of order $n^{-1/2}$ and therefore has a discretisation error that shrinks only with sample size. Prefer "information" when the Poisson assumption is tenable. "profile" is mainly provided for comparison with the literature.

What happens when the Poisson assumption fails

Real recurrent event data are usually overdispersed: some subjects are more event-prone than their covariates suggest. A gamma frailty preserves the rate model but destroys the Poisson structure.

set.seed(2)
od <- sim_pcd(200, beta = c(1, -1), frailty = 1)
odfit <- pcdreg(pcd(id, tstart, tstop, count) ~ x1 + x2, data = od)
rbind(robust = sqrt(diag(vcov(odfit, "robust"))),
      information = sqrt(diag(vcov(odfit, "information"))))

The information-based standard errors are several times too small. The paper's simulations put coverage of the resulting nominal 95% intervals around 15--22%. For this reason, "robust" is the default and other choices require care.

The baseline and predicted means

baseline() returns the estimated jumps and the cumulative baseline rate.

head(baseline(fit))

truth <- data.frame(time = fit$baseline$time,
                    value = 8 * log(1 + fit$baseline$time))
autoplot(fit) +
  ggplot2::geom_line(data = truth, colour = "#c0392b", linetype = 2)

The dashed line is the truth, $\Lambda(t) = 8\log(1+t)$.

Predicted mean counts follow each subject's covariate trajectory, $E[N(t) \mid X] = \int_0^t \exp(\beta'X(s)) \, d\hat\Lambda(s)$:

pred <- predict(fit, d, type = "mean")
head(pred)

one <- pred[pred$id == 1, ]
ggplot2::ggplot(one, ggplot2::aes(x = time, y = mean)) +
  ggplot2::geom_step(colour = "#2a78d6", linewidth = 0.6) +
  ggplot2::labs(x = "Time", y = "Predicted mean count", title = "Subject 1") +
  ggplot2::theme_minimal(base_size = 11)

The curve is non-decreasing even though x1 changes part way through follow-up. This is a structural advantage of the rate model.

Tidy output

tidy(), glance() and augment() methods are provided for the broom generics.

tidy(fit, conf.int = TRUE)
glance(fit)

exponentiate = TRUE reports $e^{\beta}$, the multiplicative effect on the rate, read as a hazard ratio is read.

tidy(fit, conf.int = TRUE, exponentiate = TRUE)

augment() attaches the fitted mean to each row, together with the observed cumulative count and their difference. .resid is the excess of events seen so far over the number the model implies, so it is the quantity to plot against time or against a covariate when checking fit.

aug <- augment(fit)
head(aug[!is.na(aug$.observed),
         c("id", "tstop", ".fitted", ".observed", ".resid")])

The means model

pcdreg(model = "mean") fits the other classical model for these data,

$$E[N(t) \mid X(t)] = \mu(t)\exp(\beta'X(t)),$$

by the estimating equation of Hu, Sun and Wei (2003). It is the comparator in the paper's application. It is inexpensive because only examination times enter; the grid expansion required by the rate model is unnecessary.

mfit <- pcdreg(pcd(id, tstart, tstop, count) ~ x1 + x2, data = d, model = "mean")
summary(mfit)

There is no likelihood behind this estimator, so there is no log likelihood to report and vcov() offers only the sandwich.

The two fits require care in interpretation.

cbind(rate = coef(fit), mean = coef(mfit))

These numbers do not estimate one quantity. With time-varying covariates, the models are not reparametrisations of each other: $\beta$ acts on the instantaneous rate in one and the cumulative mean in the other. They answer different questions and need not agree. These data were generated from the rate model, so pcdreg() recovers the simulation values and pcdreg(model = "mean") does not.

Second, nothing constrains the fitted $\hat\mu$ to increase:

mu <- baseline(mfit)$mean
c(increasing = !is.unsorted(mu), decreases_at = sum(diff(mu) < 0))

This is the difficulty that the package's rate model addresses. It appears in ordinary simulated data. A declining mean function is not a numerical failure. It is the estimator's result when covariate values move up and down. Predicted means from pcdreg(model = "mean") can therefore decrease, while those from pcdreg() cannot.

Computation

The EM algorithm converges linearly and slowly on this problem: plain EM can still be changing in the fifth decimal place after tens of thousands of passes. pcdreg() extrapolates the iterations with the SQUAREM scheme of Varadhan and Roland. An extrapolated point is kept only if a stabilising EM pass from it does at least as well on the observed-data log likelihood. The fixed point is unchanged. The required number of passes typically falls by an order of magnitude.

c(iterations = fit$iterations, EM_passes = fit$passes,
  converged = fit$converged)

Set accelerate = FALSE in pcdreg_control() to recover plain EM.

The dominant cost is the pooled grid of distinct examination times. The covariate trajectory for every subject is evaluated at every grid time within follow-up. Tied examination times across subjects, as with scheduled visits, keep the grid small.

References

Sun, D., Guo, Y., Li, Y., Tu, W. and Sun, J. (2024). A robust approach for regression analysis of panel count data. Bernoulli 30(4), 3251--3275. doi:10.3150/23-BEJ1713

Hu, X. J., Sun, J. and Wei, L. J. (2003). Regression parameter estimation from panel counts. Scandinavian Journal of Statistics 30(1), 25--43. doi:10.1111/1467-9469.00316

Varadhan, R. and Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics 35, 335--353.

pcdreg documentation built on Nov. 10, 2026, 5:09 p.m.