Spline Trajectories for Repeated Counts

Why Combine Splines and a Count Model?

Repeated counts can change nonlinearly over follow-up: they may decline, level off, and later rise. A straight-line time effect can miss such patterns, while a separate model for the event outcome does not account for the connection between a subject’s evolving count trajectory and event risk. A spline-based joint model addresses both features: it allows a flexible time course for the counts and links the event hazard to a latent feature of that trajectory.

The focus here is a count measured repeatedly over equal observation windows. For unequal exposure periods, the count model should include the log exposure duration as an offset; otherwise a longer window could appear to have a higher count simply because it is longer. The event component describes time to event, with right-censoring represented by the observed follow-up time and event indicator.

Model Specification

Let $Y_{ij}$ be the count for subject $i$ at visit time $t_{ij}$, and let $T_i$ be the event or censoring time. Write $\delta_i=1$ if the event is observed and $\delta_i=0$ if follow-up is censored. Let $B(t)$ be a chosen spline basis, and let $b_i$ be subject-specific random effects. A convenient latent log-mean trajectory is

m_i(t)=B(t)^\mathsf{T}\beta+z(t)^\mathsf{T}b_i, \qquad \mu_i(t)=\exp\{m_i(t)\}.

Here $\mu_i(t)$ is the conditional expected count at time $t$, given the subject’s random effects. The fixed spline coefficients $\beta$ describe the reference trajectory, while $z(t)^\mathsf{T}b_i$ allows subjects to differ. For instance, a random intercept makes one subject’s trajectory consistently higher or lower on the log-mean scale.

Definition: In this tutorial’s joint model, conditional on $b_i$, the repeated count $Y_{ij}$ has a negative-binomial distribution with mean $\mu_i(t_{ij})$ and size parameter $\kappa>0$: $$ E(Y_{ij}\mid b_i)=\mu_i(t_{ij}), \qquad \operatorname{Var}(Y_{ij}\mid b_i)=\mu_i(t_{ij})+\frac{\mu_i(t_{ij})^2}{\kappa}. $$ The event hazard is $$ h_i(t\mid b_i)=h_0(t)\exp\{\gamma^\mathsf{T}w_i+\alpha m_i(t)\}. $$ The parameter $\alpha$ links the latent log mean to the event hazard; $w_i$ denotes baseline covariates and $h_0(t)$ is the baseline hazard.

This negative-binomial parameterization allows conditional variance to exceed the conditional mean. As $\kappa$ grows without bound, the variance approaches the mean and the distribution approaches a Poisson distribution. The random effects are commonly assigned a mean-zero multivariate normal distribution with covariance matrix $D$. Inference then accounts for the shared random effects rather than treating repeated counts from the same subject as independent unconditionally.

The model assumes the count observations are conditionally independent given the subject’s random effects and model covariates, and that the event process and count process are conditionally independent given those same random effects and covariates. Their dependence in the joint model arises through their shared $b_i$ and the trajectory $m_i(t)$. The censoring mechanism must be noninformative given the modeled information, and the visit process must not carry unmodeled information about the count or event process. These assumptions are substantive: for example, event-driven extra visits may require modeling the visit process or otherwise accounting for its informativeness.

Interpreting the Spline Trajectory

A spline is a function built from basis terms that can change its slope or curvature across time while remaining suitably joined at knots. For a simple piecewise-linear spline with a knot at month 6, one possible basis is $B(t)=(1,t,(t-6)_+)^\mathsf{T}$, where $(t-6)_+=\max(t-6,0)$. With no random effects, the log mean is linear before and after the knot, with slopes that may differ. More complex spline bases can represent smooth curvature.

The spline is defined on the log-mean scale in this model. Exponentiating it gives the expected count. Thus, a coefficient or slope on the log-mean scale is not an additive change in the expected count. A log-mean change of $d$ corresponds to a mean-count ratio of $\exp(d)$. Also, $B(t)^\mathsf{T}\beta$ describes the trajectory for $b_i=0$, not generally the population-average expected count: averaging $\exp\{m_i(t)\}$ over random effects need not equal $\exp\{B(t)^\mathsf{T}\beta\}$.

Worked Example: Reading a Piecewise-Linear Count Trajectory

Suppose hypothetical fitted fixed effects, with time in months and a knot at month 6, give $m(t)=1.10-0.08t+0.12(t-6)_+$ for the reference subject $b_i=0$. Counts are measured over equal observation windows. Calculate the expected counts at months 0, 6, and 10, and interpret the segment slopes.

At month 0, $(0-6)_+=0$, so $m(0)=1.10$ and $\mu(0)=\exp(1.10)=3.004$ counts per window. At month 6, the knot term is also zero, giving $m(6)=1.10-0.08(6)=0.62$ and $\mu(6)=\exp(0.62)=1.859$ counts per window.

At month 10, $(10-6)_+=4$, so $m(10)=1.10-0.08(10)+0.12(4)=0.78$. Therefore, $\mu(10)=\exp(0.78)=2.181$ counts per window. The fitted reference trajectory falls from about 3.00 counts at month 0 to 1.86 at month 6, then rises to about 2.18 at month 10.

Before month 6, the log-mean slope is $-0.08$ per month, corresponding to a mean-count ratio of $\exp(-0.08)=0.9231$ for a one-month increase, or about a 7.69% decrease. After month 6, the log-mean slope is $-0.08+0.12=0.04$ per month; the corresponding ratio is $\exp(0.04)=1.041$, or about a 4.08% increase per month. These are relative changes in the conditional mean, not changes of 0.08 or 0.04 counts.

Dispersion and the Shared Event Association

The size parameter $\kappa$ governs extra-Poisson variation conditional on the random effects. A smaller $\kappa$ produces greater variance for a given mean. This addresses count variability beyond the Poisson variance, but does not automatically correct a wrong mean trajectory, unmodeled zero-generating process, or dependence remaining after conditioning on the random effects.

In the event component, $\alpha$ is the association between the latent log mean and the instantaneous hazard. At the same time $t$, for subjects with the same baseline covariates and latent log means differing by $d$, the conditional hazard ratio is $\exp(\alpha d)$. Since $d$ is a difference in log expected counts, a one-unit increase in $m_i(t)$ means the expected count is multiplied by $e$, not increased by one count. The association is conditional and is not, by itself, evidence that changing the count causes a change in event risk.

Worked Example: Translating the Association Into a Hazard Ratio

Suppose a hypothetical fitted model has $\alpha=0.70$. At month 6, two subjects with the same baseline covariates have latent log means that differ by $0.50$, with subject B higher. Calculate the expected-count ratio and the conditional hazard ratio.

The expected-count ratio for B versus A is $\exp(0.50)=1.649$. Thus B’s conditional expected count is about 1.65 times A’s at that time. The corresponding conditional hazard ratio is $\exp\{0.70(0.50)\}=\exp(0.35)=1.419$.

Under the model, B’s instantaneous event hazard at month 6 is therefore about 1.42 times A’s, conditional on the stated covariates and latent trajectories. This is a comparison of instantaneous hazards, not a 42% increase in event probability over a specified follow-up period.

What the Joint Likelihood Connects

The shared trajectory matters because the count data inform the subject-specific random effects, and those random effects also enter the event hazard. Conditional on $b_i$, the event contribution to the likelihood is the hazard at the observed time raised to the event indicator, multiplied by the conditional probability of remaining event-free through that time. The count and event contributions are then combined and integrated over the random-effects distribution. In schematic form,

L_i=\int \left\{\prod_j f_{\mathrm{NB}}\bigl(Y_{ij}\mid \mu_i(t_{ij}),\kappa\bigr)\right\} h_i(T_i\mid b_i)^{\delta_i} \exp\left\{-\int_0^{T_i}h_i(u\mid b_i)\,du\right\} f_b(b_i)\,db_i .

Here $f_{\mathrm{NB}}$ is the negative-binomial probability mass function under the stated mean-size parameterization, and $f_b$ is the random-effects density. The integral expresses the fact that $b_i$ is unobserved. It also shows why fitting a count trajectory first and then treating estimated subject-specific values as known event predictors does not generally reproduce joint-model inference: it omits uncertainty in those latent effects and can distort the association estimate.

Worked Example: Checking a Negative-Binomial Count Probability

For illustration, take a hypothetical fitted conditional mean of $\mu=2.181$ counts per window and size $\kappa=2.5$. Calculate the conditional variance and the probability of observing exactly four counts. Under this parameterization, $P(Y=y)=\frac{\Gamma(y+\kappa)}{\Gamma(\kappa)y!}p^\kappa(1-p)^y$, where $p=\kappa/(\kappa+\mu)$.

The variance is $\mu+\mu^2/\kappa=2.181+2.181^2/2.5=2.181+1.903=4.084$ counts squared (rounding differences from the unrounded mean are negligible). This exceeds the mean, as expected for finite $\kappa$.

For the probability calculation, $p=2.5/(2.5+2.181)=0.5341$. Substituting $y=4$, \[ P(Y=4)= \frac{\Gamma(6.5)}{\Gamma(2.5)4!}(0.5341)^{2.5}(0.4659)^4 \approx 9.023(0.2085)(0.0471) \approx 0.0886. \] Thus, given the specified conditional mean and size, the probability of exactly four counts in a window is about 0.0886. It is a model-based probability for that count outcome, not an estimate of the event risk.

Diagnostics, Alternatives, and Limitations

A fitted spline can be flexible without being adequate. Inspect observed and fitted count patterns over time, including whether residual discrepancies persist in particular time regions or across subjects. Check whether the fitted dispersion is plausible and whether variability remains much larger than the modeled conditional variance. Review the event component as well: a poorly supported association estimate, unstable random-effects estimates, or a trajectory that follows a few observations may indicate limited information. Sparse visits near a knot and few observed events can make flexible trajectory and association estimates especially uncertain.

Compare the spline specification with a simpler time trend and with reasonable alternative knot placement or basis complexity. A nearly linear fitted trajectory may not justify extra spline terms. Conversely, systematic curvature in the residual pattern can signal that a linear time effect is inadequate. Spline behavior near the ends of follow-up deserves particular scrutiny because there may be little data to constrain it. These checks help evaluate the chosen model; they do not establish that the random-effects distribution, censoring assumption, or visit-process assumption is correct.

Approach	Strength	Important limitation
Linear-time negative-binomial joint model	Parsimonious and easier to interpret when change is approximately linear on the log-mean scale.	Can miss bends or changes in slope.
Spline negative-binomial joint model	Represents nonlinear time trends while retaining a count likelihood and shared event association.	Requires basis and knot choices; flexible fits need adequate repeated measurements and events.
Poisson spline joint model	Provides a simpler count distribution when conditional variance is close to the mean.	Can understate uncertainty when counts are overdispersed.

A negative-binomial distribution is not a remedy for every count-data feature. If zeros arise from a distinct structural process, a zero-inflated or hurdle formulation may be more appropriate, but it changes the longitudinal model and the meaning of the shared association. Likewise, the association here is specifically with the current latent log mean. Linking the hazard to a different trajectory feature, such as a slope, defines a different scientific model and a different interpretation. No spline specification supports reliable extrapolation far beyond the observed follow-up without additional assumptions.

Common Pitfalls

Reading a log-scale coefficient as an additive count change. A change of $d$ in $m_i(t)$ multiplies the expected count by $\exp(d)$; it does not add $d$ counts.
Interpreting $\alpha$ as a hazard ratio by itself. The hazard ratio for a specified trajectory contrast $d$ is $\exp(\alpha d)$, not $\alpha$.
Calling the reference trajectory the population-average count curve. The reference curve sets random effects to zero. In general, averaging the exponentiated trajectory over random effects gives a different curve.
Using a Poisson model despite clear overdispersion. Its conditional variance equals its conditional mean; a mismatch can affect uncertainty and inference.
Treating the fitted association as causal. It is a conditional model association and may reflect unmeasured factors or model misspecification.

Key takeaway: A spline supplies a flexible time pattern for the conditional log mean, the negative-binomial distribution accommodates extra-Poisson variation, and the event model can link hazard to that latent pattern. Interpret each component on its own scale and state the conditioning assumptions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For $m(t)=1.10-0.08t+0.12(t-6)_+$, calculate the reference expected count at month 8.
What is the log-mean slope after month 6, and what one-month expected-count ratio does it imply?
If the conditional mean is 3 and $\kappa=3$, what is the conditional negative-binomial variance? What happens to the variance as $\kappa$ tends to infinity?
If $\alpha=0.70$ and one subject’s latent log mean is 0.30 lower than another’s at the same time, calculate the first subject’s conditional hazard ratio relative to the second.

Solutions. (1) At month 8, $(8-6)_+=2$, so $m(8)=1.10-0.08(8)+0.12(2)=0.70$, and $\mu(8)=\exp(0.70)=2.014$ counts per window. (2) The post-knot slope is $-0.08+0.12=0.04$ per month, so the expected-count ratio per month is $\exp(0.04)=1.041$. (3) The variance is $3+3^2/3=6$ counts squared. As $\kappa$ tends to infinity, the variance tends to 3, the Poisson variance. (4) The log-mean difference is $-0.30$, so the hazard ratio is $\exp\{0.70(-0.30)\}=\exp(-0.21)=0.8106$. The first subject’s conditional instantaneous hazard is about 0.811 times the second’s under the model.