Modelling clinical score in the TNA study
Model variation with individual effect of goat 12
## 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.
## Warning: Method 'posterior_samples' is deprecated. Please
## see ?as_draws for recommended alternatives.1 Introduction
In this variation we investigate the influence of goat number 12 in the results by introducing a specific intercept for it in the model.
Compare results with the previous model.
2 Statistical model
Same model as before with an additional indicator variable in the linear predictor for the individual 12.
3 Results
Figure 3.1: Posterior means by group. Medians and 50 % credible bands.
Figure 3.2: Predicted outcome by group. Medians and 50 % credible bands.
Figure 3.3: Average effect of the test by day.
Figure 3.4: Posterior probability that the average effect of the test is a decrease in clinical score, by day.
Figure 3.5: Predictive effect of the test by day.
Figure 3.6: Posterior probability that the test decreases the clinical score in an individual, by day.
Figure 3.7: Posterior effect of the boxes, relative to the group means, in the latent scale.
Figure 3.8: Difference between box 2 and box 5 in the test group. The posterior mean difference is of about 0.2 score units, and the probability of a positive difference is of about 98 %
Figure 3.9: Posterior distribution of the effect of the goat number 12.
Figure 3.10: Conditional effect of the boxes and the time.
3.1 Model comparison criteria
## Output of model 'tna_cs_fit':
##
## Computed from 800 by 464 log-likelihood matrix
##
## Estimate SE
## elpd_loo -657.2 23.5
## p_loo 11.8 1.7
## looic 1314.4 47.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 463 99.8% 245
## (0.5, 0.7] (ok) 1 0.2% 209
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'tna_cs_fit_id12':
##
## Computed from 800 by 464 log-likelihood matrix
##
## Estimate SE
## elpd_loo -651.4 22.8
## p_loo 12.4 1.7
## looic 1302.8 45.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## tna_cs_fit_id12 0.0 0.0
## tna_cs_fit -5.8 5.24 Conclusions
The model indeed fits the data better, which was expected since it has an additional parameter just to account for some previously unexplained variation.
However, Figure ?? shows that individual 12 explains only part of the effect of box 2. And while the outcome of goat 12 is clearly different from the rest (Fig. 3.9), the magnitude of the effect is small.
Unless some evidence of measurement of manipulation error is found, there is not sufficient support for this model in favour of the previous one.
5 Annex: model diagnostics and validation
## Family: poisson
## Links: mu = log
## Formula: sc ~ grp + (1 | box) + s(t, k = 5, by = grp, id = "grp") + id12
## Data: tna_data %>% mutate(id12 = as.numeric(id_goat == " (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
## sds(stgrpcontrol_1) 4.07 1.53 1.87 7.99
## sds(stgrptest_1) 3.70 1.54 1.60 7.50
## Rhat Bulk_ESS Tail_ESS
## sds(stgrpcontrol_1) 1.01 329 445
## sds(stgrptest_1) 1.01 359 406
##
## Group-Level Effects:
## ~box (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.24 0.20 0.04 0.76 1.03
## Bulk_ESS Tail_ESS
## sd(Intercept) 271 285
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.23 0.22 -0.23 0.68 1.02
## grptest -0.29 0.28 -0.83 0.23 1.01
## id12 0.65 0.16 0.31 0.95 1.00
## st:grpcontrol_1 1.13 1.59 -1.96 4.20 1.01
## st:grptest_1 0.59 1.55 -2.40 3.73 1.00
## Bulk_ESS Tail_ESS
## Intercept 421 405
## grptest 383 292
## id12 776 545
## st:grpcontrol_1 760 527
## st:grptest_1 985 570
##
## 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.
Figure 5.1: Observed vs posterior predictive modes, and highest-density credible intervals at 50% and 89%.
## Using 10 posterior draws for ppc type 'dens_overlay_grouped' by default.
Figure 5.2: Posterior predictive simulations from the model vs. observed empirical distribution.
## Warning: Argument 'nsamples' is deprecated. Please use
## argument 'ndraws' instead.
Figure 5.3: Posterior predictive number of zeros in the outcome by day and group.
Figure 5.4: Horizontally-jittered residuals (mode, 50 and 89% HDI credible intervals) agains time.
Figure 5.5: Quantile residuals. Is this appropriate for a discrete variable?