In this tutorial, we will have a closer look at the object returned by calling run_mcmc()
on a network model. Understanding the structure of the object containing the results of a run is important for model diagnostics and interpretation.
To quickly obtain an MCMC run output for us to examine, let’s run the simple model aquarium_mod
which is provided with the package. Feel free to read the help ?aquarium_mod
if you are curious about the model itself.
library(isotracer)
aquarium_mod run_mcmc(aquarium_mod, iter = 1000) fit <-
By default, run_mcmc()
returns an mcmc.list
object. An mcmc.list
has a simple format to store the content of parallel MCMC chains:
length(fit)
## [1] 4
str(fit[[1]])
## 'mcmc' num [1:500, 1:8] 0.1076 0.1001 0.0879 0.0999 0.083 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:8] "eta" "lambda_algae" "lambda_daphnia" "lambda_NH4" ...
## - attr(*, "mcpar")= num [1:3] 501 1000 1
The output above can seem a little obscure if you are not familiar with R data structures, but in a nutshell it tells us that the mcmc.list
is basically a list with one element per chain, each chain being stored as a matrix.
The mcmc.list
class is implemented by the coda package, and it has the advantage of being recognized by many other R packages dealing with Bayesian MCMC such as bayesplot of ggmcmc.
In the isotracer package, the returned fit
is very slightly extended compared to the base mcmc.list
class:
class(fit)
## [1] "networkModelStanfit" "mcmc.list"
By having also a networkModelStanfit
class, the output from run_mcmc()
can be recognized automatically by some methods implemented in isotracer, such as plot()
:
plot(fit)
# Note: the figure below only shows a few of the traceplots for vignette concision
An mcmc.list
object can be converted to an even simpler, flat matrix:
as.matrix(fit)
z <-head(z)
## eta lambda_algae lambda_daphnia lambda_NH4 upsilon_algae_to_daphnia
## [1,] 0.10760688 0.05522696 0.02610060 0.05495770 0.06229580
## [2,] 0.10006176 0.08957437 0.05246263 0.09519854 0.05928665
## [3,] 0.08794362 0.09731864 0.01296905 0.05948244 0.08819510
## [4,] 0.09988717 0.14601943 0.02732369 0.12801800 0.08151178
## [5,] 0.08301680 0.14333069 0.02702577 0.17689506 0.05724034
## [6,] 0.08218898 0.10842972 0.04077288 0.16627876 0.06871787
## upsilon_daphnia_to_NH4 upsilon_NH4_to_algae zeta
## [1,] 0.04736973 0.3805771 0.5649217
## [2,] 0.04795699 0.3789833 0.5012341
## [3,] 0.05089449 0.3055059 0.2947550
## [4,] 0.05185232 0.3143063 0.3159154
## [5,] 0.04046591 0.4112083 0.7108229
## [6,] 0.04015712 0.4102410 0.4765276
str(z)
## num [1:2000, 1:8] 0.1076 0.1001 0.0879 0.0999 0.083 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:8] "eta" "lambda_algae" "lambda_daphnia" "lambda_NH4" ...
or to a data frame:
as.data.frame(as.matrix(fit))
z <-head(z)
## eta lambda_algae lambda_daphnia lambda_NH4 upsilon_algae_to_daphnia
## 1 0.10760688 0.05522696 0.02610060 0.05495770 0.06229580
## 2 0.10006176 0.08957437 0.05246263 0.09519854 0.05928665
## 3 0.08794362 0.09731864 0.01296905 0.05948244 0.08819510
## 4 0.09988717 0.14601943 0.02732369 0.12801800 0.08151178
## 5 0.08301680 0.14333069 0.02702577 0.17689506 0.05724034
## 6 0.08218898 0.10842972 0.04077288 0.16627876 0.06871787
## upsilon_daphnia_to_NH4 upsilon_NH4_to_algae zeta
## 1 0.04736973 0.3805771 0.5649217
## 2 0.04795699 0.3789833 0.5012341
## 3 0.05089449 0.3055059 0.2947550
## 4 0.05185232 0.3143063 0.3159154
## 5 0.04046591 0.4112083 0.7108229
## 6 0.04015712 0.4102410 0.4765276
str(z)
## 'data.frame': 2000 obs. of 8 variables:
## $ eta : num 0.1076 0.1001 0.0879 0.0999 0.083 ...
## $ lambda_algae : num 0.0552 0.0896 0.0973 0.146 0.1433 ...
## $ lambda_daphnia : num 0.0261 0.0525 0.013 0.0273 0.027 ...
## $ lambda_NH4 : num 0.055 0.0952 0.0595 0.128 0.1769 ...
## $ upsilon_algae_to_daphnia: num 0.0623 0.0593 0.0882 0.0815 0.0572 ...
## $ upsilon_daphnia_to_NH4 : num 0.0474 0.048 0.0509 0.0519 0.0405 ...
## $ upsilon_NH4_to_algae : num 0.381 0.379 0.306 0.314 0.411 ...
## $ zeta : num 0.565 0.501 0.295 0.316 0.711 ...
or to a tibble:
tibble::as_tibble(as.matrix(fit))
z <- z
## # A tibble: 2,000 × 8
## eta lambda_algae lambda_daphnia lambda_NH4 upsilon_algae_to_d… upsilon_daphnia…
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.108 0.0552 0.0261 0.0550 0.0623 0.0474
## 2 0.100 0.0896 0.0525 0.0952 0.0593 0.0480
## 3 0.0879 0.0973 0.0130 0.0595 0.0882 0.0509
## 4 0.0999 0.146 0.0273 0.128 0.0815 0.0519
## 5 0.0830 0.143 0.0270 0.177 0.0572 0.0405
## 6 0.0822 0.108 0.0408 0.166 0.0687 0.0402
## 7 0.0815 0.0314 0.0882 0.0959 0.0568 0.0541
## 8 0.0765 0.0879 0.0630 0.0638 0.0432 0.0493
## 9 0.181 0.170 0.0999 0.127 0.0830 0.0331
## 10 0.119 0.0924 0.0156 0.196 0.0794 0.0339
## # … with 1,990 more rows, and 2 more variables: upsilon_NH4_to_algae <dbl>,
## # zeta <dbl>
Converting your output to one of those simple tabular formats can be useful if you want to manipulate and perform operations on your MCMC samples.
However, for simple manipulations, isotracer
provides convenient methods to perform calculations on parameter chains directly from the output of run_mcmc()
. You can thus produce derived parameter chains directly from the mcmc.list
object, without having to convert your output to another format:
fit[, "upsilon_algae_to_daphnia"] + fit[, "lambda_algae"]
algal_total_out <- 1 / algal_total_out
algal_turnover <-plot(algal_turnover)
You can read more about this in the vignette about calculating derived parameters.
You can combine derived parameters into a single mcmc.list
object using the usual c()
syntax. This can be convenient for more compact plotting or summary calculations:
c("out rate" = algal_total_out, "turnover" = algal_turnover)
my_derived <-plot(my_derived)
summary(my_derived)
##
## Iterations = 501:1000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## out rate 0.1666 0.05965 0.001334 0.001945
## turnover 6.9019 2.83356 0.063360 0.126940
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## out rate 0.06908 0.123 0.1615 0.2033 0.3012
## turnover 3.32046 4.920 6.1908 8.1333 14.4769
Calling run_mcmc()
will run a Stan model behind the scenes. Stan is great since it will let you know loudly when something went wrong with the run, such as problems with divergent chains or low Bayesian fraction of missing information. Such problems should not be ignored! The Stan development team has a nice page explaining Stan’s warnings.
In any case, if something went wrong with your run, you might want to have a more complete output than simply the mcmc.list
object. You can ask run_mcmc()
to return the original stanfit
object produced by Stan with:
run_mcmc(aquarium_mod, iter = 1000, stanfit = TRUE) fit2 <-
fit2
is now a regular stanfit
object:
class(fit2)
## [1] "stanfit"
## attr(,"package")
## [1] "rstan"
This is a more complicated type of object than an mcmc.list
, but it also contains much more information about the Stan run. It also comes with the benefit of the existing methods for stanfit
object, for example:
::plot(fit2) rstan
You can go through Stan documentation for more details about this format. If you are reading about solving Stan model issues on online forums and the suggested solutions require to examine some Stan output, that’s the object you want to look at!
For example, you can examine it with ShinyStan:
library(shinystan)
launch_shinystan(fit2)
Note that with the current version of isotracer the parameters are indexed but not named in the stanfit
object. That is something that will probably be improved in the future!