Computational Scaling for Large Joint Model Cohorts

Why Computational Cost Depends on the Likelihood Structure

A joint model can be statistically manageable yet computationally expensive when it is fitted to many subjects, many repeated measurements, or a complicated event-time contribution. In MCMC, a candidate parameter value must be assessed by its posterior density, and often by its gradient as well. The key scaling question is therefore not only how many parameters the model has, but how much data must be revisited for each proposed update.

The earlier material on the MCMC target in Bayesian joint models establishes the posterior target, while the earlier material on subject-specific random effects establishes subject-level updates. Here the emphasis is computational: count the evaluations required, identify which subject contributions actually change, and organize calculations so the same quantities are not repeatedly recomputed. Operation counts are estimates of relative workload, not predictions of exact elapsed time; hardware, software, memory access, and implementation also matter.

Definition: For $n$ subjects, let subject $i$ have longitudinal observations $y_{ij}$ at times $t_{ij}$, for $j=1,\ldots,m_i$, event or censoring time $T_i$, event indicator $\delta_i$, and random effects $b_i$. Write $M=\sum_{i=1}^n m_i$ for the total number of longitudinal measurements. Conditional on shared parameters $\theta$ and the random effects, a standard factorization of the log likelihood is $$ \ell(\theta,\{b_i\})=\sum_{i=1}^n \ell_i(\theta,b_i), \qquad \ell_i(\theta,b_i)=\ell_{y,i}(\theta,b_i)+\ell_{T,i}(\theta,b_i). $$ Here $\ell_{y,i}$ is the longitudinal log-likelihood contribution and $\ell_{T,i}$ is the event-time log-likelihood contribution. The posterior log density adds prior terms to this sum. This factorization assumes subjects are conditionally independent given the model parameters and random effects, with each subject's observations and event outcome represented by that subject's contribution.

For example, in a proportional-hazards formulation with hazard $h_i(t\mid b_i,\theta)$ and cumulative hazard $H_i(T_i\mid b_i,\theta)=\int_0^{T_i}h_i(u\mid b_i,\theta)\,du$, the event-time contribution is

\ell_{T,i}=\delta_i\log h_i(T_i\mid b_i,\theta)-H_i(T_i\mid b_i,\theta).

Thus an event contribution may require a hazard evaluation at the observed time and a cumulative-hazard calculation. The latter can be inexpensive in a simple model or costly when it uses numerical integration. Censoring changes the event indicator, not the need to account for survival up to the observed time.

First-Order Cost Estimates

Suppose evaluating one longitudinal observation costs roughly $c_y$ units of work, and evaluating a subject's event contribution costs roughly $c_T$. A full likelihood evaluation then has leading cost

C_{\mathrm{full}}\approx c_yM+c_Tn.

If the event term uses $R$ quadrature points per subject, its cost may itself grow with $R$: a simple approximation is $c_T\approx c_0+c_qR$, where $c_0$ covers fixed work and $c_q$ is the cost per integration point. If event-time calculations also evaluate a trajectory at several times, those evaluations belong in the cost count too. The important distinction is between a cost proportional to the number of observations and a cost proportional to the number of subjects, with additional factors for numerical integration or model complexity.

A single MCMC iteration can contain several parameter updates. If $K$ global proposals each require a full likelihood evaluation, those evaluations alone cost approximately $K C_{\mathrm{full}}$. Conversely, when updating only $b_i$, the likelihood ratio depends on the changed subject's contribution $\ell_i$, plus the prior contribution for $b_i$. Recomputing every other subject's likelihood is unnecessary. A sweep through all $n$ subject-specific effects still requires work across the cohort, but can be organized to scale like $O(M+n c_T)$, rather than $n$ full-cohort evaluations.

Worked Example: Counting the Cost of Global Evaluations

Consider hypothetical data from $n=2{,}400$ subjects, each with $m_i=8$ measurements. Suppose evaluating one longitudinal observation takes 6 model-evaluation units. For the event contribution, assume 12 integration points at 10 units each, plus 8 units for the event-time and fixed subject-level calculations. Use 2 additional units per subject for other bookkeeping.

The total number of longitudinal measurements is $M=2{,}400(8)=19{,}200$. The longitudinal work is therefore $19{,}200(6)=115{,}200$ units. Event work is $2{,}400[12(10)+8]=2{,}400(128)=307{,}200$ units. Bookkeeping adds $2{,}400(2)=4{,}800$ units, so one full evaluation costs about $115{,}200+307{,}200+4{,}800=427{,}200$ units.

If four global proposals each require a full evaluation, their combined cost is approximately $4(427{,}200)=1{,}708{,}800$ units. This calculation does not say the iteration takes a particular number of seconds; it identifies which component dominates under the stated assumptions. Here event calculations are the largest component, so reducing redundant hazard or integration work may matter more than optimizing the relatively smaller bookkeeping term.

Organizing Subject-Level Calculations

The factorization by subject suggests a practical organization. Keep data for each subject together, write a routine that returns that subject's longitudinal and event contributions, and combine the contributions when a global parameter changes. For a local proposal to $b_i$, evaluate the proposed and current values for subject $i$ only. For a shared parameter proposal, contributions generally change for every subject, so the full sum must be updated unless valid cached quantities or special algebra reduce the work.

A safe implementation can maintain the current vector of subject contributions. If a local update changes $b_i$, replace the stored $\ell_i$ with its proposed value and update the total by the difference. If the proposal is rejected, retain the old value. This bookkeeping preserves the same likelihood sum as a fresh calculation, while avoiding repeated recomputation of unchanged subjects. It is essential that every quantity affected by the proposal be included in the replaced contribution.

Precompute fixed data features.
Store observed times, design-matrix entries, interval or quadrature indices, and other quantities that do not depend on sampled parameters.

Evaluate a subject contribution.
Given current shared parameters and $b_i$, calculate the longitudinal and event terms for subject $i$, including all model-dependent pieces.

Update only what changed.
For a local random-effects proposal, recalculate that subject's contribution. For a shared-parameter proposal, recalculate all affected subjects or update them using an algebraically equivalent cached representation.

Check against a direct calculation.
For a small test data set, compare the optimized total log likelihood with a direct sum of subject contributions at the same parameter values.

Vectorization can reduce overhead by evaluating the same formula over a batch of subjects or measurements. It does not, by itself, change the leading order of work: evaluating all $M$ observations remains proportional to $M$. Parallel evaluation across subjects is also natural because the conditional likelihood contributions are separate. Parallelism can reduce wall-clock time, but adds communication, synchronization, and memory costs; it does not reduce the total number of likelihood calculations.

Worked Example: Local Updates Versus Repeated Full Scans

Suppose a hypothetical cohort has $n=1{,}000$ subjects. Evaluating the complete contribution for one subject, including its measurements and event term, takes 50 units. A proposed update to one subject's random effects changes only that subject's contribution, under the conditional-independence factorization.

A subject-wise update therefore costs about 50 units for its likelihood calculation. Updating all 1,000 random effects once each costs approximately $1{,}000(50)=50{,}000$ units, apart from prior calculations and other sampler overhead.

A wasteful implementation might recompute all 1,000 subject contributions for every subject-specific proposal. Each full scan costs $1{,}000(50)=50{,}000$ units, and 1,000 such scans cost $1{,}000(50{,}000)=50{,}000{,}000$ units. The ratio is $50{,}000{,}000/50{,}000=1{,}000$. The optimized and wasteful calculations target the same posterior if implemented correctly; the latter repeats unchanged work a thousandfold in this example.

When Caching Can Change the Workload

Some likelihood calculations admit summaries that can be computed once and reused. For instance, with a Gaussian longitudinal model, fixed design matrix $X_i$, and a mean linear in coefficients $\beta$, the residual sum of squares satisfies

(y_i-X_i\beta)^\mathsf{T}(y_i-X_i\beta) =y_i^\mathsf{T}y_i-2\beta^\mathsf{T}X_i^\mathsf{T}y_i+\beta^\mathsf{T}X_i^\mathsf{T}X_i\beta.

The quantities $y_i^\mathsf{T}y_i$, $X_i^\mathsf{T}y_i$, and $X_i^\mathsf{T}X_i$ can be precomputed if the data and design are fixed. Then a likelihood evaluation can use small matrix operations rather than revisiting every measurement. This can be useful when the coefficient dimension is small relative to $m_i$. It does not apply unchanged if the design depends on sampled quantities, the observation model is not Gaussian, or the mean is nonlinear in the parameters. And it does not automatically accelerate an event contribution that depends on the trajectory at event or integration times.

Worked Example: A Gaussian Summary for Repeated Evaluations

For one hypothetical subject, suppose a Gaussian model has 10 fixed-design measurements and coefficient vector $\beta=(1,0.5)^\mathsf{T}$. Assume the precomputed summaries are $y_i^\mathsf{T}y_i=82$, $X_i^\mathsf{T}y_i=(20,14)^\mathsf{T}$, and $X_i^\mathsf{T}X_i=\begin{pmatrix}10&5\\5&5\end{pmatrix}$. These summaries are treated as fixed inputs calculated from the subject's measurements and design.

First, $\beta^\mathsf{T}X_i^\mathsf{T}y_i=1(20)+0.5(14)=27$. Next, $\beta^\mathsf{T}X_i^\mathsf{T}X_i\beta=10(1^2)+2(5)(1)(0.5)+5(0.5^2)=10+5+1.25=16.25$. The residual sum of squares is therefore $82-2(27)+16.25=44.25$.

A direct calculation would form and square 10 residuals each time this subject's Gaussian longitudinal likelihood is needed. With these fixed-design summaries, the residual sum of squares instead requires a few operations on a two-dimensional coefficient vector. This is an algebraic shortcut for the specified Gaussian component, not a general replacement for evaluating all terms in a joint model.

Diagnostics, Alternatives, and Limitations

Timing alone can be misleading: an apparently fast implementation may omit a likelihood term, while a correct implementation may be slow because it repeatedly performs expensive integration. Test numerical agreement between a direct and optimized likelihood on small cases, including censored subjects and subjects with different measurement counts. Profile the implementation by component: count calls to longitudinal, hazard, cumulative-hazard, and gradient calculations, and measure their time separately. Check that cached values are invalidated whenever a proposal changes a quantity on which they depend.

The main organizational choices have different trade-offs. A direct full-cohort calculation is simple and useful as a reference implementation, but can repeat work during local updates. Subject-wise evaluation avoids that repetition and makes the conditional factorization explicit. Vectorization often improves throughput but may use more temporary memory. Parallel subject evaluation can help when each contribution is costly and the cohort is large, but overhead can erase the gain for small problems. Caching can substantially reduce repeated work when valid fixed summaries exist, but increases implementation complexity and creates risks of stale values.

Memory also scales with the amount of retained data and the quantities stored for computation. The observations themselves require storage proportional to $M$; subject-level parameters and contributions require storage proportional to $n$. Saving every random effect at every retained MCMC draw can require storage proportional to the product of the number of draws, $n$, and the random-effects dimension. Storing only summaries or selected draws reduces memory use but changes what can be inspected after sampling. The computational cost of running more iterations is distinct from the statistical adequacy of those iterations; a lower cost per iteration does not establish that a chain explores the posterior well.

A common mistake is to assume that a subject-wise factorization means every MCMC update costs only $O(1)$, regardless of cohort size. That is true only for a single local contribution when its data size is bounded. A full sweep across all subjects still scales with the total amount of data, and a shared-parameter update typically affects every subject. Another mistake is to treat precomputation as permission to omit terms that depend on sampled parameters. The correction is to identify exactly what changes under each proposal, preserve the full posterior target, and verify optimized calculations against a direct implementation.

Key takeaway: Estimate cost by separating longitudinal work, event-time work, and sampler overhead; then count how often each is recomputed. Subject-level factorization makes local updates efficient, while global updates still require all affected contributions. Caching, vectorization, and parallelism help only when they preserve the exact target calculation.

Check Your Understanding

Answer each question, then compare with the solutions.

A cohort has 600 subjects with 5 measurements each. How many longitudinal observations are there?
If one subject's random-effects proposal changes no other subject's likelihood contribution, which likelihood terms need to be recalculated for its acceptance ratio?
A full likelihood evaluation costs 80,000 units. What is the approximate cost of three independent global proposals, each requiring a full evaluation?
Name one condition needed for the Gaussian residual-summary calculation to remain valid.
Why does parallel evaluation across subjects not reduce the total number of likelihood calculations?

Solutions. (1) $M=600(5)=3{,}000$ observations. (2) Recalculate the proposed and current likelihood contribution for that subject, together with the relevant random-effects prior terms; unchanged subjects cancel. (3) $3(80{,}000)=240{,}000$ units. (4) For example, the Gaussian observation model and fixed design matrix must support the stated quadratic residual identity, and the cached summaries must not depend on the proposed coefficients. (5) Parallelism distributes the same subject contributions across processors; it may reduce elapsed time, but the total work remains the sum of those calculations.