Examples of Frequentist vs. Bayesian Longitudinal Proportional Odds Models
Nathan James
nathan.t.james@vanderbilt.edu
2020-03-31
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.
Schizophrenia Example from 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 ' ' 1Now 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).
# random intercept and random slope for SqrtWeek
#! takes a few minutes to compile & run
brms.fit2 <- brm(imps79o ~ TxDrug + SqrtWeek + TxSWeek + (1+SqrtWeek|id),
data = schizophrenia, family=cumulative("logit"))plot(brms.fit2)
summary(brms.fit2) Family: cumulative
Links: mu = logit; disc = identity
Formula: imps79o ~ TxDrug + SqrtWeek + TxSWeek + (1 + SqrtWeek | 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
sd(Intercept) 2.70 0.25 2.22 3.22 1.00 932
sd(SqrtWeek) 1.46 0.15 1.18 1.77 1.01 643
cor(Intercept,SqrtWeek) -0.40 0.09 -0.55 -0.22 1.01 933
Tail_ESS
sd(Intercept) 1693
sd(SqrtWeek) 1074
cor(Intercept,SqrtWeek) 2137
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1] -7.41 0.48 -8.39 -6.49 1.00 1369 2470
Intercept[2] -3.46 0.39 -4.25 -2.70 1.00 1696 2343
Intercept[3] -0.83 0.36 -1.52 -0.14 1.00 1970 2789
TxDrug 0.06 0.40 -0.73 0.82 1.00 2043 2487
SqrtWeek -0.89 0.22 -1.35 -0.45 1.00 1975 2265
TxSWeek -1.73 0.26 -2.26 -1.22 1.00 1737 2323
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).
Computing Environment
R version 3.6.3 (2020-02-29) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Pop!_OS 19.10 Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] brms_2.12.0 Rcpp_1.0.4 mixor_1.0.4 Hmisc_4.4-0 [5] ggplot2_3.3.0 Formula_1.2-3 survival_3.1-11 lattice_0.20-40To cite R in publication use:
R Core Team (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.