Running Fits Only When Something Changes

You’ll see file='...' in the longer running of the blrm calls below. If the file already exists and none of the data nor the options sent to rstan nor the underlying Stan code have changed from what were used to create the fit object stored in that file (as judged by their md5 hash), that saved fit object is returned immediately without running the rstan code, often saving a great deal of execution time. This works well with the workflow of long R markdown reports making it so that only portions of Bayesian analysis that have changed are run. Note that using the knitr cache=TRUE option does not work well as the cache files for this script were about a GB in size, and the system does not accurately recognize when a model fit hasn’t changed and doesn’t need to be run when rstan is involved.

A restored fit object does not contain the rstan object, saving tens of megabytes of storage. Standard Stan diagnostics are stored in the fit object separately, and it is assumed that if the user wanted to run rstan::loo loo=TRUE would have been specified to blrm so that loo is run and its (small) result stored in the blrm fit object. If you want to run pairs to get more graphical diagnostics, intead of relying on the rstan object always being available, specify pairs=TRUE or pairs='filename.png' to blrm to graph the pair plots. The latter is recommended, and one can put knitr::include_graphics('filename.png') in the R code chunk to render the graph in the report even if blrm was not re-run.

When file='filename.rds' is specified to blrm and the file does not exist or analysis inputs have changed since the file was created, the blrm fit object will be saved in saveRDS .rds binary format in your current working directory at the completion of blrm’s work. The rstan component is omitted from the saved file.

The utility function fitIf available from here is another way to efficiently manage Bayesian simulation workflow. The fitIf approach runs the analysis only if the file doesn’t already exist. If it exists but the data or the model have changed, the fitIf approach is not intelligent enough to re-run the analysis (unlike the file='...' approach above). Whether using fitIf or file=, you lose the ability to run the rms::stanDxplot(..., rstan=TRUE) on restored objects, so the stanDxplot function tries to find an existing trace plot image file corresponding to the current R markdown chunk when the fit object no longer has the rstan component. For most purposes this doesn’t matter, because running using the defaults shows the needed non-burnin samples which are always stored in fit objects. The file=filename.rds approach is preferred.

Example: 10-level Ordinal Outcome

Simulate a dataset with three predictors and one 10-level ordinal outcome. Run the frequentist proportional odds model then the Bayesian one.

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

Logistic Regression Model

 lrm(formula = y ~ x1 + pol(x2, 2) + x3, eps = 1e-07)
 

Frequencies of Responses

  1  2  3  4  5  6  7  8  9 10 
 50 50 50 50 50 50 50 50 50 50 
 

Before getting posterior distributions of parameters, use rstan to just get maximum likelihood estimates and compare them with those from lrm. Do this for increasingly flat priors for the \(\beta\)s. The default prior SD for blrm is 100. Running method='optimizing' is a quick way to study the effect of priors on the posterior modes for non-intercepts when there are no random effects in the model. It is also a way to get a sense of the scaling of parameters used in the actual calculations (but see a later section for how to show standard deviations of variables on the transformed scale given to Stan). Except for variables listed in the keepsep argument to blrm, prior standard deviations apply to the QR orthonormalized version of the design matrix.

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)
}

Prior SD: 0.25 
-2 log likelihood: 2269.378 
        y>=2         y>=3         y>=4         y>=5         y>=6         y>=7 
 2.191046301  1.378643586  0.838324839  0.395207927 -0.011602895 -0.418570617 
        y>=8         y>=9        y>=10           x1           x2         x2^2 
-0.862048892 -1.402835097 -2.215654082  0.024804025  0.013114821 -0.003458598 
          x3 
 0.026136325 

Prior SD: 1 
-2 log likelihood: 2234.618 
       y>=2        y>=3        y>=4        y>=5        y>=6        y>=7 
 2.12797711  1.29550040  0.73699395  0.27627257 -0.14878702 -0.57642046 
       y>=8        y>=9       y>=10          x1          x2        x2^2 
-1.04203046 -1.60705055 -2.44522185  0.31969267  0.16828142 -0.04299895 
         x3 
 0.33452830 

Prior SD: 10 
-2 log likelihood: 2076.618 
      y>=2       y>=3       y>=4       y>=5       y>=6       y>=7       y>=8 
 2.0518358  1.0895535  0.4153739 -0.1563173 -0.6939620 -1.2471658 -1.8490486 
      y>=9      y>=10         x1         x2       x2^2         x3 
-2.5683076 -3.5681608  1.5016729  0.7577512 -0.2091988  1.5632871 

Prior SD: 100 
-2 log likelihood: 2086.232 
      y>=2       y>=3       y>=4       y>=5       y>=6       y>=7       y>=8 
 2.0546496  1.0827924  0.4002833 -0.1792460 -0.7247696 -1.2867285 -1.8981788 
      y>=9      y>=10         x1         x2       x2^2         x3 
-2.6286338 -3.6405325  1.5685655  0.7895744 -0.2202702  1.6335934 

Prior SD: 10000 
-2 log likelihood: 2122.98 
      y>=2       y>=3       y>=4       y>=5       y>=6       y>=7       y>=8 
 2.0551518  1.0831610  0.4005549 -0.1790191 -0.7246383 -1.2866789 -1.8982289 
      y>=9      y>=10         x1         x2       x2^2         x3 
-2.6288108 -3.6408433  1.5695440  0.7896154 -0.2213424  1.6341990 

# Compare with ordinary MLEs and deviance
f$deviance

[1] 2302.585 2067.549

coef(f)

      y>=2       y>=3       y>=4       y>=5       y>=6       y>=7       y>=8 
 2.0551529  1.0831900  0.4005868 -0.1790073 -0.7246249 -1.2866643 -1.8982180 
      y>=9      y>=10         x1         x2       x2^2         x3 
-2.6287996 -3.6408392  1.5694325  0.7900802 -0.2212900  1.6340800 

Fit the model with Dirichlet priors on intercepts and wide normal priors on the \(\beta\)s. Show the model fit summary. Note that the indexes of predictive discrimination/accuracy include 0.95 highest posterior density intervals. In frequentist inference we pretend that quantities such as AUROC and \(R^2\) are estimated without error, which is far from the case.

In several places you will see an index named Symmetry. This is a measure of the symmetry of a posterior distribution. Values farther from 1.0 indicate asymmetry, which indicates that the use of standard errors and the use of a normal approximation for the posterior distribution are not justified. The symmetry index is the ratio of the gap between the posterior mean and the 0.95 quantile of the posterior distribution to the gap between the 0.05 quantile and the mean.

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

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.233 for Intercepts

 blrm(formula = y ~ x1 + pol(x2, 2) + x3, file = "bs.rds")
 

Frequencies of Responses

  1  2  3  4  5  6  7  8  9 10 
 50 50 50 50 50 50 50 50 50 50 
 

# Show more detailed analysis of model performance measures
blrmStats(bs, pl=TRUE)

Indexes computed for a random sample of 400 of 4000 posterior draws

gp, B, EV, and vp are for intercept 1 out of 9 intercepts

           Dxy     C     g    gp     B    EV     v    vp
Mean     0.477 0.738 1.571 0.114 0.077 0.147 1.897 0.014
SE       0.003 0.002 0.112 0.014 0.001 0.023 0.271 0.003
Lower    0.471 0.735 1.383 0.088 0.076 0.110 1.416 0.008
Upper    0.481 0.740 1.806 0.141 0.079 0.196 2.432 0.020
Symmetry 0.423 0.423 1.097 1.167 1.482 1.226 1.186 1.314

Dxy: 2*(C - 0.5)   C: concordance probability
g: Gini mean |difference| on linear predictor (lp)
gp: Gini on predicted probability        B: Brier score
EV: explained variation on prob. scale   v: var(lp)   vp: var(prob)

Show basic Stan diagnostics. Had stanDxplots(bs, rstan=TRUE) been used, intercepts would have been shifted from what is in g because of subtractions of covariate means before passing data to rstan.

stanDxplot(bs)
stanDx(bs)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)

      n_eff  Rhat
y>=2   7163 1.000
y>=3   5962 1.000
y>=4   7114 1.000
y>=5   8307 0.999
y>=6  10963 0.999
y>=7   9538 0.999
y>=8   8619 0.999
y>=9   8206 0.999
y>=10  7256 0.999
x1     6646 0.999
x2     9568 0.999
x2^2   8848 1.000
x3     5968 1.000

Here are the posterior distributions, calculated using kernel density estimates from posterior draws. Posterior models, shown as vertical lines, are parameter values that maximize the log posterior density (using rstan::optimizing in the original model fit) so do not necessarily coincide with the peak of the kernel density estimates.

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))

             MLE       mode       mean     median
y>=2   2.0551529  2.0546496  2.0608013  2.0583224
y>=3   1.0831900  1.0827924  1.0835423  1.0835161
y>=4   0.4005868  0.4002833  0.3980935  0.4010254
y>=5  -0.1790073 -0.1792460 -0.1825232 -0.1832366
y>=6  -0.7246249 -0.7247696 -0.7302207 -0.7303855
y>=7  -1.2866643 -1.2867285 -1.2945701 -1.2941131
y>=8  -1.8982180 -1.8981788 -1.9063342 -1.9034973
y>=9  -2.6287996 -2.6286338 -2.6387192 -2.6354563
y>=10 -3.6408392 -3.6405325 -3.6543224 -3.6510901
x1     1.5694325  1.5685655  1.5746271  1.5724682
x2     0.7900802  0.7895744  0.7916821  0.7942466
x2^2  -0.2212900 -0.2202702 -0.2220469 -0.2225900
x3     1.6340800  1.6335934  1.6404096  1.6372988

# Compare covariance matrix of posterior draws with MLE
round(diag(vcov(f)) / diag(vcov(bs)), 2)

 y>=2  y>=3  y>=4  y>=5  y>=6  y>=7  y>=8  y>=9 y>=10    x1    x2  x2^2    x3 
 1.01  1.04  1.04  0.98  0.99  1.01  1.04  1.00  0.96  0.96  0.98  1.06  0.95 

range(vcov(f) / vcov(bs))

[1] -0.8151463  1.5075811