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)

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]  1.636743 10.106585  3.513336 11.023126  8.338745 17.326182

head(list_of_draws$tau)

[1]  9.572745  1.990020  5.071862  5.896383 21.024283  2.244242

head(list_of_draws$theta)

          
iterations       [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
      [1,]  1.8772295 11.129563  1.580120 -4.820752  2.910281 -3.193219
      [2,] 11.1700082 10.736322  7.949672  9.920369  7.585557 10.378372
      [3,]  0.5476151  6.849613  2.806382  3.297504 12.099958  8.499055
      [4,] 15.9757297 16.946180 16.896805 18.205883  5.337248 10.420042
      [5,] 28.6057365 12.585173 -5.922589 -1.343499 11.377180 22.199806
      [6,] 21.3575933 15.047661 16.108431 16.484844 13.623674 14.455254
          
iterations      [,7]      [,8]
      [1,] 16.759992  9.074250
      [2,]  9.831815 13.308429
      [3,]  5.489958 -6.017911
      [4,] 13.393023 18.655358
      [5,]  6.437247  3.896331
      [6,] 17.731023 15.110924

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,]  7.935440  9.625193
      [2,]  6.636731  7.714632
      [3,]  2.077142  5.929779
      [4,] 13.197995 11.888387
      [5,] 11.083340  6.367072
      [6,] 10.186964 10.764070

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.080801408 0.11205740 5.2075056  -2.0749510   4.5926859
tau        6.468872365 0.14314978 5.5361712   0.2182851   2.4211382
eta[1]     0.364483616 0.01589541 0.9397065  -1.5142246  -0.2529054
eta[2]     0.002946661 0.01451296 0.8814964  -1.7588001  -0.5833192
eta[3]    -0.198930820 0.01382374 0.9278652  -2.0142361  -0.8138516
eta[4]    -0.048173578 0.01397762 0.8957642  -1.8170455  -0.6409908
eta[5]    -0.365129872 0.01415966 0.8797499  -2.0837888  -0.9449596
eta[6]    -0.212589063 0.01578280 0.9153735  -1.9886084  -0.8376942
eta[7]     0.328016836 0.01455812 0.9175026  -1.5678663  -0.2498630
eta[8]     0.046709776 0.01409150 0.9253261  -1.8179771  -0.5679017
theta[1]  11.296863420 0.15946130 8.2286601  -1.8944536   5.9899087
theta[2]   7.936680090 0.09610888 6.4394188  -4.8882366   3.8449608
theta[3]   6.132254370 0.13245875 7.8714581 -11.9830147   2.0490261
theta[4]   7.695826305 0.10425387 6.6798046  -5.6843890   3.5682426
theta[5]   5.137745151 0.10251973 6.4321704  -9.2304525   1.4408935
theta[6]   6.392899799 0.10370934 6.7181595  -8.4261449   2.5326189
theta[7]  10.658276186 0.11513742 6.9407663  -2.0018650   6.1114540
theta[8]   8.444759974 0.13676735 7.6950799  -6.5072939   3.9432244
lp__     -39.664060295 0.07628426 2.7152740 -45.8245041 -41.2916710
                   50%         75%      97.5%    n_eff      Rhat
mu         7.999463336  11.4383465  18.609302 2159.625 1.0012826
tau        5.182528278   8.8589141  20.108089 1495.677 1.0004633
eta[1]     0.396314473   0.9967559   2.144416 3494.949 0.9996598
eta[2]    -0.002458377   0.5935754   1.717753 3689.173 0.9999393
eta[3]    -0.213413961   0.4125764   1.627728 4505.245 0.9996014
eta[4]    -0.051268014   0.5497469   1.694867 4106.966 0.9997956
eta[5]    -0.366101365   0.1993965   1.384328 3860.228 1.0003907
eta[6]    -0.212427122   0.4001244   1.570076 3363.789 1.0013953
eta[7]     0.342566320   0.9224976   2.114650 3971.950 0.9997084
eta[8]     0.050910504   0.6819199   1.834986 4311.964 0.9996240
theta[1]  10.361269440  15.3604216  31.235137 2662.856 1.0000660
theta[2]   7.966544858  12.0469553  20.726894 4489.173 1.0000285
theta[3]   6.673599321  11.0517115  20.033778 3531.418 0.9995364
theta[4]   7.772971015  11.7123420  21.177812 4105.283 1.0001019
theta[5]   5.569444083   9.4745957  16.695643 3936.409 1.0008361
theta[6]   6.764191552  10.5941905  18.886780 4196.284 0.9997453
theta[7]  10.190632594  14.4894225  26.060456 3633.974 0.9995850
theta[8]   8.252979095  12.8373335  25.323362 3165.639 0.9993903
lp__     -39.359995581 -37.7498956 -34.961036 1266.944 1.0007259

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.080801 0.1120574 5.207506 1.582639 14.43934 2159.625 1.001283
tau 6.468872 0.1431498 5.536171 1.033303 13.80869 1495.677 1.000463

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.582639 14.43934
tau 1.033303 13.80869

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.8871174 0.9019111 0.8683277 0.8614057

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

$tau
[1] 0.9667098

$eta
[1] -1.2508485  0.7433085 -0.3798570 -1.0922688 -0.7170482 -0.2388567 -1.5208663
[8]  0.1987448

$theta
[1] -2.6561980 -0.7284269 -1.8142020 -2.5028974 -2.1401680 -1.6778956 -2.9172269
[8] -1.2548620

(P)RNG seed

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

print(get_seed(fit))

[1] 1899244817

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.024  0.027
chain:2  0.024  0.029
chain:3  0.023  0.021
chain:4  0.025  0.023