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

    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] 14.143929  8.478591  6.504845 11.970911  8.524122  8.160192

    head(list_of_draws$tau)

    [1]  4.3662195  0.4272034  4.1657276 12.2439305  4.3656120  4.5503272

    head(list_of_draws$theta)

    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:

    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,]  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