The default format: mcmc.list
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
fit <- run_mcmc(aquarium_mod, iter = 1000)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] 4str(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 1The 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
How to convert the default output to a table?
An mcmc.list object can be converted to an even simpler, flat matrix:
z <- as.matrix(fit)
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.4765276str(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:
z <- as.data.frame(as.matrix(fit))
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.4765276str(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:
z <- tibble::as_tibble(as.matrix(fit))
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>How to calculate derived parameters?
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:
algal_total_out <- fit[, "upsilon_algae_to_daphnia"] + fit[, "lambda_algae"]
algal_turnover <- 1 / algal_total_out
plot(algal_turnover)
You can read more about this in the vignette about calculating derived parameters.
How to combine derived parameters?
You can combine derived parameters into a single mcmc.listobject using the usual c() syntax. This can be convenient for more compact plotting or summary calculations:
my_derived <- c("out rate" = algal_total_out, "turnover" = algal_turnover)
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