Accessing the contents of a stanfit object

Stan Development Team

2023-09-08

Source: vignettes/stanfit_objects.Rmd

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 1 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]  8.158595  7.880213 14.086734  1.992156  7.283674 15.916944
head(list_of_draws$tau)
[1] 11.5467890  3.1630415  8.9115022  0.7782174  3.3101477  4.7479295
head(list_of_draws$theta)
          
iterations       [,1]      [,2]      [,3]      [,4]       [,5]      [,6]
      [1,] 28.4498512 13.812518 27.569254  9.143677  15.323212 18.270046
      [2,]  9.4066896  9.594654  5.298375  5.605437   7.827925  6.899140
      [3,] 14.0018236 18.129483 27.020226 15.679345 -10.604681 -3.521928
      [4,]  0.8767176  2.344286  1.289285  1.479599   1.725346  2.312681
      [5,]  4.9130196  5.703014  9.707277  9.056205   4.472401  9.378470
      [6,] 15.0274075  5.513797 16.454255 14.355015  21.074407  9.402672
          
iterations      [,7]      [,8]
      [1,]  1.479218 15.372751
      [2,] 14.243631  1.350251
      [3,] 25.729921 12.187797
      [4,]  2.758654  2.501818
      [5,]  7.335756  9.839808
      [6,] 16.360368 11.608590


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,] 12.822571 15.297461
      [2,]  9.392738 11.477776
      [3,]  6.081793  5.126646
      [4,]  7.517083  7.173908
      [5,]  8.536018 24.396196
      [6,] 12.179750  9.792373


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.88112638 0.11886981 5.1400521  -2.1874849   4.5862578
tau        6.45264757 0.14793579 5.6059239   0.2438317   2.4334292
eta[1]     0.36243921 0.01534587 0.9281426  -1.4718403  -0.2628066
eta[2]     0.01886754 0.01341454 0.8668440  -1.6882721  -0.5607231
eta[3]    -0.16085954 0.01555826 0.9356425  -1.9333095  -0.8083447
eta[4]    -0.02161591 0.01528133 0.8929354  -1.7499464  -0.6141818
eta[5]    -0.37753813 0.01462871 0.8545046  -2.0589638  -0.9574080
eta[6]    -0.21296151 0.01452506 0.9063613  -1.9594580  -0.8191162
eta[7]     0.34961133 0.01448700 0.9062853  -1.5323017  -0.2308520
eta[8]     0.06684927 0.01482872 0.9363868  -1.8603479  -0.5547307
theta[1]  11.14131968 0.15570247 8.2709442  -2.6106520   6.0436778
theta[2]   7.95859334 0.09304677 6.2312790  -4.3871431   3.9185803
theta[3]   6.41387106 0.12989815 7.6799788 -11.0819752   2.2992616
theta[4]   7.60884404 0.10202428 6.5850034  -6.0371411   3.5611151
theta[5]   4.94911661 0.10464727 6.2732871  -9.3887854   1.2949771
theta[6]   6.20658841 0.10725834 6.8355631  -8.3747947   2.1210493
theta[7]  10.60462091 0.12250274 6.6961993  -0.6030978   6.1064082
theta[8]   8.47344823 0.14454695 8.0284266  -7.4842285   3.8346619
lp__     -39.58643914 0.06863020 2.6044560 -45.5215532 -41.1448851
                   50%         75%      97.5%    n_eff      Rhat
mu         7.787782492  11.1431031  18.388916 1869.786 0.9995491
tau        5.088594901   8.9653879  20.987474 1435.978 1.0029909
eta[1]     0.373245147   0.9923146   2.129075 3658.022 1.0004282
eta[2]    -0.004334584   0.5905507   1.717842 4175.711 0.9997220
eta[3]    -0.181754770   0.4567332   1.691041 3616.576 1.0000403
eta[4]    -0.024249343   0.5771657   1.765282 3414.425 0.9997084
eta[5]    -0.390550370   0.1578636   1.371040 3412.063 0.9995389
eta[6]    -0.217677125   0.3668237   1.612039 3893.739 0.9999304
eta[7]     0.357785797   0.9556616   2.131229 3913.571 0.9998760
eta[8]     0.095595700   0.7077347   1.812303 3987.523 1.0010843
theta[1]   9.983158055  15.3383770  30.678479 2821.754 1.0005669
theta[2]   7.915022564  11.7913960  20.958119 4484.890 0.9995185
theta[3]   6.792102722  10.9150764  21.121257 3495.539 1.0018001
theta[4]   7.645471822  11.8346127  20.775881 4165.862 0.9996712
theta[5]   5.393538462   9.1365562  16.129178 3593.640 0.9997582
theta[6]   6.585073185  10.5238483  19.289236 4061.500 0.9999494
theta[7]   9.962639496  14.5929002  25.730154 2987.894 1.0005973
theta[8]   8.139927461  12.8687269  25.952480 3084.914 1.0010763
lp__     -39.331531763 -37.7405538 -35.308693 1440.136 1.0009663

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.881126 0.1188698 5.140052 1.5573901 14.50957 1869.786 0.9995491
tau 6.452648 0.1479358 5.605924 0.9572862 13.24043 1435.978 1.0029909

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.5573901 14.50957
tau 0.9572862 13.24043


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.8514222 0.8829814 0.8919147 0.8416380
max_treedepth_by_chain <- sapply(sampler_params, function(x) max(x[, "treedepth__"]))
print(max_treedepth_by_chain)
[1] 4 5 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] 1.830195

$tau
[1] 1.48977

$eta
[1] -1.9644185  1.2900084 -0.5047292  1.3198299 -1.9966361 -1.9750366  1.8860891
[8]  1.8957488

$theta
[1] -1.096337  3.752011  1.078264  3.796438 -1.144334 -1.112156  4.640034
[8]  4.654425


(P)RNG seed

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

[1] 1466491891


Warmup and sampling times

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

        warmup sample
chain:1  0.050  0.042
chain:2  0.047  0.051
chain:3  0.051  0.049
chain:4  0.046  0.043