Demo set-up

For demonstration purposes, let’s start by writing a minimal Stan model file lm.stan:

data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
y ~ normal(alpha + beta * x, sigma);
}

This Stan file encodes a simple linear regression model of the form:

\[ y_i \ \overset{\text{iid}}{\sim} \ N(\alpha + \beta \cdot x_i, \sigma^2), \quad i = 1,\ldots,N \]

Next, we create a small shiny-app contained in a single file app.R (in the same directory as lm.stan) that draws posterior samples from the Stan model in lm.stan with calls to rstan::stan_model() and rstan::sampling():

library(shiny)
library(rstan)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
numericInput("N", label = "N", value = 10)
),
mainPanel(
plotOutput("posteriors")
)
)
)
server <- function(input, output, session) {
## compile stan model
model <- stan_model(file = "lm.stan")
## draw samples
draws <- reactive({
N <- input$N
sampling(
object = model,
data = list(N = N, x = seq_len(N), y = rnorm(N, seq_len(N), 0.1)),
chains = 2,
iter = 1000
)
})
## plot histograms
output$posteriors <- renderPlot({
req(draws())
op <- par(mfrow = c(1, 2), cex = 1.25)
hist(extract(draws(), "alpha")[[1]], main = bquote("Posterior samples"~alpha), xlab = expression(alpha))
hist(extract(draws(), "beta")[[1]], main = bquote("Posterior samples"~beta), xlab = expression(beta))
par(op)
})
}
shinyApp(ui = ui, server = server)

The contents of this shiny-app can be summarized in a simple reactive graph:

New noisy responses \(y_1, \ldots, y_N\) are generated according to \(y_i \overset{\text{iid}}{\sim} N(\alpha + \beta \cdot x_i, \sigma^2)\) with \(x_i = i\) for each \(i = 1,\ldots, N\) and fixed parameters \(\alpha = 0\), \(\beta = 1\) and \(\sigma = 0.1\).

Slow model compilation

A fixed number of 2000 posterior samples for \(\alpha\) and \(\beta\) (and \(\sigma\)) is drawn across two individual chains (i.e. 1000 draws per chain), which should not take more than a few seconds to complete on an ordinary laptop computer, especially if \(N\) is small. However, when launching the shiny-app, it becomes evident that it takes significantly longer to complete drawing any initial posterior samples.

The reason for this is (obviously) that the Stan model has to be recompiled from the lm.stan file whenever we launch the shiny-app in a new R-session due to the call to rstan::stan_model(). Depending on the compiler settings, it takes up to ~1 minute to compile this single Stan model on my laptop computer, which more or less defeats the purpose of creating a shiny-app for interactive use.

Luckily, it is quite simple to avoid this unnecessary computational effort: we just have to compile our Stan models beforehand so that we can sample directly from the already compiled Stan models and skip the compilation step in rstan::stan_model(). Before describing a general R-package approach using rstantools, we start with a simpler approach –suitable for a large set of standard regression models– which is to take advantage of the pre-compiled Stan models in rstanarm.

Models created with brms

Besides rstanarm, the brms-package also provides a flexible interface to build Stan models directly using R syntax. The difference between rstanarm and brms, however, is that brms does not rely on pre-compiled Stan models and compiles generated .stan files on-the-fly. This provides additional flexibility with respect to rstanarm, but also means that calling brms::brm() directly in an interactive shiny-app suffers from the same unresponsiveness as rstan::stan_model().

As a workaround, we can call brms::make_stancode() to return the Stan program generated by brms:

brms::make_stancode(
formula = y ~ x,
data = data.frame(x = numeric(1), y = numeric(1)),
family = "gaussian"
)
#> // generated with brms 2.14.4
#> functions {
#> }
#> data {
#> int<lower=1> N; // total number of observations
#> vector[N] Y; // response variable
#> int<lower=1> K; // number of population-level effects
#> matrix[N, K] X; // population-level design matrix
#> int prior_only; // should the likelihood be ignored?
#> }
#> transformed data {
#> int Kc = K - 1;
#> matrix[N, Kc] Xc; // centered version of X without an intercept
#> vector[Kc] means_X; // column means of X before centering
#> for (i in 2:K) {
#> means_X[i - 1] = mean(X[, i]);
#> Xc[, i - 1] = X[, i] - means_X[i - 1];
#> }
#> }
#> parameters {
#> vector[Kc] b; // population-level effects
#> real Intercept; // temporary intercept for centered predictors
#> real<lower=0> sigma; // residual SD
#> }
#> transformed parameters {
#> }
#> model {
#> // likelihood including all constants
#> if (!prior_only) {
#> target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
#> }
#> // priors including all constants
#> target += student_t_lpdf(Intercept | 3, 0, 2.5);
#> target += student_t_lpdf(sigma | 3, 0, 2.5)
#> - 1 * student_t_lccdf(0 | 3, 0, 2.5);
#> }
#> generated quantities {
#> // actual population-level intercept
#> real b_Intercept = Intercept - dot_product(means_X, b);
#> }

By including this Stan code in the inst/stan folder and rebuilding the R-package, we circumvent the compilation step in brms::brm() and can directly sample from the compiled Stan model with rstan::sampling() as in the previous section. Note that the model input is slightly different, since brms has generated the Stan code for a more general multiple linear model:

system.time({
brms_fit <- rstan::sampling(
object = shinyStanModels:::stanmodels[["brms_lm"]],
data = list(N = 10L, ## number of observations
Y = rnorm(10, seq_len(10), 0.1), ## response vector
K = 2L, ## number of predictors
X = cbind(alpha = rep(1, 10), beta = seq_len(10)), ## predictor matrix
prior_only = FALSE ## set to TRUE to evaluate only the priors
),
chains = 2,
iter = 1000
)
})
#>
#> SAMPLING FOR MODEL 'brms_lm' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 6e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.010256 seconds (Warm-up)
#> Chain 1: 0.007216 seconds (Sampling)
#> Chain 1: 0.017472 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'brms_lm' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 1e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.010079 seconds (Warm-up)
#> Chain 2: 0.00762 seconds (Sampling)
#> Chain 2: 0.017699 seconds (Total)
#> Chain 2:
#> user system elapsed
#> 0.256 0.000 0.256
## alpha is b_Intercept
## beta is b[1]
rstan::summary(brms_fit, pars = c("b_Intercept", "b", "sigma"))[["summary"]]
#> mean se_mean sd 2.5% 25% 50%
#> b_Intercept 0.1362800 0.0025833278 0.07209531 -0.01553664 0.09773925 0.13826663
#> b[1] 0.9790349 0.0004261517 0.01127780 0.95785082 0.97213418 0.97865059
#> sigma 0.1014935 0.0019946898 0.03384615 0.05923768 0.07988739 0.09474943
#> 75% 97.5% n_eff Rhat
#> b_Intercept 0.1801695 0.2614545 778.8522 1.000521
#> b[1] 0.9853540 1.0031303 700.3586 1.001592
#> sigma 0.1133669 0.1869625 287.9173 1.005812

Remark: the Stan code generated by brms contains a bit of unnecessary complexity to sample from a simple linear model, but it does not take any effort to generate this Stan code, as we only need to provide the correct brms model syntax.

A note on deployment

When deploying the shiny-app to e.g. shinyapps.io or RStudio Connect using the rsconnect-package, the R-package generated with rstantools can be made available in a git repository on e.g. github or some other public repository, from which rsconnect is able to fetch and install the R-package when deploying the shiny-app.

Another solution to deploy and host shiny-apps on a server is ShinyProxy, which launches shiny-apps from individual Docker containers. By installing the R-package generated by rstantools when building the Docker image of the shiny-app, we ensure that we can directly sample from our compiled Stan models whenever a new Docker container is started. The following Dockerfile provides a minimal template to install the rstantools-generated R-package from a bundled (.tar.gz) package file and serve the shiny-app at port 3838:

FROM rocker/r-ver:4.0.3
# install system dependencies
RUN apt-get update && \
apt-get install -y --no-install-recommends \
libv8-dev && \
apt-get clean && \
rm -rf /var/lib/apt/lists/
# install R packages (using littler)
# this assumes .tar.gz exists in same folder as Dockerfile
COPY shinyStanModels_0.1.tar.gz ./
RUN install2.r --error shiny rstan rstantools && \
install2.r --error shinyStanModels_0.1.tar.gz && \
rm *.tar.gz
EXPOSE 3838
CMD ["R", "-e", "shiny::runApp(appDir = system.file('app', package = 'shinyStanModels'), port = 3838, host = '0.0.0.0')"]