Summary method for stanreg objects

Source: R/print-and-summary.R

Summaries of parameter estimates and MCMC convergence diagnostics (Monte Carlo error, effective sample size, Rhat).

# S3 method for stanreg
summary(
  object,
  pars = NULL,
  regex_pars = NULL,
  probs = c(0.1, 0.5, 0.9),
  ...,
  digits = 1
)

# S3 method for summary.stanreg
print(x, digits = max(1, attr(x, "print.digits")), ...)

# S3 method for summary.stanreg
as.data.frame(x, ...)

# S3 method for stanmvreg
summary(object, pars = NULL, regex_pars = NULL, probs = NULL, ..., digits = 3)

# S3 method for summary.stanmvreg
print(x, digits = max(1, attr(x, "print.digits")), ...)

Arguments

object

A fitted model object returned by one of the rstanarm modeling functions. See stanreg-objects.
pars

An optional character vector specifying a subset of parameters to display. Parameters can be specified by name or several shortcuts can be used. Using pars="beta" will restrict the displayed parameters to only the regression coefficients (without the intercept). "alpha" can also be used as a shortcut for "(Intercept)". If the model has varying intercepts and/or slopes they can be selected using pars = "varying".

In addition, for stanmvreg objects there are some additional shortcuts available. Using pars = "long" will display the parameter estimates for the longitudinal submodels only (excluding group-specific pparameters, but including auxiliary parameters). Using pars = "event" will display the parameter estimates for the event submodel only, including any association parameters. Using pars = "assoc" will display only the association parameters. Using pars = "fixef" will display all fixed effects, but not the random effects or the auxiliary parameters. pars and regex_pars are set to NULL then all fixed effect regression coefficients are selected, as well as any auxiliary parameters and the log posterior.

If pars is NULL all parameters are selected for a stanreg object, while for a stanmvreg object all fixed effect regression coefficients are selected as well as any auxiliary parameters and the log posterior. See Examples.
regex_pars

An optional character vector of regular expressions to use for parameter selection. regex_pars can be used in place of pars or in addition to pars. Currently, all functions that accept a regex_pars argument ignore it for models fit using optimization.
probs

For models fit using MCMC or one of the variational algorithms, an optional numeric vector of probabilities passed to quantile.
...

Currently ignored.
digits

Number of digits to use for formatting numbers when printing. When calling summary, the value of digits is stored as the "print.digits" attribute of the returned object.
x

An object of class "summary.stanreg".

Value

The summary method returns an object of class "summary.stanreg" (or "summary.stanmvreg", inheriting "summary.stanreg"), which is a matrix of summary statistics and diagnostics, with attributes storing information for use by the print method. The print method for summary.stanreg or summary.stanmvreg objects is called for its side effect and just returns its input. The as.data.frame method for summary.stanreg objects converts the matrix to a data.frame, preserving row and column names but dropping the print-related attributes.

Details

mean_PPD diagnostic

Summary statistics are also reported for mean_PPD, the sample average posterior predictive distribution of the outcome. This is useful as a quick diagnostic. A useful heuristic is to check if mean_PPD is plausible when compared to mean(y). If it is plausible then this does not mean that the model is good in general (only that it can reproduce the sample mean), however if mean_PPD is implausible then it is a sign that something is wrong (severe model misspecification, problems with the data, computational issues, etc.).

See also

prior_summary to extract or print a summary of the priors used for a particular model.

Examples

if (!exists("example_model")) example(example_model)
summary(example_model, probs = c(0.1, 0.9))
#> 
#> Model Info:
#>  function:     stan_glmer
#>  family:       binomial [logit]
#>  formula:      cbind(incidence, size - incidence) ~ size + period + (1 | herd)
#>  algorithm:    sampling
#>  sample:       1000 (posterior sample size)
#>  priors:       see help('prior_summary')
#>  observations: 56
#>  groups:       herd (15)
#> 
#> Estimates:
#>                                       mean   sd   10%   90%
#> (Intercept)                         -1.5    0.6 -2.3  -0.8 
#> size                                 0.0    0.0  0.0   0.0 
#> period2                             -1.0    0.3 -1.4  -0.6 
#> period3                             -1.1    0.3 -1.6  -0.7 
#> period4                             -1.6    0.5 -2.2  -1.0 
#> b[(Intercept) herd:1]                0.6    0.4  0.1   1.2 
#> b[(Intercept) herd:2]               -0.4    0.4 -0.9   0.1 
#> b[(Intercept) herd:3]                0.4    0.4 -0.1   0.8 
#> b[(Intercept) herd:4]                0.0    0.5 -0.6   0.6 
#> b[(Intercept) herd:5]               -0.3    0.4 -0.8   0.2 
#> b[(Intercept) herd:6]               -0.5    0.4 -1.1   0.0 
#> b[(Intercept) herd:7]                0.9    0.4  0.4   1.4 
#> b[(Intercept) herd:8]                0.5    0.5 -0.1   1.1 
#> b[(Intercept) herd:9]               -0.3    0.5 -1.0   0.3 
#> b[(Intercept) herd:10]              -0.7    0.5 -1.3  -0.1 
#> b[(Intercept) herd:11]              -0.2    0.4 -0.7   0.3 
#> b[(Intercept) herd:12]              -0.1    0.5 -0.8   0.6 
#> b[(Intercept) herd:13]              -0.8    0.5 -1.5  -0.3 
#> b[(Intercept) herd:14]               1.0    0.5  0.4   1.6 
#> b[(Intercept) herd:15]              -0.6    0.5 -1.3  -0.1 
#> Sigma[herd:(Intercept),(Intercept)]  0.6    0.4  0.2   1.0 
#> 
#> Fit Diagnostics:
#>            mean   sd   10%   90%
#> mean_PPD 1.8    0.2  1.5   2.1  
#> 
#> The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
#> 
#> MCMC diagnostics
#>                                     mcse Rhat n_eff
#> (Intercept)                         0.0  1.0   692 
#> size                                0.0  1.0   668 
#> period2                             0.0  1.0  1653 
#> period3                             0.0  1.0  1192 
#> period4                             0.0  1.0   867 
#> b[(Intercept) herd:1]               0.0  1.0   877 
#> b[(Intercept) herd:2]               0.0  1.0   850 
#> b[(Intercept) herd:3]               0.0  1.0   569 
#> b[(Intercept) herd:4]               0.0  1.0   910 
#> b[(Intercept) herd:5]               0.0  1.0   985 
#> b[(Intercept) herd:6]               0.0  1.0   804 
#> b[(Intercept) herd:7]               0.0  1.0   695 
#> b[(Intercept) herd:8]               0.0  1.0   817 
#> b[(Intercept) herd:9]               0.0  1.0   853 
#> b[(Intercept) herd:10]              0.0  1.0   698 
#> b[(Intercept) herd:11]              0.0  1.0   681 
#> b[(Intercept) herd:12]              0.0  1.0   957 
#> b[(Intercept) herd:13]              0.0  1.0   761 
#> b[(Intercept) herd:14]              0.0  1.0   936 
#> b[(Intercept) herd:15]              0.0  1.0   951 
#> Sigma[herd:(Intercept),(Intercept)] 0.0  1.0   424 
#> mean_PPD                            0.0  1.0  1115 
#> log-posterior                       0.3  1.0   290 
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

# These produce the same output for this example, 
# but the second method can be used for any model
summary(example_model, pars = c("(Intercept)", "size",
                                paste0("period", 2:4)))
#> 
#> Model Info:
#>  function:     stan_glmer
#>  family:       binomial [logit]
#>  formula:      cbind(incidence, size - incidence) ~ size + period + (1 | herd)
#>  algorithm:    sampling
#>  sample:       1000 (posterior sample size)
#>  priors:       see help('prior_summary')
#>  observations: 56
#>  groups:       herd (15)
#> 
#> Estimates:
#>               mean   sd   10%   50%   90%
#> (Intercept) -1.5    0.6 -2.3  -1.5  -0.8 
#> size         0.0    0.0  0.0   0.0   0.0 
#> period2     -1.0    0.3 -1.4  -1.0  -0.6 
#> period3     -1.1    0.3 -1.6  -1.1  -0.7 
#> period4     -1.6    0.5 -2.2  -1.5  -1.0 
#> 
#> MCMC diagnostics
#>             mcse Rhat n_eff
#> (Intercept) 0.0  1.0   692 
#> size        0.0  1.0   668 
#> period2     0.0  1.0  1653 
#> period3     0.0  1.0  1192 
#> period4     0.0  1.0   867 
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(example_model, pars = c("alpha", "beta"))
#> 
#> Model Info:
#>  function:     stan_glmer
#>  family:       binomial [logit]
#>  formula:      cbind(incidence, size - incidence) ~ size + period + (1 | herd)
#>  algorithm:    sampling
#>  sample:       1000 (posterior sample size)
#>  priors:       see help('prior_summary')
#>  observations: 56
#>  groups:       herd (15)
#> 
#> Estimates:
#>               mean   sd   10%   50%   90%
#> (Intercept) -1.5    0.6 -2.3  -1.5  -0.8 
#> size         0.0    0.0  0.0   0.0   0.0 
#> period2     -1.0    0.3 -1.4  -1.0  -0.6 
#> period3     -1.1    0.3 -1.6  -1.1  -0.7 
#> period4     -1.6    0.5 -2.2  -1.5  -1.0 
#> 
#> MCMC diagnostics
#>             mcse Rhat n_eff
#> (Intercept) 0.0  1.0   692 
#> size        0.0  1.0   668 
#> period2     0.0  1.0  1653 
#> period3     0.0  1.0  1192 
#> period4     0.0  1.0   867 
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

# Only show parameters varying by group
summary(example_model, pars = "varying")
#> 
#> Model Info:
#>  function:     stan_glmer
#>  family:       binomial [logit]
#>  formula:      cbind(incidence, size - incidence) ~ size + period + (1 | herd)
#>  algorithm:    sampling
#>  sample:       1000 (posterior sample size)
#>  priors:       see help('prior_summary')
#>  observations: 56
#>  groups:       herd (15)
#> 
#> Estimates:
#>                          mean   sd   10%   50%   90%
#> b[(Intercept) herd:1]   0.6    0.4  0.1   0.6   1.2 
#> b[(Intercept) herd:2]  -0.4    0.4 -0.9  -0.4   0.1 
#> b[(Intercept) herd:3]   0.4    0.4 -0.1   0.4   0.8 
#> b[(Intercept) herd:4]   0.0    0.5 -0.6   0.0   0.6 
#> b[(Intercept) herd:5]  -0.3    0.4 -0.8  -0.3   0.2 
#> b[(Intercept) herd:6]  -0.5    0.4 -1.1  -0.5   0.0 
#> b[(Intercept) herd:7]   0.9    0.4  0.4   0.9   1.4 
#> b[(Intercept) herd:8]   0.5    0.5 -0.1   0.5   1.1 
#> b[(Intercept) herd:9]  -0.3    0.5 -1.0  -0.3   0.3 
#> b[(Intercept) herd:10] -0.7    0.5 -1.3  -0.6  -0.1 
#> b[(Intercept) herd:11] -0.2    0.4 -0.7  -0.2   0.3 
#> b[(Intercept) herd:12] -0.1    0.5 -0.8  -0.1   0.6 
#> b[(Intercept) herd:13] -0.8    0.5 -1.5  -0.8  -0.3 
#> b[(Intercept) herd:14]  1.0    0.5  0.4   1.0   1.6 
#> b[(Intercept) herd:15] -0.6    0.5 -1.3  -0.6  -0.1 
#> 
#> MCMC diagnostics
#>                        mcse Rhat n_eff
#> b[(Intercept) herd:1]  0.0  1.0  877  
#> b[(Intercept) herd:2]  0.0  1.0  850  
#> b[(Intercept) herd:3]  0.0  1.0  569  
#> b[(Intercept) herd:4]  0.0  1.0  910  
#> b[(Intercept) herd:5]  0.0  1.0  985  
#> b[(Intercept) herd:6]  0.0  1.0  804  
#> b[(Intercept) herd:7]  0.0  1.0  695  
#> b[(Intercept) herd:8]  0.0  1.0  817  
#> b[(Intercept) herd:9]  0.0  1.0  853  
#> b[(Intercept) herd:10] 0.0  1.0  698  
#> b[(Intercept) herd:11] 0.0  1.0  681  
#> b[(Intercept) herd:12] 0.0  1.0  957  
#> b[(Intercept) herd:13] 0.0  1.0  761  
#> b[(Intercept) herd:14] 0.0  1.0  936  
#> b[(Intercept) herd:15] 0.0  1.0  951  
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
as.data.frame(summary(example_model, pars = "varying"))
#>                                mean       mcse        sd         10%
#> b[(Intercept) herd:1]   0.604783025 0.01471882 0.4359469  0.06148133
#> b[(Intercept) herd:2]  -0.382994268 0.01462581 0.4264550 -0.91458443
#> b[(Intercept) herd:3]   0.370266800 0.01516369 0.3618129 -0.09723500
#> b[(Intercept) herd:4]   0.005244565 0.01594533 0.4810427 -0.58123388
#> b[(Intercept) herd:5]  -0.289911929 0.01328420 0.4168151 -0.82669390
#> b[(Intercept) herd:6]  -0.527324356 0.01579077 0.4476147 -1.09942183
#> b[(Intercept) herd:7]   0.902459532 0.01627785 0.4290710  0.37536777
#> b[(Intercept) herd:8]   0.513016497 0.01767929 0.5053240 -0.06661935
#> b[(Intercept) herd:9]  -0.295294019 0.01809977 0.5287524 -0.98924943
#> b[(Intercept) herd:10] -0.667694946 0.01776687 0.4693281 -1.29119377
#> b[(Intercept) herd:11] -0.179119207 0.01607633 0.4194779 -0.70084966
#> b[(Intercept) herd:12] -0.094495552 0.01732486 0.5359092 -0.78562181
#> b[(Intercept) herd:13] -0.846300347 0.01770006 0.4882524 -1.51910581
#> b[(Intercept) herd:14]  0.999953352 0.01479055 0.4525484  0.43870004
#> b[(Intercept) herd:15] -0.645605381 0.01594687 0.4916736 -1.30940428
#>                                 50%         90% n_eff      Rhat
#> b[(Intercept) herd:1]   0.601862454  1.16625406   877 1.0014746
#> b[(Intercept) herd:2]  -0.360913301  0.13469303   850 1.0011829
#> b[(Intercept) herd:3]   0.383693821  0.81024339   569 0.9994316
#> b[(Intercept) herd:4]  -0.007906003  0.59977894   910 0.9985698
#> b[(Intercept) herd:5]  -0.269964967  0.23529695   985 1.0036898
#> b[(Intercept) herd:6]  -0.498233143  0.01616496   804 0.9992756
#> b[(Intercept) herd:7]   0.881424738  1.43156904   695 0.9991341
#> b[(Intercept) herd:8]   0.507587177  1.14087504   817 0.9985972
#> b[(Intercept) herd:9]  -0.256824944  0.32911594   853 1.0016025
#> b[(Intercept) herd:10] -0.641624524 -0.08721296   698 1.0047409
#> b[(Intercept) herd:11] -0.164954759  0.32468865   681 1.0017397
#> b[(Intercept) herd:12] -0.069369761  0.56870495   957 1.0010215
#> b[(Intercept) herd:13] -0.791944334 -0.27120170   761 1.0024421
#> b[(Intercept) herd:14]  0.981693175  1.59229011   936 1.0008012
#> b[(Intercept) herd:15] -0.606239187 -0.05207457   951 0.9986671