Accessing the contents of a stanfit object

This vignette demonstrates how to access most of data stored in a stanfit object. A stanfit object (an object of class "stanfit") contains the output derived from fitting a Stan model using Markov chain Monte Carlo or one of Stan’s variational approximations (meanfield or full-rank). Throughout the document we’ll use the stanfit object obtained from fitting the Eight Schools example model:

library(rstan)
fit <- stan_demo("eight_schools", refresh = 0)

Warning: There were 5 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

class(fit)

[1] "stanfit"
attr(,"package")
[1] "rstan"

Posterior draws

There are several functions that can be used to access the draws from the posterior distribution stored in a stanfit object. These are extract, as.matrix, as.data.frame, and as.array, each of which returns the draws in a different format.

extract()

The extract function (with its default arguments) returns a list with named components corresponding to the model parameters.

list_of_draws <- extract(fit)
print(names(list_of_draws))

[1] "mu"    "tau"   "eta"   "theta" "lp__"

In this model the parameters mu and tau are scalars and theta is a vector with eight elements. This means that the draws for mu and tau will be vectors (with length equal to the number of post-warmup iterations times the number of chains) and the draws for theta will be a matrix, with each column corresponding to one of the eight components:

head(list_of_draws$mu)

[1]  7.429953 14.004288  6.808272  7.770409 14.035968  6.316280

head(list_of_draws$tau)

[1] 0.5625588 5.6883492 8.2123765 3.7490770 4.7368366 2.7055457

head(list_of_draws$theta)

          
iterations      [,1]      [,2]      [,3]       [,4]      [,5]      [,6]
      [1,]  8.051981  8.067713  6.974470  7.4043922  7.203581  6.929355
      [2,] 22.194709  8.919106  6.271111 18.6575729 12.443539 10.287069
      [3,] 16.968542 13.076664  6.326293  0.5240619  4.054372  4.065661
      [4,]  4.609341  2.651477 13.655380  7.5270801  5.192615 13.136100
      [5,] 21.278883  8.022200 13.907137 11.9597867 14.160836 16.948794
      [6,]  5.701759  8.459602  9.437653  7.9950718  6.874026  4.403695
          
iterations      [,7]       [,8]
      [1,]  8.126263  8.0232804
      [2,] 14.817760 18.1436859
      [3,] 19.097667  3.3680835
      [4,] 16.347559 -0.2611905
      [5,] 16.379834 18.8709982
      [6,]  8.003112  5.8278579

as.matrix(), as.data.frame(), as.array()

The as.matrix, as.data.frame, and as.array functions can also be used to retrieve the posterior draws from a stanfit object:

matrix_of_draws <- as.matrix(fit)
print(colnames(matrix_of_draws))

 [1] "mu"       "tau"      "eta[1]"   "eta[2]"   "eta[3]"   "eta[4]"  
 [7] "eta[5]"   "eta[6]"   "eta[7]"   "eta[8]"   "theta[1]" "theta[2]"
[13] "theta[3]" "theta[4]" "theta[5]" "theta[6]" "theta[7]" "theta[8]"
[19] "lp__"

df_of_draws <- as.data.frame(fit)
print(colnames(df_of_draws))

 [1] "mu"       "tau"      "eta[1]"   "eta[2]"   "eta[3]"   "eta[4]"  
 [7] "eta[5]"   "eta[6]"   "eta[7]"   "eta[8]"   "theta[1]" "theta[2]"
[13] "theta[3]" "theta[4]" "theta[5]" "theta[6]" "theta[7]" "theta[8]"
[19] "lp__"

array_of_draws <- as.array(fit)
print(dimnames(array_of_draws))

$iterations
NULL

$chains
[1] "chain:1" "chain:2" "chain:3" "chain:4"

$parameters
 [1] "mu"       "tau"      "eta[1]"   "eta[2]"   "eta[3]"   "eta[4]"  
 [7] "eta[5]"   "eta[6]"   "eta[7]"   "eta[8]"   "theta[1]" "theta[2]"
[13] "theta[3]" "theta[4]" "theta[5]" "theta[6]" "theta[7]" "theta[8]"
[19] "lp__"

The as.matrix and as.data.frame methods essentially return the same thing except in matrix and data frame form, respectively. The as.array method returns the draws from each chain separately and so has an additional dimension:

print(dim(matrix_of_draws))
print(dim(df_of_draws))
print(dim(array_of_draws))

[1] 4000   19
[1] 4000   19
[1] 1000    4   19

By default all of the functions for retrieving the posterior draws return the draws for all parameters (and generated quantities). The optional argument pars (a character vector) can be used if only a subset of the parameters is desired, for example:

mu_and_theta1 <- as.matrix(fit, pars = c("mu", "theta[1]"))
head(mu_and_theta1)

          parameters
iterations         mu   theta[1]
      [1,] 13.1290441 15.2248631
      [2,]  3.5303972  0.5025017
      [3,] 13.0267306 13.3777897
      [4,]  2.4144074  2.2301669
      [5,] 16.3292762 16.7689955
      [6,]  0.1242002 -5.9148130

Posterior summary statistics and convergence diagnostics

Summary statistics are obtained using the summary function. The object returned is a list with two components:

fit_summary <- summary(fit)
print(names(fit_summary))

[1] "summary"   "c_summary"

In fit_summary$summary all chains are merged whereas fit_summary$c_summary contains summaries for each chain individually. Typically we want the summary for all chains merged, which is what we’ll focus on here.

The summary is a matrix with rows corresponding to parameters and columns to the various summary quantities. These include the posterior mean, the posterior standard deviation, and various quantiles computed from the draws. The probs argument can be used to specify which quantiles to compute and pars can be used to specify a subset of parameters to include in the summary.

For models fit using MCMC, also included in the summary are the Monte Carlo standard error (se_mean), the effective sample size (n_eff), and the R-hat statistic (Rhat).

print(fit_summary$summary)

                  mean    se_mean        sd        2.5%         25%
mu         7.845399067 0.09862177 4.9862077  -1.7798694   4.5126697
tau        6.439273843 0.12174256 5.3068458   0.2149051   2.4354665
eta[1]     0.404820903 0.01430157 0.9345049  -1.4731718  -0.1996558
eta[2]    -0.002291299 0.01374855 0.8726460  -1.7932811  -0.5893252
eta[3]    -0.178468622 0.01384425 0.9505377  -2.0458986  -0.8151339
eta[4]    -0.030199852 0.01308819 0.9024600  -1.7684224  -0.6168620
eta[5]    -0.326338960 0.01418823 0.8898505  -2.0885471  -0.8832006
eta[6]    -0.211740930 0.01444361 0.8977131  -1.9907735  -0.7989340
eta[7]     0.347082918 0.01387133 0.8893795  -1.4910397  -0.2069065
eta[8]     0.065895475 0.01409353 0.9097974  -1.6909845  -0.5464345
theta[1]  11.319899741 0.14200463 8.1831643  -2.2511892   5.9402038
theta[2]   7.904337692 0.09193808 6.3857133  -4.5651925   3.8349507
theta[3]   6.184096063 0.11762367 7.8765378 -11.4252763   1.9758095
theta[4]   7.574931220 0.09515011 6.6082253  -5.7350551   3.6562614
theta[5]   5.125211867 0.09799325 6.3481408  -9.1552087   1.4180764
theta[6]   6.159564566 0.10122425 6.7833773  -8.8154700   2.3718718
theta[7]  10.510641417 0.10563894 6.7783407  -1.6235518   6.0548067
theta[8]   8.370052438 0.12676547 7.7310331  -7.1068530   3.9089512
lp__     -39.614874594 0.08104245 2.6988489 -45.5156888 -41.2336905
                  50%         75%      97.5%    n_eff      Rhat
mu         7.90179850  11.1001664  17.443809 2556.202 0.9997484
tau        5.29806075   9.0885628  19.878020 1900.151 1.0023480
eta[1]     0.42883352   1.0336297   2.243369 4269.685 0.9994841
eta[2]     0.01343915   0.5815423   1.658351 4028.677 0.9995530
eta[3]    -0.18078654   0.4486399   1.731511 4714.112 0.9997082
eta[4]    -0.04269481   0.5653958   1.842109 4754.410 0.9993020
eta[5]    -0.35307052   0.2085373   1.566812 3933.486 0.9999215
eta[6]    -0.21243408   0.3630094   1.549097 3862.989 1.0000187
eta[7]     0.36904597   0.9251363   2.095560 4110.910 1.0001905
eta[8]     0.07064734   0.6689094   1.789397 4167.255 1.0006128
theta[1]  10.28882575  15.5428855  31.054846 3320.760 0.9997089
theta[2]   7.88548987  11.9397159  20.743545 4824.229 0.9999034
theta[3]   6.72258700  11.1234670  20.425373 4484.161 0.9993256
theta[4]   7.74122164  11.7051569  20.409569 4823.376 0.9995895
theta[5]   5.52357562   9.4161936  16.164821 4196.631 0.9997148
theta[6]   6.60817792  10.4530636  18.568481 4490.791 1.0001590
theta[7]  10.03977790  14.5248710  25.478101 4117.169 0.9999521
theta[8]   8.11909890  12.6710419  25.031445 3719.402 1.0007129
lp__     -39.33491698 -37.6442584 -35.047399 1109.001 1.0052792

If, for example, we wanted the only quantiles included to be 10% and 90%, and for only the parameters included to be mu and tau, we would specify that like this:

mu_tau_summary <- summary(fit, pars = c("mu", "tau"), probs = c(0.1, 0.9))$summary
print(mu_tau_summary)

        mean    se_mean       sd       10%      90%    n_eff      Rhat
mu  7.845399 0.09862177 4.986208 1.4446008 14.04145 2556.202 0.9997484
tau 6.439274 0.12174256 5.306846 0.9349537 13.35709 1900.151 1.0023480

Since mu_tau_summary is a matrix we can pull out columns using their names:

mu_tau_80pct <- mu_tau_summary[, c("10%", "90%")]
print(mu_tau_80pct)

          10%      90%
mu  1.4446008 14.04145
tau 0.9349537 13.35709

Sampler diagnostics

For models fit using MCMC the stanfit object will also contain the values of parameters used for the sampler. The get_sampler_params function can be used to access this information.

The object returned by get_sampler_params is a list with one component (a matrix) per chain. Each of the matrices has number of columns corresponding to the number of sampler parameters and the column names provide the parameter names. The optional argument inc_warmup (defaulting to TRUE) indicates whether to include the warmup period.

sampler_params <- get_sampler_params(fit, inc_warmup = FALSE)
sampler_params_chain1 <- sampler_params[[1]]
colnames(sampler_params_chain1)

[1] "accept_stat__" "stepsize__"    "treedepth__"   "n_leapfrog__" 
[5] "divergent__"   "energy__"

To do things like calculate the average value of accept_stat__ for each chain (or the maximum value of treedepth__ for each chain if using the NUTS algorithm, etc.) the sapply function is useful as it will apply the same function to each component of sampler_params:

mean_accept_stat_by_chain <- sapply(sampler_params, function(x) mean(x[, "accept_stat__"]))
print(mean_accept_stat_by_chain)

[1] 0.8384776 0.7451438 0.8929761 0.9130661

max_treedepth_by_chain <- sapply(sampler_params, function(x) max(x[, "treedepth__"]))
print(max_treedepth_by_chain)

[1] 4 4 4 4

Model code

The Stan program itself is also stored in the stanfit object and can be accessed using get_stancode:

code <- get_stancode(fit)

The object code is a single string and is not very intelligible when printed:

print(code)

[1] "data {\n  int<lower=0> J;          // number of schools\n  real y[J];               // estimated treatment effects\n  real<lower=0> sigma[J];  // s.e. of effect estimates\n}\nparameters {\n  real mu;\n  real<lower=0> tau;\n  vector[J] eta;\n}\ntransformed parameters {\n  vector[J] theta;\n  theta = mu + tau * eta;\n}\nmodel {\n  target += normal_lpdf(eta | 0, 1);\n  target += normal_lpdf(y | theta, sigma);\n}"
attr(,"model_name2")
[1] "schools"

A readable version can be printed using cat:

cat(code)

data {
  int<lower=0> J;          // number of schools
  real y[J];               // estimated treatment effects
  real<lower=0> sigma[J];  // s.e. of effect estimates
}
parameters {
  real mu;
  real<lower=0> tau;
  vector[J] eta;
}
transformed parameters {
  vector[J] theta;
  theta = mu + tau * eta;
}
model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

Initial values

The get_inits function returns initial values as a list with one component per chain. Each component is itself a (named) list containing the initial values for each parameter for the corresponding chain:

inits <- get_inits(fit)
inits_chain1 <- inits[[1]]
print(inits_chain1)

$mu
[1] 0.7149569

$tau
[1] 4.408078

$eta
[1] -1.1341982  0.5930025 -0.3489150  1.8952199  0.4201683  1.5565210  1.9530570
[8]  1.6872806

$theta
[1] -4.2846770  3.3289582 -0.8230876  9.0692342  2.5670914  7.5762231  9.3241848
[8]  8.1526217

(P)RNG seed

The get_seed function returns the (P)RNG seed as an integer:

print(get_seed(fit))

[1] 1824311018

Warmup and sampling times

The get_elapsed_time function returns a matrix with the warmup and sampling times for each chain:

print(get_elapsed_time(fit))

        warmup sample
chain:1  0.019  0.019
chain:2  0.016  0.011
chain:3  0.017  0.019
chain:4  0.016  0.019