Summary of the Stan Model • epikinetics

Introduction

This document provides a summary of a hierarchical statistical model of antibody kinetics, implemented in Stan. The model is designed to analyse longitudinal titre data, accounting for boosting and waning effects over time. It incorporates individual-level random effects and covariate influences on key model parameters. A full description of the model can be found in the supplementary material in the published paper using the methods described here: Russell TW et al. Real-time estimation of immunological responses against emerging SARS-CoV-2 variants in the UK: a mathematical modelling study. Lancet Infect Dis. 2024 Sep 11:S1473-3099(24)00484-5.

Usage Notes

This model is suitable for analysing longitudinal titre data with either lower, upper, both (or no) censoring and incorporates both individual-level variability and arbitrary regression structure to adjust for covariates. Users can adapt the model by specifying appropriate covariates via a R style linear model formula, priors, and data inputs relevant to their study.

Model Specification

Overview

The model describes the expected log-transformed titre value $\mu_{n}$ for individual $n$ at time $t_n$ and titre type $k$ , using a piecewise linear function to capture boosting and waning phases:

Boosting Phase: For $t_n < t_{p,n}$ , the titre increases at rate $m_{1,n}$ .
Plateau Phase: For $t_{p,n} \leq t_n \leq t_{s,n}$ , the titre remains elevated.
Waning Phase: For $t_n > t_{s,n}$ , the titre decreases at rate $m_{3,n}$ .

Mathematical Formulation

Expected Titre Value

The expected log-transformed titre value $\mu_{n}$ is given by:

$\mu_{n} = t_{0,n} + \begin{cases} m_{1,n} \cdot t_n, & \text{if } t_n < t_{p,n} \\ m_{1,n} \cdot t_{p,n} + m_{2,n} \cdot (t_n - t_{p,n}), & \text{if } t_{p,n} \leq t_n \leq t_{s,n} \\ m_{1,n} \cdot t_{p,n} + m_{2,n} \cdot (t_{s,n} - t_{p,n}) + m_{3,n} \cdot (t_n - t_{s,n}), & \text{if } t_n > t_{s,n} \end{cases}$

where:

$t_{0,n}$ : Baseline titre level for individual $n$ and titre type $k$ .
$t_{p,n}$ : Time to peak titre for individual $n$ and titre type $k$ .
$t_{s,n}$ : Time to start of waning for individual $n$ and titre type $k$ .
$m_{1,n}$ : Boosting rate for individual $n$ and titre type $k$ .
$m_{2,n}$ : Plateau rate (assumed to be zero in this model).
$m_{3,n}$ : Waning rate for individual $n$ and titre type $k$ .
$t_n$ : Observation time for individual $n$ .
$\mu_{n} \geq 0$ : Ensured by taking the maximum with zero.

Observation Model

The observed log-transformed titre values $y_n$ are modeled as:

$y_n \sim \text{Normal}(\mu_{n}, \sigma)$

where $\sigma$ is the measurement error standard deviation.

Censoring

The model accounts for left-censoring and right-censoring:

Left-Censoring: For observations below detection limit $L$ , the likelihood contribution is:

$P(y_n \leq L) = \Phi\left(\dfrac{L - \mu_{n}}{\sigma}\right)$
Right-Censoring: For observations above detection limit $U$ , the likelihood contribution is:

$P(y_n \geq U) = 1 - \Phi\left(\dfrac{U - \mu_{n}}{\sigma}\right)$

where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution.

Hierarchical Structure

Individual-Level Parameters

For each individual $n$ and titre type $k$ , the parameters are modeled as:

$\begin{aligned} t_{0,n} &= t_{0,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{t_0} + \sigma_{t_0,k} \cdot z_{t_0,n} \\ t_{p,n} &= t_{p,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{t_p} + \sigma_{t_p,k} \cdot z_{t_p,n} \\ t_{s,n} &= t_{p,n} + \Delta t_{s,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{t_s} + \sigma_{t_s,k} \cdot z_{t_s,n} \\ m_{1,n} &= m_{1,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{m_1} + \sigma_{m_1,k} \cdot z_{m_1,n} \\ m_{2,n} &= m_{2,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{m_2} + \sigma_{m_2,k} \cdot z_{m_2,n} \\ m_{3,n} &= m_{3,k} + \mathbf{x}_n^\top \boldsymbol{\beta}_{m_3} + \sigma_{m_3,k} \cdot z_{m_3,n} \end{aligned}$

where:

$\mathbf{x}_n$ : Covariate vector for individual $n$ .
$\boldsymbol{\beta}_{\cdot}$ : Regression coefficients for the corresponding parameter.
$\sigma_{\cdot,k}$ : Standard deviation of individual-level random effects for titre type $k$ .
$z_{\cdot,n} \sim \text{Normal}(0, 1)$ : Standard normal random variables.

Population-Level Parameters

The population-level parameters for each titre type $k$ have the following priors:

$\begin{aligned} T_0^{p,k} &\sim \text{Normal}(\mu_{t_0}, \sigma_{t_0}) \\ t_p^{p,k} &\sim \text{Normal}(\mu_{t_p}, \sigma_{t_p}) \\ \Delta t_{p,k} &\sim \text{Normal}(\mu_{t_s} - \mu_{t_p}, \sigma_{t_s}) \\ m_1^{p,k} &\sim \text{Normal}(\mu_{m_1}, \sigma_{m_1}) \\ m_2^{p,k} &\sim \text{Normal}(\mu_{m_2}, \sigma_{m_2}) \\ m_3^{p,k} &\sim \text{Normal}(\mu_{m_3}, \sigma_{m_3}) \end{aligned}$

The standard deviations of the individual-level random effects have priors:

$\sigma_{k} \sim \text{Normal}(0, \sigma_{p})$

Regression Coefficients

The regression coefficients have the following priors:

$\boldsymbol{\beta}_{\cdot} \sim \text{Normal}(\mathbf{0}, \sigma_{\beta_{\cdot}})$

with appropriate constraints for parameters that must be positive or negative.

Data and Parameters

Data Inputs

$N$ : Total number of observations.
$K$ : Number of titre types.
$t_n$ : Observation times.
$y_n$ : Observed log-transformed titre values.
$\mathbf{x}_n$ : Covariate vector for each individual.
Censoring indicators for left and right censoring.

Parameters to Estimate

Population-level parameters: $t_{0,k}$ , $t_{p,k}$ , $\Delta t_{s,k}$ , $m_{1,k}$ , $m_{2,k}$ , $m_{3,k}$ .
Individual-level random effects: $z_{\cdot,n}$ .
Regression coefficients: $\boldsymbol{\beta}_{\cdot}$ .
Measurement error standard deviation: $\sigma$ .

Priors

The prior distributions are specified based on previous studies and domain knowledge:

Population Means: Specified using normal distributions with means $\mu_{\cdot}$ and standard deviations $\sigma_{\cdot}$ .
Random Effect Standard Deviations: Weakly informative normal priors centered at zero.
Regression Coefficients: Weakly informative normal priors centered at zero.
Measurement Error: $\sigma \sim \text{Normal}(0, 2)$ , constrained to be positive.

Likelihood

The likelihood function combines the observation model with the censoring mechanisms:

$\begin{aligned} &\text{For uncensored observations:} \quad y_n \sim \text{Normal}(\mu_{n}, \sigma) \\ &\text{For left-censored observations:} \quad P(y_n \leq L) = \Phi\left(\dfrac{L - \mu_{n}}{\sigma}\right) \\ &\text{For right-censored observations:} \quad P(y_n \geq U) = 1 - \Phi\left(\dfrac{U - \mu_{n}}{\sigma}\right) \end{aligned}$