In this vignette, we discuss how to specify multilevel models with
compositional outcomes using multilevelcoda. In addition to
multilevelcoda, we will use brms package (to
fit models) and bayestestR package (to compute useful
indices and compare models). We will also attach built in datasets
mcompd (simulated compositional sleep and wake variables)
and sbp (sequential binary partition).
library(multilevelcoda)
library(brms)
library(bayestestR)
data("mcompd")
data("sbp")
options(digits = 3)The ILR coordinates outcomes can be calculated using the
complr() functions.
A model with multilevel compositional outcomes is multivariate, as it
has multiple ILR coordinate outcomes,each of which is predicted by a set
of predictors. Our brms model can be then fitted using the
brmcoda() function.
mv <- brmcoda(complr = cilr,
formula = mvbind(ilr1, ilr2, ilr3, ilr4) ~ Stress + (1 | ID),
cores = 8, seed = 123, backend = "cmdstanr")
#> Warning: In the future, 'rescor' will be set to FALSE by default for all models. It is thus
#> recommended to explicitely set 'rescor' via 'set_rescor' instead of using the default.Here is a summary() of the model. We can see that stress
significantly predicted ilr1 and ilr2.
summary(mv)
#> Family: MV(gaussian, gaussian, gaussian, gaussian)
#> Links: mu = identity; sigma = identity
#> mu = identity; sigma = identity
#> mu = identity; sigma = identity
#> mu = identity; sigma = identity
#> Formula: ilr1 ~ Stress + (1 | ID)
#> ilr2 ~ Stress + (1 | ID)
#> ilr3 ~ Stress + (1 | ID)
#> ilr4 ~ Stress + (1 | ID)
#> Data: tmp (Number of observations: 3540)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Multilevel Hyperparameters:
#> ~ID (Number of levels: 266)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(ilr1_Intercept) 0.33 0.02 0.30 0.37 1.00 1228 2060
#> sd(ilr2_Intercept) 0.30 0.01 0.28 0.33 1.00 1134 1929
#> sd(ilr3_Intercept) 0.39 0.02 0.35 0.43 1.00 1731 2763
#> sd(ilr4_Intercept) 0.30 0.02 0.27 0.33 1.00 1710 2483
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> ilr1_Intercept -0.44 0.02 -0.48 -0.39 1.00 942 1598
#> ilr2_Intercept 1.47 0.02 1.42 1.51 1.01 876 1608
#> ilr3_Intercept -0.88 0.03 -0.94 -0.82 1.00 1740 2565
#> ilr4_Intercept 0.65 0.02 0.60 0.69 1.00 1648 2562
#> ilr1_Stress -0.01 0.00 -0.02 -0.00 1.00 5967 3302
#> ilr2_Stress 0.01 0.00 0.00 0.01 1.00 5446 3533
#> ilr3_Stress 0.00 0.01 -0.01 0.01 1.00 6975 3101
#> ilr4_Stress 0.01 0.00 -0.00 0.01 1.00 6794 3387
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma_ilr1 0.44 0.01 0.43 0.45 1.00 6228 3088
#> sigma_ilr2 0.38 0.00 0.37 0.39 1.00 6496 3182
#> sigma_ilr3 0.70 0.01 0.68 0.71 1.00 5403 3385
#> sigma_ilr4 0.53 0.01 0.51 0.54 1.00 5484 3418
#>
#> Residual Correlations:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> rescor(ilr1,ilr2) -0.54 0.01 -0.57 -0.52 1.00 5503 3659
#> rescor(ilr1,ilr3) -0.18 0.02 -0.21 -0.14 1.00 5662 3335
#> rescor(ilr2,ilr3) -0.05 0.02 -0.08 -0.02 1.00 4974 3475
#> rescor(ilr1,ilr4) 0.11 0.02 0.07 0.14 1.00 6096 3480
#> rescor(ilr2,ilr4) -0.05 0.02 -0.08 -0.01 1.00 5544 3305
#> rescor(ilr3,ilr4) 0.56 0.01 0.54 0.58 1.00 5224 3252
#>
#> Draws were sampled using sample(hmc). 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).We are often interested in whether a predictor significantly predict the overall composition, in addition to the individual ILR coordinates. In Bayesian, this can be done by comparing the marginal likelihoods of two models. Bayes Factors