Set up

• setup print options
• parallel::detectCores() detect number of cores
• simulate a three predictors and one outcome (10 level ranks) data
• and run a frequentist logistic regression (ordinal)

require(rmsb)
options(prType='html')
## options(mc.cores = parallel::detectCores())   # use max # CPUs

set.seed(1)
n <- 500
x1 <- runif(n, -1, 1)
x2 <- runif(n, -1, 1)
x3 <- sample(0 : 1, n, TRUE)
y <- x1 + 0.5 * x2 + x3 + rnorm(n)
y <- as.integer(cut2(y, g=10))
dd <- datadist(x1, x2, x3); options(datadist='dd')
f <- lrm(y ~ x1 + pol(x2, 2) + x3, eps=1e-7) # eps to check against rstan
f
## table(y)

summary(f)

Flat Priors

• Flat normal priors for the betas
• default prior SD for blrm is 100;

for(psd in c(0.25, 1, 10, 100, 10000)) {
    cat('\nPrior SD:', psd, '\n')
    g <- blrm(y ~ x1 + pol(x2, 2) + x3, method='optimizing', priorsd=psd)
    cat('-2 log likelihood:', g$deviance, '\n')
    print(g$coefficients)
}

Bayesian Proportional Odds Ordinal Logistic Model

• Dirichlet priors on intercepts
• 0.95 highest posterior density intervals
• AUROC and \(R^2\) should be estimated with error
• The symmetry of a posterior distribution. The value of 1.0 indicates symmetry. The symmetry index is the ratio of distance from mean to 0.95 quantile and the distance from mean to 0.05 quantile.
• The proportional odds ordinal (PO) logistic model is a generalization of Wilcoxon/Kruskal-Wallis tests.

bs <- blrm(y ~ x1 + pol(x2, 2) + x3, file='bs.rds')
bs

Model Performance

• Stand diagnostic

blrmStats(bs, pl=TRUE)
stanDxplot(bs)
stanDx(bs)

Posterior Distribution

• The posterior distributions are calculated using kernel density estimates from posterior draws
• Using rstan::optimizing

plot(bs)
# Also show 2-d posterior density contour for two collinear terms
plot(bs, c('x2', 'x2^2'), bivar=TRUE)   # assumes ellipse
plot(bs, c('x2', 'x2^2'), bivar=TRUE, bivarmethod='kernel')   # kernel density

# Print frequentist side-by-side with Bayesian posterior mean, median, mode

cbind(MLE=coef(f), t(bs$param))

Contrasts

• Bayesian contrast’s point estimate is the posterior mean and the 0.95 posterior density interval
• instead of p-value, the posterior probability that the constrast is positive is computed

contrast(f, list(x1=0,x3=1),list(x1=.25,x3=0))
contrast(bs,list(x1=0:1,x3=1),list(x1=.25,x3=0))

Posterior Probability

P <- PostF(bs, pr=TRUE)   # show new short legal R names

P(b3 > 0 & b1 > 1.5)
P(b3 > 0)
P(abs(b3) < 0.25)  ## evidence for small nonlinearity

Contrained Partial PO

• Use the constrained partial proportional odds model to assess the proportional odds assumption.
• Assume departures from proportional odds (constant increments in log odds) are modeled as linear in square root of the outcome level.
• Group-stratified empirical CDFs to see visual evidence for this
• qlogis is the logit function logit(p)=log(p/(1-p))
• Relative explained variation (REV) is the Wald \(\chi^2\) statistics divided by the Wald statistics for the whole model.

## bcalp <- blrm(calpro ~ endo, ~ endo, cppo=sqrt)
## bcalp
## Ecdf(~ calpro, group=endo, fun=qlogis)

Missing Values

When possible, full joint Bayesian modeling of possible missing covariates and the outcome variable should be used to get exact inference in the presence of missing covariate values.

Then do posterior inference on the full stacked posterior distribution.

Imputation before model fitting

Extract those datasets from mice imputed as a list of data frames, and then pass them to the model fitting. The returned fitted model (from brm_multiple) is an ordinary brmsfit object. Therefore, the post-processing methods are straightforward without having to worry about pooling at all.

data("nhanes", package="mice")
head(nhanes)

library(mice)
library(brms)
imp <- mice(nhanes, m=5, print=FALSE)
fit_imp1 <- brm_multiple(bmi ~ age * chl, data=imp, chains=2)

summary(fit_imp1)
plot(fit_imp1, pars="^b")
## conditional_effects(fit_imp1, "age:chl")

Imputation during model fitting

Which variables contain missing values and how they should be predicted, and which of these imputed variables should be used as predictors.

bform <- bf(bmi | mi() ~ age * mi(chl)) +
  bf(chl | mi() ~ age) + set_rescor(FALSE)
fit_imp2 <- brm(bform, data = nhanes)

Multivariate Models

brms Multivariate Model

The term (1|p|fosternest) indicates a varying intercept over fosternest. Compute and store the LOO information criterion of the model for latter use of model comparisons.

fit1 <- brm(
  mvbind(tarsus, back) ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam),
  data = BTdata, chains = 2, cores = 2
)

fit1 <- add_criterion(fit1, "loo")
summary(fit1)

pp_check(fit1, resp = "tarsus") # posterior-predictive checks
pp_check(fit1, resp = "back")
bayes_R2(fit1) # Bayesian generalization of R Squared