Modelling clinical score in the TNA study
## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## Use [add_]epred_draws() to get the expectation of the posterior predictive.
## Use [add_]linpred_draws() to get the distribution of the linear predictor.
## For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## Use [add_]epred_draws() to get the expectation of the posterior predictive.
## Use [add_]linpred_draws() to get the distribution of the linear predictor.
## For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
1 Introduction
Trial for a therapeutic vaccine (i.e. to treat the disease, rather than prevent it). Inserted the treatment in attenuated virus from classic PPR vaccine made resistant. This virus will hopefully transport the treatment to infected cells with the real (wild) virus.
25-30 individual animals (goats) split into 5 boxes
Boxes 1, 2, 5: treatment group with vaccine administered at days -1, 0, 1 with respect to the challenge (infection with real virus).
Boxes 6, 8: negative control group, with classical PPR vaccine (not expected to have any therapeutic effect) and no treatment at all.
Response variable:
- Clinical score (SC): categorical variable from 1 (no symptoms) to 5 (death)
Measurements:
- Each animal is evaluated for both response variables almost every day from 5 to 14. There are some missing values.
2 Statistical model
Let \(y(t) | b\) be the observed clinical score at day \(t\) for a individual in box \(b \in \{1, 2, 5, 6, 8\}\) from group \(g\) (either control or test). We assume a Poisson model, whose log-expected value is described additively by four components:
\(\mu_0\), a constant baseline mean value for the control group;
\(\alpha_t\), the mean difference of test group, relative to the control;
\(\theta_b\), the mean difference of box \(b\), relative to the corresponding group;
\(f_g(t)\), the expected evolution over time of the group \(g\) centred at 0. The non-parametric curves \(f_g(t)\) were implemented using thin-plate regression splines (Wood 2017).
The model equations are as follows:
\[\begin{equation} \begin{aligned} y(t) | b & \sim \mathcal{P}(\mu_b(t)) \\ \log(\mu_b(t)) & = \mu_0 + \alpha_t \mathbb{1}_{\{g(b) = \text{test}\}} + \theta_b + f_g(t) \\ \theta_b & \sim \mathcal{N}(0,\, \sigma_\theta) \\ f_g(t) & = \beta_g \phi(t) + \sum_{k=1}^K \gamma^{(g)}_k B_k(i) \\ \gamma_k & \sim \mathcal{N}(0, \sigma_\gamma), \quad k = 1, \ldots, K \end{aligned} \label{eq:model} \end{equation}\] where \(B_k\) are the thin-plate regression spline basis functions of the location \(i\) for the current data, with \(K = 6\) and penalty order \(2\) and \(\phi\) is a centred linear function of \(t\).
The model was fit using brms
(Bürkner 2017).
3 Results
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `stat(x)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(x)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Day | Median effect | Lower 80% CrI | Upper 80% CrI |
---|---|---|---|
0 | 0.01 | -0.10 | 0.11 |
1 | 0.01 | -0.10 | 0.12 |
2 | 0.01 | -0.12 | 0.14 |
3 | 0.01 | -0.15 | 0.19 |
4 | -0.02 | -0.22 | 0.21 |
5 | -0.08 | -0.37 | 0.23 |
6 | -0.23 | -0.69 | 0.19 |
7 | -0.46 | -1.13 | 0.31 |
8 | -0.80 | -1.79 | 0.35 |
9 | -1.19 | -2.79 | 0.18 |
10 | -1.50 | -3.20 | 0.41 |
11 | -1.66 | -3.52 | 0.47 |
12 | -1.69 | -3.30 | 0.75 |
13 | -1.61 | -3.33 | 0.34 |
14 | -1.49 | -2.98 | 0.27 |
15 | -1.41 | -2.79 | 0.20 |
4 Conclusions
The model works quite well. It proved robust with respect to changes in the prior specifications and to the use of an over-dispersed likelihood like the negative binomial.
The test group clearly performs better on average, with a posterior probability of about 90 % for a beneficial effect starting from day 8 (Fig. 3.4), and a estimated (posterior mean) reduction of up to about 2.1 clinical score units on day 11 (Fig. 3.3.
The benefit for an individual is less certain, due to the relatively high individual variability (Fig. 3.6.
In particular during the first week, were clinical scores are very low even in the control group, is quite unlikely that the test improves the outcome.
5 Annex: model diagnostics and validation
## Family: poisson
## Links: mu = log
## Formula: sc ~ 1 + grp + (1 | box) + +s(t, k = 5, by = grp, id = "grp")
## Data: tna_data (Number of observations: 464)
## Draws: 4 chains, each with iter = 400; warmup = 200; thin = 1;
## total post-warmup draws = 800
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(stgrpcontrol_1) 3.99 1.47 1.85 7.42 1.01 512 625
## sds(stgrptest_1) 3.86 1.46 1.90 7.36 1.00 412 462
##
## Group-Level Effects:
## ~box (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.39 0.11 1.53 1.00 236 281
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.22 0.42 -0.53 1.14 1.01 449 225
## grptest -0.26 0.47 -1.19 0.66 1.01 368 359
## st:grpcontrol_1 1.12 1.60 -2.01 4.35 1.01 1084 395
## st:grptest_1 0.50 1.57 -2.55 3.63 1.00 1117 660
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.
## Using 10 posterior draws for ppc type 'dens_overlay_grouped' by default.
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
## Warning:
## 4 (0.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.
## Warning:
## 4 (0.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.
## Output of model 'tna_cs_fit':
##
## Computed from 800 by 464 log-likelihood matrix
##
## Estimate SE
## elpd_waic -657.1 23.5
## p_waic 11.8 1.7
## waic 1314.3 47.0
##
## 4 (0.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.
##
## Output of model 'tna_cs_fit_negbinomial':
##
## Computed from 800 by 464 log-likelihood matrix
##
## Estimate SE
## elpd_waic -656.0 22.7
## p_waic 12.3 2.6
## waic 1312.1 45.4
##
## 4 (0.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.
##
## Model comparisons:
## elpd_diff se_diff
## tna_cs_fit_negbinomial 0.0 0.0
## tna_cs_fit -1.1 1.8