Insufficient Statistics - Neutralizing Antibody Titres

Reading the Pre-Prints

A prominant scientist has been posting his highlights of pre-prints on twitter during the course of the pandemic and occasionally, I come across one that makes me scratch my head. I think that the move to open and rapid science through pre-prints has been good, but not without issues. Review (and open review especially) makes work better and is an important part of science and progress.

One paper that stood out to me was by Suthar et al. (2021) in which the authors look at the neutralizing antibody tires (nAbs) amongst T-cell response in those persons who were vaccinated with the Pfizer SARS-CoV-2 Vaccines (BNT162b2 vaccine). The authors found: “substantial waning of antibody response and T-cell Immunity” in those who were vaccinated. Naturally, this could be concerning because antibodies and T-cells are important in fighting infection and could help us understand the endgame of this first epidemic of SARS-CoV-2 (e.g., vaccine induced reduction of transmission, symptomatic disease, severe outcomes, and associated periodicity given time since last vaccine). However, there is the open question of what is the best correlate of protection (e.g., we know that antibodies wane with time, but the ability to produce antibodies in the face of a new infection from memory B-cells is also very important) (Corbett et al. 2021)

So Look at the Data

As mentioned, along with the bad things about pre-prints (media jumping on sensationalist findings, people seeking acclaim by making bold statements that are not rigorously founded in fact, general review help, etc) is that much of the data are available for secondary analysis. With the article in question, I could find some of the original data analyzed by the authors from a prior publication. We can start to look at these data ourselves:

library(data.table)
library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

✓ ggplot2 3.3.5.9000     ✓ purrr   0.3.4.9000
✓ tibble  3.1.6.9001     ✓ dplyr   1.0.8.9000
✓ tidyr   1.1.4.9000     ✓ stringr 1.4.0.9000
✓ readr   2.1.2          ✓ forcats 0.5.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::between()   masks data.table::between()
x dplyr::filter()    masks stats::filter()
x dplyr::first()     masks data.table::first()
x dplyr::lag()       masks stats::lag()
x dplyr::last()      masks data.table::last()
x purrr::transpose() masks data.table::transpose()

theme_set(ggdist::theme_ggdist())
dat <- readxl::read_excel("nAbs.xlsx")

dat %>%
  gather(day, value, -PID, -Gender) %>%
  mutate(day = as.numeric(str_extract(string = day, "\\d+")))->dat_long

knitr::kable(dat, digits = 1)

PID	Gender	D_0	D_21	D_42	D_90_120
2000	Male	10.0	62.3	1464.8	666.0
2001	Male	10.0	10.0	964.3	273.0
2007	Female	10.0	10.0	195.3	101.5
2009	Female	10.0	10.0	224.0	155.7
2011	Female	10.0	51.0	679.9	322.5
2012	Male	10.0	73.5	343.4	173.5
2013	Female	10.0	10.0	229.1	145.1
2018	Female	10.0	12.3	243.6	218.8
2020	Female	10.0	10.0	73.2	88.0
2028	Female	10.0	40.8	814.5	326.0
2030	Female	10.0	47.0	286.9	155.1
2031	Female	10.0	10.0	166.9	68.8
2034	Male	10.0	10.0	474.1	237.0
2036	Male	10.0	42.9	306.4	NA
2040	Male	10.0	10.0	249.3	188.0
2042	Female	10.0	34.8	484.1	200.0
2043	Male	10.0	10.0	37.4	108.9
2044	Male	10.0	30.5	927.8	341.1
2045	Female	10.0	10.0	266.8	139.7
2046	Male	10.0	10.0	262.2	75.6
2047	Male	10.0	279.4	313.1	245.2
2048	Male	10.0	10.0	678.0	297.4
2049	Male	10.0	10.0	182.7	167.5
2050	Male	10.0	10.0	431.5	217.9
2051	Female	10.0	79.2	435.8	NA
2052	Male	10.0	10.0	252.1	151.3
2053	Female	10.0	27.5	319.3	184.0
2054	Female	10.0	35.1	179.2	117.1
2055	Female	10.0	28.2	493.1	342.0
2056	Male	10.0	33.1	387.3	466.0
2002	Female	10.0	77.0	801.1	617.0
2005	Female	10.0	19.4	349.3	NA
2006	Female	10.0	44.6	423.2	217.0
2008	Female	10.0	61.5	291.7	101.7
2010	Female	10.0	34.5	473.0	NA
2014	Female	10.0	10.0	15.8	NA
2016	Female	10.0	44.1	483.2	NA
2019	Female	10.0	10.0	112.0	NA
2022	Female	10.0	10.0	131.4	108.0
2023	Female	10.0	13.0	203.3	NA
2025	Female	58.9	856.3	2095.1	769.0
2026	Female	10.0	21.4	387.2	NA
2027	Female	10.0	10.0	NA	80.1
2003	Male	10.0	10.0	267.0	NA
2004	Male	10.0	12.6	1268.4	523.0
2015	Male	10.0	2876.5	4379.5	2492.0
2017	Male	10.0	122.9	587.6	546.5
2024	Male	10.0	10.0	10.0	NA
2029	Male	10.0	10.0	341.9	159.8
2032	Male	10.0	1143.4	1845.0	247.0
2033	Male	10.0	15.4	NA	NA
2035	Male	10.0	35.3	422.1	NA
2037	Male	10.0	26.9	NA	NA
2038	Male	10.0	45.9	391.5	NA
2039	Male	10.0	10.0	143.6	NA
2041	Male	47.8	2131.4	3044.5	NA
2021	Female	10.0	10.0	10.0	NA

So it appears that there are 57 unique participants. Some participants were lost to follow-up and the final date range is kind of wide (90-120 days).

One of the most important steps in any analysis is to graph the data to gain a general understanding of the trends and shape of the data.

dat_long %>%
  ggplot(aes(day, value, group = PID))+
  geom_line()+
  facet_wrap(~PID)+
  labs(
    title = "Antibody Titres by Participate"
  )

Warning: Removed 20 row(s) containing missing values (geom_path).

Antibody Titres Over Time by Participant

One things that stands out to me is participant 2015. This individual has nAbs that are much higher than the other participants.

dat_long %>% 
  mutate(tst = ifelse(PID==2015,"PID: 2015","All Others")) %>% 
  group_by(tst, day) %>% 
  summarise(mu = mean(value, na.rm=TRUE)) %>% 
  spread(day,mu) %>% 
  knitr::kable(digits = 1, caption = "Mean Value of nAbs by Group")

`summarise()` has grouped output by 'tst'. You can override using the `.groups`
argument.

Mean Value of nAbs by Group
tst	0	21	42	90
All Others	11.5	104.0	499.3	251.1
PID: 2015	10.0	2876.5	4379.5	2492.0

If this were me, I would go back and check those data just in case someone got a little too happy when titrating (been there). Regardless, accepting this level of variability in the individual, these kind of data lend themselves well to hierarchical modeling.

So What’s the Half Life?

We can model the half life of the nAbs in the participants using the multilevel model framework with lme4. I will subset the data to include only the levels after the peak antibody level and include a random effect for the participant (since we know that there is within-person variability and serial dependence).

library(lme4)

Loading required package: Matrix


Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

fit_lme <- lmer(log(value) ~ day + (1|PID), data = subset(dat_long, day >21))

(half_life <- -log(2)/confint(fit_lme)[4,])

Computing profile confidence intervals ...

   2.5 %   97.5 % 
46.31701 84.74999

So based on this model, the half-life for nAbs is between 46.3 and 84.7 days.

nab_out <- function(t,hl){
  2^(-t/hl)
}

sim_dat <- tibble(day = 1:365) %>% 
  mutate(lower_bound = nab_out(day, half_life[1]),
         upper_bound = nab_out(day, half_life[2]))
sim_dat %>%   
  ggplot(aes(x = day))+
  geom_ribbon(aes(ymin = lower_bound,
                  ymax = upper_bound), 
              fill = "orange", alpha = .5)+
  geom_hline(yintercept = .25, col = "red", lty = "dashed")+
  labs(
    title = "Simulated Decay of nAbs",
    y  = "Percent of Maximum nAbs",
    x = "Day After Peak nAbs"
  )

So it appears that in roughly 3-6 months after maximum antibody levels a given vaccinated person would have 1/4 of nAb titres.

days_until_perc <- function(target, hl){
  -log(target)*hl/log(2)
}

round(days_until_perc(.25, half_life)/30,1)

 2.5 % 97.5 % 
   3.1    5.6

The big question of course is first, does this matter and in what context (protection vs infection where high levels of nAbs make sense vs symptomatic disease where a trained immune system is all that is required (thinking memory B and T cells)) and what concentration is enough at each level.

Is 5% of peak is enough then we have between 6 months and a year, perhaps in time for a teaser from infection to boost nAbs.

round(days_until_perc(.05, half_life)/30,1)

 2.5 % 97.5 % 
   6.7   12.2

What about compartmental models?

This is not my forte, but it would seem that compartmental models would be a good way to approach this problem. We would have to fix some parameters due to the data available but one could imagine that a decent data generating process could look like (at least as far as I understand how mRNA vaccines work):

Inoculum is added to the blood stream via vaccination
A latent phase for vaccinated individual’s cells to produce the spike protein which would be a function of the volume of inoculum and the number of cells encountered. The inoculum also has some decay function.
A phase where the encoded protein is in the bloodstream and recognized by the immune system as a foreign body which ramps up the innate and humoral immune response.
Eventually the spike proteins would be cleared at some rate as the immune response increases

I started down this path and then realized that I would need to make some assumptions for some of the parameters (transition rates between compartments). Being that I do not have reasonable priors for these rates, I decided not to code up this mdoel.

Other Method

I read another paper on waning immunity on coronaviruses that had much more longitudinal data (Townsend et al. 2021). The key model in the paper was locked away in a proprietary Mathematica Notebook that couldn’t open without buying the software. I coded the below model based on their methods, but we will find is that we do not have enough data to fit a very good model.

rw("decay.stan")

data {
    int<lower=1> N; // Number of Samples
    int<lower=1> J; //Number of participants
    vector [N] nAbs; //The concentration of the antibodies
    vector [N] time; //The concentration of the antibodies
}

parameters{
    real omega;
    real hl;
    real<lower=0> sigma;
}

transformed parameters {
    vector[N] theta;  //vac effect
    theta[1:N] = omega + (1- omega)  * exp(-hl*time[1:N]);
}

model {
    hl ~ normal(0,2);
    omega ~ normal(0,1);
    nAbs ~ normal(theta, sigma);
}

library(cmdstanr)

This is cmdstanr version 0.4.0

- Online documentation and vignettes at mc-stan.org/cmdstanr

- CmdStan path set to: /Users/michael/.cmdstanr/cmdstan-2.28.2

- Use set_cmdstan_path() to change the path


A newer version of CmdStan is available. See ?install_cmdstan() to install it.
To disable this check set option or environment variable CMDSTANR_NO_VER_CHECK=TRUE.

mod_ex <- cmdstan_model("decay.stan")

dat_stan_base <- dat_long %>% 
  filter(!is.na(value)) %>% 
  filter(day>21)

stan_dat <- list(
  N = nrow(dat_stan_base),
  J = length(unique(dat_stan_base$PID)),
  nAbs = dat_stan_base$value,
  time = dat_stan_base$day
)

fit_ex <- mod_ex$sample(stan_dat, 
                        refresh = 0, iter_sampling = 500,
                                                iter_warmup = 200,
                        adapt_delta = .98,
                        max_treedepth = 15)

Running MCMC with 4 parallel chains...

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: normal_lpdf: Location parameter[55] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3

Chain 2 finished in 0.2 seconds.
Chain 4 finished in 0.1 seconds.

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[55] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/Rtmp8xeEVQ/model-4cb07f1e8483.stan', line 22, column 1 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 3 finished in 0.3 seconds.
Chain 1 finished in 42.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 10.7 seconds.
Total execution time: 42.7 seconds.

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'


Warning: 232 of 2000 (12.0%) transitions ended with a divergence.
This may indicate insufficient exploration of the posterior distribution.
Possible remedies include: 
  * Increasing adapt_delta closer to 1 (default is 0.8) 
  * Reparameterizing the model (e.g. using a non-centered parameterization)
  * Using informative or weakly informative prior distributions

187 of 2000 (9.0%) transitions hit the maximum treedepth limit of 15 or 2^15-1 leapfrog steps.
Trajectories that are prematurely terminated due to this limit will result in slow exploration.
Increasing the max_treedepth limit can avoid this at the expense of more computation.
If increasing max_treedepth does not remove warnings, try to reparameterize the model.

Not a great fit as we are trying to draw a complicated curve through basically two points and there are (infinitely) many ways to do this.

fit_ex$summary(variables = c("omega", "hl", "sigma"))

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

# A tibble: 3 × 10
  variable    mean  median      sd     mad     q5    q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 omega     -0.523 -0.851    0.937   0.851 -1.43    1.20  1.72     7.97    34.0 
2 hl         0.139 -0.0548   1.52    1.43  -1.81    2.98  2.49     4.90    11.5 
3 sigma    389.    71.6    401.    105.     0.547 878.    2.37     4.98     7.61

Viral Titre?

A model that is somewhat of a compromise between the compartmental model which captures more of the underlying biological process and the hierarchical model which fit the within-host variability is to model the viral decay with a direct functional form. This is similar to the approach performed by Singanayagam et al. (2021) who use this functional form for estimated RNA copies from Cycle Threshold values. The issues, as before, is that we only have a few observations per person and I don’t have strong priors, but this model is closer to what we see in the graphs where the titre increases to some point and then decay afterwards. Ideally we would have some additional measures for the participants.

This can be modeled by the following equations:

\[ v(\tau) = v_{max} \frac{(a+b)}{be^{-a(\tau-\tau_{max})}+ae^{b(\tau-\tau_{max})}} \]

The a and b parameters can be fit and influence the rate of increase and decrease of nAb tire (or virus titre for that matter).

nab_kinetics <- function(v_max, a,b,t,tmax){
  v_max * (a+b)/(b*exp(-a*(t-tmax))+a*exp(b*(t-tmax)))
}

plot(nab_kinetics(50, a = .8, b= .1, t = 1:70, tmax = 20), type = "l",
     ylab = "nAb Titre", main = "Simulated Results")

We need to comport our data to provide the ability to identify the individual maximums and times for each participant.

long_model_dat <- dat_long %>% 
    filter(!is.na(value)) %>% 
    group_by(PID) %>% 
    mutate(v_max = max(value)) %>% 
    mutate(t_max =  max((v_max == value)*day))%>% 
    ungroup() %>% 
    mutate(id = as.integer(as.factor(PID)))

And our model:

rw("rise.stan")

Warning in readLines(x): incomplete final line found on 'rise.stan'

functions{
    real viral_kinetics (real v_max, real a, real b, real t, real t_max){
        return(v_max * (a + b) / (b * exp(-a*(t - t_max)) + a * exp(b*(t-t_max))));
    }
}

data {
    int<lower=1> N; // Number of records
    int<lower=0> id[N]; //Participant IDs
    vector<lower=0> [N] nAbs_max; // Maximum Ab concentration
    vector<lower=0> [N] nAbs; //The concentration of the antibodies
    vector [N] time; //The concentration of the antibodies
    vector [N] t_max; //The concentration of the antibodies
    real prior_a_mu;
    real prior_b_mu;
    real prior_a_sd;
    real prior_b_sd;
    int pred_window;
}

parameters {
    real<lower=0> a;
    real<lower=0> b;
    real<lower=0> sigma;
}

transformed parameters{
    vector<lower=0>[N] mu;
    
    for( i in 1:N){
        mu[i] = viral_kinetics(nAbs_max[id[i]], a, b, time[id[i]], t_max[id[i]]);
    }
    
}
model {
    a ~ normal(prior_a_mu, prior_a_sd);
    b ~ normal(prior_b_mu, prior_b_sd);
    sigma ~ normal(0,1);
    
    nAbs ~ normal(mu, sigma);
}

generated quantities{
    vector[pred_window] yhat;
    vector[pred_window] mu_pred;
    real avg_nabs_max = mean(nAbs_max);
    real avg_tmax = mean(t_max);
    
    for(i in 1:pred_window){
        mu_pred[i] = viral_kinetics(avg_nabs_max, a, b, i, avg_tmax);
        yhat[i] = normal_rng(mu_pred[i], sigma);
    }
    
}

We can fit the model in the usual ways:

mod_titre <- cmdstan_model("rise.stan")

titre_dat <- list(
    N = nrow(long_model_dat),
    id = long_model_dat$id,
    nAbs_max = long_model_dat$v_max,
nAbs = long_model_dat$value,
    time= long_model_dat$day,
    t_max = long_model_dat$t_max,
    prior_a_mu = .8,
    prior_b_mu = .1,
    prior_a_sd = .05,
    prior_b_sd = .05,
pred_window = 180
)

set.seed(336)

fit_titre <- mod_titre$sample(data = titre_dat, 
                                                            refresh = 0, adapt_delta = .99,
                                                            max_treedepth = 15)

Running MCMC with 4 parallel chains...

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: rise_model_namespace::log_prob: mu[sym1__] is nan, but must be greater than or equal to 0 (in '/var/folders/0x/6bnjy4n15kz8nbfbbwbyk3pr0000gn/T/RtmpJ4OPS5/model-dbcf188d05cc.stan', line 28, column 1 to column 23)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4

Chain 1 finished in 2.2 seconds.
Chain 2 finished in 2.3 seconds.
Chain 3 finished in 2.3 seconds.
Chain 4 finished in 2.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.2 seconds.
Total execution time: 2.6 seconds.

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

fit_titre$summary(variables = c("a", "b", "sigma"))

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

# A tibble: 3 × 10
  variable   mean  median     sd    mad      q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>   <dbl>  <dbl>  <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 a         0.801  0.801  0.0511 0.0501  0.714   0.885  1.00    2165.    1948.
2 b         0.102  0.0996 0.0479 0.0505  0.0252  0.182  1.00    1447.    1034.
3 sigma    91.9   91.9    0.499  0.494  91.1    92.7    1.01    2466.    2321.

Note that there is a fair bit of variance in our model represented by zero. Additionally, we have some values of the predicted titre below zero. We could resolve this by changing the distribution to use a distribution that is always positive like a lognormal distribution. This would require some careful transformation of our data and associated parameters. We’ll stick with this formulation for now.

We can now visualize our results from this model.

library(posterior)

This is posterior version 1.2.1


Attaching package: 'posterior'

The following objects are masked from 'package:stats':

    mad, sd, var

library(ggdist)
 pred_nab <- tidybayes::gather_draws(fit_titre$draws(variables = "yhat"),yhat[i] ) %>%
    dplyr::ungroup() %>%
    dplyr::group_by(i) %>% 
    ggdist::median_qi(.value, .width = c(.5, .8, .95))

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(from = 1, len = along - 1): partial argument match of
'len' to 'length.out'

Warning in seq.default(to = N - 1, len = N - along): partial argument match of
'len' to 'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.list)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = ncol(arg.dim)): partial argument match of 'len' to
'length.out'

Warning in seq.default(len = N): partial argument match of 'len' to 'length.out'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(along = arg.names): partial argument match of 'along' to
'along.with'

Warning in seq.default(len = length(arg.names)): partial argument match of 'len'
to 'length.out'

Warning in seq.default(along = perm): partial argument match of 'along' to
'along.with'

pred_nab %>% 
    ggplot(aes(i, .value))+
    geom_lineribbon(aes(ymin = .lower, ymax = .upper))+
    scale_fill_brewer()+
    labs(
        title = "Estimated Quantity of nAbs with Viral Evolution Model",
        y = "Quantity",
        x = "Days",
        fill = "Proportion of Scenarios"
    )+
    geom_hline(yintercept = 0, col = "orange", lty = "dashed")+
    theme(legend.position = "top")

Again, we see evidence of detectable tire above 180 days, but there is a decent bit of variability.

Conclusions

If we had additional measurements per person and some additional information on the kinetics, In the time it took me to write this post several other papers have come out with much more information that would be fun to analyze that find that immunity is likely long lasting (Mateus et al., n.d.; Grifoni et al. 2021).

References

Corbett, Kizzmekia S., Martha C. Nason, Britta Flach, Matthew Gagne, Sarah O’Connell, Timothy S. Johnston, Shruti N. Shah, et al. 2021. “Immune Correlates of Protection by mRNA-1273 Vaccine Against SARS-CoV-2 in Nonhuman Primates.” Science, July. https://doi.org/10.1126/science.abj0299.

Grifoni, Alba, John Sidney, Randi Vita, Bjoern Peters, Shane Crotty, Daniela Weiskopf, and Alessandro Sette. 2021. “SARS-CoV-2 Human t Cell Epitopes: Adaptive Immune Response Against COVID-19.” Cell Host & Microbe 29 (7): 1076–92. https://doi.org/10.1016/j.chom.2021.05.010.

Mateus, Jose, Jennifer M. Dan, Zeli Zhang, Carolyn Rydyznski Moderbacher, Marshall Lammers, Benjamin Goodwin, Alessandro Sette, Shane Crotty, and Daniela Weiskopf. n.d. “Low-Dose mRNA-1273 COVID-19 Vaccine Generates Durable Memory Enhanced by Cross-Reactive t Cells.” Science 0 (0): eabj9853. https://doi.org/10.1126/science.abj9853.

Singanayagam, Anika, Seran Hakki, Jake Dunning, Kieran J. Madon, Michael A. Crone, Aleksandra Koycheva, Nieves Derqui-Fernandez, et al. 2021. “Community Transmission and Viral Load Kinetics of the SARS-CoV-2 Delta (b.1.617.2) Variant in Vaccinated and Unvaccinated Individuals in the UK: A Prospective, Longitudinal, Cohort Study.” The Lancet Infectious Diseases 0 (0). https://doi.org/10.1016/S1473-3099(21)00648-4.

Suthar, Mehul S., Prabhu S. Arunachalam, Mengyun Hu, Noah Reis, Meera Trisal, Olivia Raeber, Sharon Chinthrajah, et al. 2021. “Durability of Immune Responses to the BNT162b2 mRNA Vaccine.” https://doi.org/10.1101/2021.09.30.462488.

Townsend, Jeffrey P., Hayley B. Hassler, Zheng Wang, Sayaka Miura, Jaiveer Singh, Sudhir Kumar, Nancy H. Ruddle, Alison P. Galvani, and Alex Dornburg. 2021. “The Durability of Immunity Against Reinfection by SARS-CoV-2: A Comparative Evolutionary Study.” The Lancet Microbe 0 (0). https://doi.org/10.1016/S2666-5247(21)00219-6.

Reuse

https://creativecommons.org/licenses/by/4.0/

Citation

BibTeX citation:

@online{dewitt2021,
  author = {Michael DeWitt},
  title = {Neutralizing {Antibody} {Titres}},
  date = {2021-11-11},
  url = {https://michaeldewittjr.com/programming/2021-10-11-neutralizing-antibody-titres},
  langid = {en}
}

For attribution, please cite this work as:

Michael DeWitt. 2021. “Neutralizing Antibody Titres.” November 11, 2021. https://michaeldewittjr.com/programming/2021-10-11-neutralizing-antibody-titres.