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]  3.4728672  7.4003441 -0.2154997  6.6868831 -0.3536323  7.6920400

head(list_of_draws$tau)

[1] 0.6237294 8.9290545 8.2189972 2.6435660 7.4717684 4.7409298

head(list_of_draws$theta)

          
iterations      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
      [1,]  2.392530  3.644290  3.769977  3.238493  3.901345  3.862098
      [2,] 11.253007 21.237775 -4.567886  7.363705  5.836193 -5.177758
      [3,] 16.285581  8.200160 -5.122407 -1.033780 -4.245903  4.161598
      [4,] 10.768829  7.812361  4.076096  2.966960  8.573051  1.156534
      [5,]  2.970755 -3.206895 -7.676241 11.804634  6.190454 -3.529938
      [6,] 11.348168 10.103474  6.960276 13.159084 11.674984  6.815160
          
iterations      [,7]     [,8]
      [1,]  3.673559 2.801205
      [2,]  6.936138 8.452885
      [3,]  6.050967 9.368861
      [4,]  5.536694 1.646799
      [5,]  5.327677 2.050719
      [6,] 12.373056 7.667725

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,] 10.537607 10.522905
      [2,]  2.938166 29.334814
      [3,] -4.615421  8.238651
      [4,]  4.535026 23.065188
      [5,]  2.819068 15.310006
      [6,] 10.908570  6.805632

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.920884945 0.10274229 5.1140662  -2.3401189   4.7003808
tau        6.544862144 0.15080833 5.6565933   0.2736641   2.5736824
eta[1]     0.375819859 0.01378907 0.9387065  -1.5278079  -0.2380295
eta[2]     0.002680126 0.01334088 0.8832060  -1.7400620  -0.5692828
eta[3]    -0.201181770 0.01390818 0.9368744  -2.0626873  -0.8239312
eta[4]    -0.040855277 0.01316900 0.8836881  -1.8101482  -0.6128079
eta[5]    -0.342098614 0.01499935 0.8772424  -2.0185329  -0.9185974
eta[6]    -0.224052232 0.01341102 0.9039113  -2.0275318  -0.8243524
eta[7]     0.349908532 0.01401329 0.8762909  -1.4378839  -0.2209485
eta[8]     0.037580943 0.01368377 0.9324548  -1.8193377  -0.5997747
theta[1]  11.392648244 0.14496326 8.4841745  -2.5758970   5.9964760
theta[2]   7.973278662 0.09627116 6.3227710  -4.4708498   3.8679721
theta[3]   6.227297739 0.11909702 7.7401863 -10.6180879   1.9880647
theta[4]   7.515028977 0.09028588 6.5144721  -6.0853936   3.4975795
theta[5]   5.192334344 0.09608942 6.3468044  -8.9111106   1.2921103
theta[6]   6.104247312 0.10037443 6.7633929  -8.4533432   2.0964797
theta[7]  10.745684739 0.10101460 6.8247826  -1.0349905   6.1380260
theta[8]   8.319483021 0.12166044 7.7499806  -7.0897167   3.7634656
lp__     -39.547088378 0.07638554 2.6347329 -45.3724286 -41.1467983
                   50%         75%      97.5%    n_eff      Rhat
mu         7.890351882  11.2164803  18.139187 2477.617 1.0000284
tau        5.299792576   9.0538264  20.052104 1406.887 1.0029480
eta[1]     0.387191645   1.0340843   2.138936 4634.357 0.9996945
eta[2]     0.007291953   0.5722985   1.773130 4382.834 0.9992282
eta[3]    -0.225814402   0.4402071   1.630638 4537.557 0.9996257
eta[4]    -0.044828212   0.5299942   1.706805 4502.904 0.9998282
eta[5]    -0.362917648   0.2141584   1.423233 3420.540 0.9995929
eta[6]    -0.240951928   0.3588670   1.580491 4542.844 0.9992013
eta[7]     0.374800559   0.9454167   1.978462 3910.357 1.0002004
eta[8]     0.049685654   0.6668390   1.806490 4643.488 1.0000520
theta[1]  10.260849371  15.7546436  31.981301 3425.337 1.0015243
theta[2]   7.928647134  12.0401938  20.586894 4313.427 0.9999035
theta[3]   6.663974701  10.9799013  20.479693 4223.778 0.9996042
theta[4]   7.582698685  11.6440264  20.100192 5206.175 0.9995314
theta[5]   5.506990823   9.5279596  16.474048 4362.738 1.0002054
theta[6]   6.476566877  10.4806817  18.613434 4540.284 0.9994761
theta[7]  10.165768543  14.7294775  26.127286 4564.670 0.9996166
theta[8]   8.058250671  12.5827117  24.795244 4057.910 1.0004818
lp__     -39.227276776 -37.6781744 -35.067125 1189.737 1.0038634

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.920885 0.1027423 5.114066 1.574383 14.15416 2477.617 1.000028
tau 6.544862 0.1508083 5.656593 1.002575 13.28935 1406.887 1.002948

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.574383 14.15416
tau 1.002575 13.28935

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.8698983 0.8916800 0.9125469 0.8893986

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

[1] 4 4 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  array[J] real y;               // estimated treatment effects\n  array[J] real<lower=0> sigma;  // 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
  array[J] real y;               // estimated treatment effects
  array[J] real<lower=0> sigma;  // 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.427139

$tau
[1] 1.408791

$eta
[1]  0.1002901 -1.4719556 -1.4158338  0.6120305  1.7836682 -0.2654251  1.7289236
[8] -0.6173598

$theta
[1]  1.5684269 -0.6465389 -0.5674749  2.2893621  3.9399549  1.0532106  3.8628312
[8]  0.5574081

(P)RNG seed

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

print(get_seed(fit))

[1] 644592279

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.022  0.022
chain:2  0.019  0.021
chain:3  0.020  0.024
chain:4  0.020  0.023