I am conducting a mutlivariate regression model in brms, modeling the effect of ampehtamine use at start of treatment on three outcomes - psychological health, physical health and quality of life - over the first year of treatment. These outcomes three outcomes are all modelled on a 0-10 scale where higher scores indicate better health. My goal is to compare a Gaussian version of the model to an ordinal version. Both models use the same outcome data. To enable comparison we add 1 to all scores, making it a 1-11 rather than a 0-10 scale (see here).
The Gaussian model is belowats_baseline is number of days of amphetamine use in the 28 days prior to start of treatment. yearsFromStart indicates the time in years each measurement was made. Note the set_rescor(FALSE)argument here, preventing correlations between residuals of outcomes to be modelled.
bform_gaussian_mvr_ppq_rF <- bf(mvbind(psych, phys, qol) ~ ats_baseline*yearsFromStart + (1|p|encID)) + set_rescor(FALSE)
The same model with set_rescor(TRUE)(i.e. allowing correlations between residuals of outcomes) looks like this
bform_gaussian_mvr_ppq <- bf(mvbind(psych, phys, qol) ~ ats_baseline*yearsFromStart + (1|p|encID)) + set_rescor(TRUE)
And this is the linear ordinal
bform_ord_mvr_ppq <- bf(mvbind(psych, phys, qol) ~ ats_baseline*yearsFromStart + (1|p|encID))
I performed loo-cv on each model and then compared them using the loo_compare()function
These are the results. The ordinal is far better than the Gaussian without correlated residuals. (Apologies but the fit names are different to the bf()function names. I didn’t include the whole brms model specification to save space).
elpd_diff se_diff
fit_ordinal_mvr_ppq_ao_nm 0.0 0.0
fit_mvr_lin_ppq_ao_nm_rF -279.6 30.1
The ordinal model clearly fits the data better than the Gaussian sans correlations.
But when we compare it with the Gaussian model with correlations
elpd_diff se_diff
fit_mvr_lin_ppq_ao_nm 0.0 0.0
fit_ordinal_mvr_ppq_ao_nm -478.3 59.2
The results have flipped: Gaussian now clearly outperforms the ordinal.
My question is is this sort of behaviour expected or unusual? Does modelling the residual correlations typically make such a big difference
For completeness below is the output from the Gaussian with residual correlations
Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity
mu = identity
mu = identity
Formula: psych ~ ats_baseline * yearsFromStart + (1 | p | encID)
phys ~ ats_baseline * yearsFromStart + (1 | p | encID)
qol ~ ats_baseline * yearsFromStart + (1 | p | encID)
Data: mvrppq (Number of observations: 1906)
Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
total post-warmup draws = 8000
Multilevel Hyperparameters:
~encID (Number of levels: 598)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(psych_Intercept) 1.10 0.05 1.00 1.21 1.00 3119
sd(phys_Intercept) 1.08 0.05 0.98 1.19 1.00 2054
sd(qol_Intercept) 1.14 0.06 1.03 1.25 1.00 2903
cor(psych_Intercept,phys_Intercept) 0.67 0.04 0.58 0.75 1.00 2257
cor(psych_Intercept,qol_Intercept) 0.81 0.03 0.74 0.86 1.00 1989
cor(phys_Intercept,qol_Intercept) 0.66 0.04 0.58 0.74 1.00 2774
Tail_ESS
sd(psych_Intercept) 5083
sd(phys_Intercept) 4017
sd(qol_Intercept) 4011
cor(psych_Intercept,phys_Intercept) 4542
cor(psych_Intercept,qol_Intercept) 4300
cor(phys_Intercept,qol_Intercept) 4540
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
psych_Intercept 7.10 0.07 6.96 7.25 1.00 6509
phys_Intercept 7.27 0.07 7.13 7.41 1.00 6712
qol_Intercept 7.28 0.07 7.14 7.43 1.00 6420
psych_ats_baseline -0.06 0.01 -0.08 -0.03 1.00 5971
psych_yearsFromStart 0.74 0.14 0.47 1.01 1.00 9366
psych_ats_baseline:yearsFromStart 0.08 0.03 0.03 0.14 1.00 8399
phys_ats_baseline -0.03 0.01 -0.06 -0.00 1.00 6368
phys_yearsFromStart 0.58 0.13 0.32 0.83 1.00 9716
phys_ats_baseline:yearsFromStart 0.02 0.03 -0.03 0.08 1.00 8529
qol_ats_baseline -0.06 0.01 -0.09 -0.03 1.00 6516
qol_yearsFromStart 0.77 0.14 0.50 1.04 1.00 9024
qol_ats_baseline:yearsFromStart 0.08 0.03 0.03 0.14 1.00 8265
Tail_ESS
psych_Intercept 6354
phys_Intercept 6307
qol_Intercept 6362
psych_ats_baseline 6554
psych_yearsFromStart 6877
psych_ats_baseline:yearsFromStart 6710
phys_ats_baseline 6255
phys_yearsFromStart 6609
phys_ats_baseline:yearsFromStart 6818
qol_ats_baseline 6208
qol_yearsFromStart 6442
qol_ats_baseline:yearsFromStart 6469
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_psych 1.56 0.03 1.50 1.62 1.00 6359 5599
sigma_phys 1.50 0.03 1.44 1.56 1.00 5400 6244
sigma_qol 1.56 0.03 1.50 1.62 1.00 5713 5911
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(psych,phys) 0.51 0.02 0.47 0.55 1.00 5500 5697
rescor(psych,qol) 0.58 0.02 0.55 0.62 1.00 4810 5915
rescor(phys,qol) 0.54 0.02 0.51 0.58 1.00 5835 6318
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).
Worth noting brms does not support modelling residuals in ordinal models.