Accessing the contents of a stanfit object

Stan Development Team

2020-02-10

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 7 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
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.

[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:

[1] 14.143929  8.478591  6.504845 11.970911  8.524122  8.160192
[1]  4.3662195  0.4272034  4.1657276 12.2439305  4.3656120  4.5503272
          
iterations      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
      [1,] 18.368073 13.625011 13.460702 13.675479 15.789049 12.997669
      [2,]  9.054741  7.636498  8.678515  8.200096  8.767222  8.219560
      [3,] 13.039802  1.696952  3.759262  9.809671 13.098795 -1.065921
      [4,] 25.466085 18.326098  9.691086  4.186567  1.483841 11.843660
      [5,] 19.571931 11.780854 10.197163  8.152090  9.837362 11.410876
      [6,] 11.027558  5.586653  9.651510  7.027319  6.748893 12.877602
          
iterations      [,7]       [,8]
      [1,] 16.599195 19.8205363
      [2,]  8.049113  8.6893515
      [3,]  5.376531 -0.1256793
      [4,] 15.633438  2.6258156
      [5,] 19.492270  1.3516516
      [6,] 16.822057  5.9936361


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:

 [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__"    
 [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__"    
$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:

[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:

          parameters
iterations         mu  theta[1]
      [1,]  8.4438307  7.354780
      [2,] -0.8464116  3.003915
      [3,]  6.2475111  9.006303
      [4,]  6.7778817  6.760668
      [5,] 11.1856605 11.255970
      [6,] 10.4500162 10.216112


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         8.152373213 0.11286510 5.1358386  -1.393623   4.7579790
tau        6.597015226 0.16025384 5.8470179   0.182007   2.3324600
eta[1]     0.375105745 0.01521985 0.9423738  -1.539866  -0.2349581
eta[2]     0.004162496 0.01533488 0.8969768  -1.743882  -0.5786994
eta[3]    -0.200013365 0.01569546 0.9320128  -1.958753  -0.8200778
eta[4]    -0.041432261 0.01457409 0.8762457  -1.767024  -0.6309151
eta[5]    -0.363582954 0.01624522 0.9065110  -2.147410  -0.9557522
eta[6]    -0.224676543 0.01546635 0.8838848  -1.984021  -0.7946552
eta[7]     0.330342622 0.01516175 0.9056724  -1.477945  -0.2371773
eta[8]     0.037542729 0.01527968 0.9381113  -1.863492  -0.5866037
theta[1]  11.573989322 0.18028290 8.3234373  -2.026758   6.1323593
theta[2]   8.084274400 0.10084032 6.4022741  -4.421910   4.1664885
theta[3]   6.213624955 0.13581285 7.5143841 -11.118186   2.3809189
theta[4]   7.670352873 0.10107076 6.5205295  -5.721917   3.7855647
theta[5]   5.104999308 0.10547928 6.5009347  -8.850779   1.2159498
theta[6]   6.162412608 0.11207592 6.8322897  -8.568943   2.3651420
theta[7]  10.773490319 0.12874520 7.0191371  -1.067950   5.9848201
theta[8]   8.560590595 0.13489542 7.8117042  -7.062939   4.1957742
lp__     -39.659879168 0.08250758 2.7173736 -45.778398 -41.2920493
                   50%         75%      97.5%    n_eff      Rhat
mu         7.960462031  11.3173363  18.801856 2070.634 1.0014977
tau        5.202786669   9.0768754  21.980606 1331.227 1.0023140
eta[1]     0.406653997   1.0027417   2.189448 3833.767 1.0002345
eta[2]    -0.001200586   0.5974635   1.812073 3421.381 1.0010783
eta[3]    -0.229506631   0.4032847   1.678803 3526.108 1.0000286
eta[4]    -0.053394801   0.5349828   1.687687 3614.839 0.9995517
eta[5]    -0.380153750   0.2156665   1.439561 3113.831 0.9999216
eta[6]    -0.216989555   0.3607663   1.536258 3265.996 0.9994409
eta[7]     0.334079525   0.8977536   2.115716 3568.156 0.9992701
eta[8]     0.030370859   0.6495645   1.871860 3769.471 1.0002459
theta[1]  10.259084055  16.0305452  31.492743 2131.554 1.0002961
theta[2]   7.882515381  11.9434375  21.265143 4030.882 1.0004209
theta[3]   6.844580685  10.6571860  19.964577 3061.294 1.0007327
theta[4]   7.674786247  11.7692145  20.694074 4162.121 0.9992592
theta[5]   5.748314773   9.3156138  16.621352 3798.545 0.9994299
theta[6]   6.486363702  10.3469681  18.692734 3716.275 0.9998919
theta[7]  10.024560826  14.7663594  26.756726 2972.387 1.0007889
theta[8]   8.190506095  12.5719631  25.825598 3353.491 1.0005303
lp__     -39.386757916 -37.7745133 -35.060336 1084.703 1.0019144

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  8.152373 0.1128651 5.135839 1.9258178 14.62122 2070.634 1.001498
tau 6.597015 0.1602538 5.847018 0.8341973 13.77241 1331.227 1.002314

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.9258178 14.62122
tau 0.8341973 13.77241


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.7495132 0.8524633 0.8615635 0.9528063
max_treedepth_by_chain <- sapply(sampler_params, function(x) max(x[, "treedepth__"]))
print(max_treedepth_by_chain)
[1] 6 5 5 5


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.7013894

$tau
[1] 1.791801

$eta
[1]  1.19252640  1.06597020 -1.22324298  0.02262201  1.57796061  1.94667109
[7] -1.99503156 -1.50724646

$theta
[1]  1.4353801  1.2086166 -2.8931967 -0.6608552  2.1260013  2.7866569 -4.2760879
[8] -3.4020743


(P)RNG seed

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

print(get_seed(fit))
[1] 2081004574


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.030285 0.026706
chain:2 0.032524 0.028686
chain:3 0.031837 0.028443
chain:4 0.034213 0.036218