If all the nuts you ever have to tighten are the same one or two sizes then you buy the best spanners in those sizes. You don’t use a monkey wrench that needs to be set each time. But if you’re encountering different sizes of nuts, fittings, and other objects that need tightening then you get a monkey wrench. In statistical modelling, Stan is the monkey wrench.

According to https://mc-Stan.org/: “Stan is a state-of-the-art platform for statistical modelling and high-performance statistical computation.”

You type out your model in Stan language and then Stan does all the work, including: simulating the whole posterior distribution, simulating the approximate posterior, and/or Maximum Likelihood.

• It compiles your model into C++ code, so it’s very fast
• Latest simulation methods
• Can do usual Bayes, Variational Bayes or Maximum Likelihood as you like - You don’t change your model, just one line of code
• Most importantly, if you want to change your model then you change only the part of the model that needs changing. No changing software, no installing new packages, no changing notation for a different ‘PROC’, etc.

Stan is currently available for R, Python, and C++ as far as I can tell. We will focus on rstan as it’s the easiest for most statisticians.

There are many ways to use Stan in R:

• You can type your Stan model in a separate file, then tell rstan to compile that file
• You can type your Stan model in a Stan block in your Rmarkdown document and click the play button to compile it
• You can save the compiled model in a .Rds or Rdata file and then load the compiled model each time you need it
• You can mix and match the above as it suits you

I will use Rmarkdown to demonstrate.

I use Rmarkdown for almost everything, including teaching, assessment, consultation, some research, and even making presentations like this one.

To use Rmarkdown you open RStudio, click File -> New File -> R Markdown…

To put in a Stan model you click on the Insert button at the top right of the code window, then click Stan.

On the next slide I’ll show a very basic Stan model in RStudio.

Notice the play button in the top right corner (circled in red).

The critical components are:

1. The data
2. The parameters
3. The likelihood

The optional components are:

• Priors for some or all parameters - default is uniform on the domain of the parameter
• Transformed Parameters - for when your model gets complicated
• Transformed Data - better to do in R in advance
• Generated Quantities - for creating things you want later

To be able to use rstan the way I do:

1. Install R
2. Install Rtools
3. Install RStudio
4. In the RStudio Console window, type install.packages(‘rstan’)
5. In your code, type library(rstan)

• Optionally, set the number of cores for Stan to use when simulating in parallel
• Easiest is to pick a number: options(mc.cores = 3)
• On most computers you can get a big speed increase easily if you match the number of cores and chains
• So calculate and store a good number: library(parallel); mycores <- ceiling(0.75*detectCores(logical = FALSE))
• Then use that number everywhere: options(mc.cores = mycores) and later on when simulating set chains = mycores
• I like 75% because then I can still use my computer while simulating (i.e. watch YouTube while I’m technically working).

Let’s apply the t model to some data.

library(openxlsx)
sourcedata <- read.xlsx("IBM19891998.xlsx")
lret <- sourcedata$LogReturns
dates <- sourcedata$year + sourcedata$month/12 + sourcedata$day/365

library(viridisLite)
cols <- viridis(3)
par(mfrow = c(2, 1), mar = c(4, 4, 2, 0.2))
plot(dates, lret, type = "l", col = cols[1], lwd = 1, xlab = "Date", 
    ylab = "Log return", main = "IBM")
hist(lret, 50, col = cols[2], density = 20, ylab = "Frequency", 
    xlab = "Log return", main = "IBM")

stan_data <- list(x = lret, n = length(lret))  # Create the data list

fitML1 <- optimizing(studentt, stan_data)  # ML estimation

# This is the main simulation function:
fitBayes1 <- sampling(studentt, stan_data, chains = mycores)
# Type ?sampling for a lot of additional options

saveRDS(fitBayes1, "studenttIBM.Rds")  # Save model fit

summary(fitBayes1)$summary  # See results of simulation

traceplot(fitBayes1)  # See simulation paths to check for stability

list_of_draws1 <- extract(fitBayes1)  # Extract simulations

The Rhat measure should always be close to 1.

Rhat compares the variation between chains to the variation within chains. If these are the same then that is a good indication of convergence in distribution.

We don’t need to only use standard distributions in Stan.

For example, we could use a better prior on the degrees of freedom (provided by Prof Abrie van der Merwe):

\[p(\nu)\propto \left[\psi'(0.5\nu)-\psi'(0.5(\nu+1))-2(\nu+3)\nu^{-1}(\nu+1)^{-2}\right]^{0.5}\]

where \(\psi'\) is the trigamma function.

A standard regression model is just fitting a normal distribution with a linear function for the mean parameter. From Part 1.1 of the Stan user manual:

But with Stan we can fit any functions on any parameters.

                 (Within reason obviously, simpler is better.)

Opportunistic potatoes are a pest that affects farms in many negative ways and need to be controlled. This experiment attempts to find the economically optimal dose of a specific treatment for this purpose.

In the original experiment 8 measures were taken for various cultivars and doses, but we will restrict ourselves to only one measure and combine the cultivars for this discussion.

This research was done by Talana Cronje, UFS, on behalf of Potatoes South Africa.

They applied the treatment at various proportions of the recommended dose. Lower proportions are cheaper to implement.

You must enable Javascript to view this page properly.

• Looking at the data, it is clear that the economically optimal dose is near 20%.
• More than that costs more money but doesn’t yield much improvement.
• After 40% there doesn’t appear to be any real improvement.
• No need for a complex cost function analysis in this case.

Can we get the same answer more formally?

I propose the following model:

Given an unknown change point \(\gamma\), the model is formally defined as \(y_i\sim N(\beta_0+\beta_1(x_i-\gamma)I(x_i>\gamma)+\beta_2(x_i-\gamma)I(x_i<\gamma),\ \sigma^2)\).

• By far the easiest way to fit the model (in my opinion) is using STAN, via rstan
• The p-value for the change point test is: 0.0047 and for the test of a flat slope after the change point the p-value is: 0.6583
• Which means we assume there is a change point
  • with expected value 24% of recommended treatment
  • and 95% HPD interval of 10% to 42%

After fitting the model, we can obtain intervals for any quantity of interest as follows:

• Calculate the quantity for each simulated parameter vector, to get a vector of quantities
• This vector of quantities can be seen as an approximation to the posterior distribution of the quantity of interest, let us call it sims.
• Calculate a symmetric interval, say 95%, using the built-in quantile function: quantile(sims,c(0.025,0.975))
• Calculate an equivalent to the traditional two-sided p-value using this function:

pvalfunc <- function(sims, target = 0) {
    2 * min(mean(sims < target), mean(sims > target))
}

If we take logs then the zeros are a problem. We can consider zeros as censored below 1 or 0.5.

Here we want to do a regression on the scale parameter as it changes with time and claimant properties.

We were asked to help rank 200 types of wheat according to some measures of how well they resist different pathogens.

While we did various standard analyses, the measurements were quite skew and seemed like they should fit a gamma model. The measurements were taken under varying conditions, such as site, year, and pathogen. Ultimately though we need a single measure of resistance per wheat strain under general conditions.

One possible approach is to model the measurements jointly, with wheat strain as fixed effect and the nuisance factors as random effects. We also have to account for differing variances under differing conditions.