The R brms
package uses the same model syntax as the lme4
package so a basic random intercept ordinal model is fit with:
brm(outcome ~ treatment*time + (1|id), data=d, family=cumulative("logit"))
Here is a short script with an ordinal longitudinal model fit using both mixor
(frequentist) and brms
based on an example in the mixor
vignette.
For mixor see this and especially the package vignette.
mixor
Vignette## Using mixor
library(mixor)
data(schizophrenia)
# random intercept
mixor.fit1 <- mixor(imps79o ~ TxDrug*SqrtWeek, data = schizophrenia,
id = id, link = "logit")
summary(mixor.fit1)
Call:
mixor(formula = imps79o ~ TxDrug * SqrtWeek, data = schizophrenia,
id = id, link = "logit")
Deviance = 3402.758
Log-likelihood = -1701.379
RIDGEMAX = 0.2
AIC = -1708.379
SBC = -1722.659
Estimate Std. Error z value P(>|z|)
(Intercept) 5.85924 0.34288 17.0881 < 2.2e-16 ***
TxDrug -0.05843 0.31086 -0.1880 0.8509
SqrtWeek -0.76577 0.11975 -6.3950 1.606e-10 ***
TxDrug:SqrtWeek -1.20615 0.13314 -9.0595 < 2.2e-16 ***
Random.(Intercept) 3.77378 0.49543 7.6172 2.598e-14 ***
Threshold2 3.03282 0.13238 22.9103 < 2.2e-16 ***
Threshold3 5.15077 0.17925 28.7345 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# get cutpoints (instead of default departures from intercept)
cm <- matrix(c(-1, 0, 0, 0, 0, 0, 0,
-1, 0, 0, 0, 0, 1, 0,
-1, 0, 0, 0, 0, 0, 1), ncol = 3)
Contrasts(mixor.fit1, contrast.matrix = cm)
1 2 3
(Intercept) -1 -1 -1
TxDrug 0 0 0
SqrtWeek 0 0 0
TxDrug:SqrtWeek 0 0 0
Random.(Intercept) 0 0 0
Threshold2 0 1 0
Threshold3 0 0 1
Estimate Std. Error z value P(>|z|)
1 -5.85924 0.34288 -17.0881 < 2.2e-16 ***
2 -2.82642 0.29451 -9.5970 < 2.2e-16 ***
3 -0.70848 0.26989 -2.6251 0.008663 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# random intercept and random slope for SqrtWeek
mixor.fit2 <- mixor(imps79o ~ TxDrug + SqrtWeek + TxSWeek, data = schizophrenia,
id = id, which.random.slope = 2, link= "logit")
summary(mixor.fit2)
Call:
mixor(formula = imps79o ~ TxDrug + SqrtWeek + TxSWeek, data = schizophrenia,
id = id, which.random.slope = 2, link = "logit")
Deviance = 3325.486
Log-likelihood = -1662.743
RIDGEMAX = 0
AIC = -1671.743
SBC = -1690.103
Estimate Std. Error z value P(>|z|)
(Intercept) 7.318831 0.480778 15.2229 < 2.2e-16 ***
SqrtWeek -0.882261 0.234568 -3.7612 0.0001691 ***
TxDrug 0.057917 0.399102 0.1451 0.8846178
TxSWeek -1.694861 0.268131 -6.3210 2.598e-10 ***
(Intercept) (Intercept) 6.997646 1.369273 5.1105 3.213e-07 ***
(Intercept) SqrtWeek -1.508514 0.536023 -2.8143 0.0048888 **
SqrtWeek SqrtWeek 2.008916 0.453587 4.4290 9.469e-06 ***
Threshold2 3.901260 0.213257 18.2937 < 2.2e-16 ***
Threshold3 6.507172 0.289991 22.4392 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Now the corresponding Bayesian analysis, using brms
default priors.
## Using brms
library(brms)
options(mc.cores = parallel::detectCores()) # use multiple cores
# random intercept
#! takes a few minutes to compile & run
brms.fit1 <- brm(imps79o ~ TxDrug*SqrtWeek + (1|id),
data = schizophrenia, family=cumulative("logit"))
plot(brms.fit1) # check posterior densities and traceplots
summary(brms.fit1)
Family: cumulative
Links: mu = logit; disc = identity
Formula: imps79o ~ TxDrug * SqrtWeek + (1 | id)
Data: schizophrenia (Number of observations: 1603)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~id (Number of levels: 437)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.97 0.12 1.74 2.23 1.01 1120 1833
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1] -5.90 0.34 -6.59 -5.25 1.00 1416 1771
Intercept[2] -2.84 0.30 -3.45 -2.29 1.01 1419 1657
Intercept[3] -0.71 0.28 -1.26 -0.18 1.01 1368 1770
TxDrug -0.06 0.32 -0.68 0.56 1.00 1373 2010
SqrtWeek -0.77 0.13 -1.03 -0.51 1.00 2453 2486
TxDrug:SqrtWeek -1.21 0.15 -1.51 -0.92 1.00 2304 2541
Samples were drawn 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).