diff --git a/.gitignore b/.gitignore index cfbae8c6e..3fe07d6fe 100644 --- a/.gitignore +++ b/.gitignore @@ -62,3 +62,7 @@ content/**/*.html # Local folders with binary data data_gisclub/ + +# png images for the brms tutorial +/content/tutorials/r_brms/brms_eng/*.png +/content/tutorials/r_brms/brms_nl/*.png diff --git a/content/tutorials/r_brms/brms_eng/bayesian_statistics_1_script_eng.R b/content/tutorials/r_brms/brms_eng/bayesian_statistics_1_script_eng.R index 21e47a43d..21d93568d 100644 --- a/content/tutorials/r_brms/brms_eng/bayesian_statistics_1_script_eng.R +++ b/content/tutorials/r_brms/brms_eng/bayesian_statistics_1_script_eng.R @@ -48,7 +48,7 @@ confint(lm1, level = 0.9) ## ------------------------------------------------------------------------------------------------------------------------------------- -source(file = "./source/mcmc_functions.R") +source(file = here("code", "brms_modeling", "mcmc_functions.R")) ## ------------------------------------------------------------------------------------------------------------------------------------- @@ -230,7 +230,7 @@ fit_normal1 <- brm( family = gaussian(), # we use the Normal distribution data = ants_df, # specify data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -335,7 +335,7 @@ mcmc_neff(ratios_fit_normal1) + yaxis_text(hjust = 1) ## ----simple-model-fit1---------------------------------------------------------------------------------------------------------------- # Visualise model fit via bayesplot package -pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -346,7 +346,7 @@ fit_poisson1 <- brm( family = poisson(), # we use the Poisson distribution data = ants_df, # specify the data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -366,7 +366,7 @@ mcmc_rhat(rhats_fit_poisson1) + yaxis_text(hjust = 1) ## ----poisson-model-fit-vis------------------------------------------------------------------------------------------------------------ # Visualise model fit via bayesplot package -pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -377,7 +377,7 @@ fit_poisson2 <- brm( family = poisson(), data = ants_df, chains = nchains, - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -399,7 +399,7 @@ mcmc_rhat(rhats_fit_poisson2) + yaxis_text(hjust = 1) ## ----rand-intercept-model-fit-vis----------------------------------------------------------------------------------------------------- # Visualise model fit of the Poisson model with random intercept via # bayesplot package -pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -518,7 +518,7 @@ fit_poisson2 %>% # calculate average numbers and convert to long format for visualisation mutate(bog = exp(b_Intercept), forest = exp(b_Intercept + b_habitatForest)) %>% - pivot_longer(cols = c("bog", "forest"), names_to = "habitat", + pivot_longer(cols = c("bog", "forest"), names_to = "habitat", values_to = "sp_rich") %>% # visualise via ggplot() ggplot(aes(y = sp_rich, x = habitat)) + diff --git a/content/tutorials/r_brms/brms_eng/workshop_1_mcmc_en_brms_eng.Rmd b/content/tutorials/r_brms/brms_eng/workshop_1_mcmc_en_brms_eng.Rmd index 2dc1e613b..7e2930c72 100644 --- a/content/tutorials/r_brms/brms_eng/workshop_1_mcmc_en_brms_eng.Rmd +++ b/content/tutorials/r_brms/brms_eng/workshop_1_mcmc_en_brms_eng.Rmd @@ -36,12 +36,18 @@ library(tidybayes) # wrangling and visualising Bayesian models # Conflicten tussen packages conflicted::conflicts_prefer(dplyr::filter) conflicted::conflicts_prefer(dplyr::lag) +conflicted::conflicts_prefer(tidyr::extract) conflicted::conflicts_prefer(brms::ar) conflicted::conflicts_prefer(brms::dstudent_t) conflicted::conflicts_prefer(brms::pstudent_t) conflicted::conflicts_prefer(brms::qstudent_t) conflicted::conflicts_prefer(brms::rstudent_t) conflicted::conflict_prefer("rhat", "brms") + +# reduce html export file size by specifying graphics size +knitr::opts_chunk$set(dev = "jpeg", dpi = 56, + fig.width = 4, fig.height = 3, + out.width = "60%") ``` # Theoretical background @@ -239,7 +245,7 @@ When taking a step uphill, the MCMC always continues. If the step is downhill, t The absolute height of the mountain (posterior likelihood) in itself does not matter. Only the relative differences or shape of the mountain matters. This means that we do not need to calculate the denominator of Bayes' rule, which is the reason why the algorithm is frequently used in Bayesian statistics. -![Metropolis rules](../media/MCMC4b.png) +![Metropolis rules](./mcmc4b.png) A simple MCMC algorithm is the Metropolis algorithm. It is described by these formal rules: @@ -297,7 +303,7 @@ confint(lm1, level = 0.9) We source some short functions to calculate the (log) likelihood and the prior and to execute the MCMC metropolis algorithm. ```{r} -source(file = "./source/mcmc_functions.R") +source(file = here("static", "code", "brms_modeling", "mcmc_functions.R")) ``` For this simple model with a small data set, we can calculate and plot the posterior for a large number of combinations of 'beta_0' and 'beta_1'. @@ -345,7 +351,7 @@ ggplot(df, aes(x = beta_0, y = beta_1, z = post)) + The starting value is quite far from the maximum likelihood. It takes a while for the MCMC to stabilize in the vicinity of the top of the mountain. -Therefore, the first part of the MCMC is never used. This is called the 'burn-in' or 'warm-up'. +Therefore, the first part of the MCMC is never used. This is called the 'warmup', 'burn-in' or 'tuning'. We now run the same model, but with a much longer MCMC (more iterations): @@ -392,7 +398,7 @@ The MCMC for each parameter can be displayed in a so called 'trace plot' where w ann_text <- tibble(param = c("beta_0", "beta_1"), x = 210, y = c(max(mcmc_l$beta_0), max(mcmc_l$beta_1)), - lab = "burn-in") + lab = "warmup") mcmc_l %>% dplyr::select(iter, beta_0, beta_1) %>% pivot_longer(cols = c("beta_0", "beta_1"), names_to = "param") %>% @@ -519,11 +525,15 @@ Some alternatives also exist: - [**rstan**](https://cran.r-project.org/web/packages/rstan/vignettes/rstan.html) is the basis of **brms**. You need **rstan** to use **brms** because it communicates directly with Stan. A model composed in "pure" Stan code can be fitted with the **rstan** package in R (so without using **brms**), but can be quite complex. - [**rstanarm**](https://mc-stan.org/rstanarm/articles/index.html) is similar to **brms** in terms of complexity of the syntax. Moreover, with the **rstanarm** package you can use some pre-compiled Stan models. This means that fitting your model takes less time. A large part of the process time of **brms** goes into compiling the Stan program. -![Overview of various Stan software (source: https://jtimonen.github.io/posts/post-01/)](software.png) +![Overview of various Stan software (source: https://jtimonen.github.io/posts/post-01/)](./stan_software.png) + + +Although the package **brms** depends on **rstan** as a default backend, there is another software tool worth mentioning: **cmdstanr**. +[Compared to **rstan**, **cmdstanr**]( https://mc-stan.org/cmdstanr/articles/cmdstanr.html#comparison-with-rstan) benefits from more frequent updates and some performance advantages, yet it requires [extra installation steps](https://mc-stan.org/cmdstanr/articles/cmdstanr.html#installing-cmdstan). +In the `brm(...)` function introduced below, you can switch these two with the `backend = ...` keyword. -# Fitting a model with brms -## Loading the dataset and data exploration +# Loading the dataset and data exploration We load a dataset on the number of ant species in New England (USA). Type `?ants` into the console for more info. @@ -593,6 +603,9 @@ ants_df %>% As an exercise, we will create a model to compare the number of species between both habitats. From the data exploration, we already saw that the number of species seems to be higher in forests and that sites with a higher number in bogs often also have a higher number in forests. + +# Fitting a model with `brms` + ## Specification of a linear regression ### Model specification @@ -617,8 +630,8 @@ First of all we decide which MCMC parameters we will use. Type `?brm` to see wha ```{r simple-model-mcmc-par} # Set MCMC parameters nchains <- 3 # number of chains -niter <- 2000 # number of iterations (incl. burn-in, see next) -burnin <- niter / 4 # number of initial samples to remove (= burn-in) +niter <- 2000 # number of iterations (incl. warmup, see next) +warmup <- niter / 4 # number of initial samples to remove (= warmup) nparallel <- nchains # number of cores for parallel computing thinning <- 1 # thinning factor (here 1 = no thinning) ``` @@ -630,14 +643,14 @@ The model is fitted using the `brm()` function. The syntax is very similar to fu - `file` and `file_refit` to save the model object after it has been fitted. If you run the code again and the model has already been saved, `brm()` will simply load this model instead of refitting it. -```{r simple-model-fit-poisson} +```{r simple-model-fit-poisson, class.source='fold-show'} # Fit Normal model fit_normal1 <- brm( formula = sp_rich ~ habitat, # specify the model family = gaussian(), # we use the Normal distribution data = ants_df, # specify data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -649,7 +662,7 @@ Before we look at the results, we first check whether the model converges well. ### MCMC convergence -There are several ways to check convergence of the MCMC algorithm for each parameter. The burn-in samples are not taken into account. First and foremost, you have *visual controls*. +There are several ways to check convergence of the MCMC algorithm for each parameter. The warmup samples are not taken into account. First and foremost, you have *visual controls*. We can obtain the MCMC samples with the `as_draws()` functions or visualise them at once via the [**bayesplot**](https://mc-stan.org/bayesplot/) package that is compatible with brmsfit objects. @@ -714,7 +727,7 @@ as_draws_df(fit_normal1, variable = parameters) %>% pivot_longer(cols = starts_with(parameters), names_to = "name", values_to = "value") %>% # split columns names at the last underscore - extract(name, into = c("parameter", "statistic"), "(.*)_([^_]+)$") %>% + tidyr::extract(name, into = c("parameter", "statistic"), "(.*)_([^_]+)$") %>% ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) + geom_line(aes(linetype = statistic), linewidth = 0.8) + facet_wrap(~parameter, nrow = 3, scales = "free") @@ -836,14 +849,14 @@ $$ So we need to estimate two parameters: $\beta_0$ and $\beta_1$ We use the same MCMC parameters as before. The only thing we need to adjust is the choice `family = poisson()`. -```{r poisson-model-fit} +```{r poisson-model-fit, class.source='fold-show'} # Fit Poisson model fit_poisson1 <- brm( formula = sp_rich ~ habitat, # specify the model family = poisson(), # we use the Poisson distribution data = ants_df, # specify the data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -898,14 +911,14 @@ $$ b_0 \sim N(0, \sigma_b) $$ -```{r rand-intercept-model-fit} +```{r rand-intercept-model-fit, class.source='fold-show'} # Fit Poisson model with random intercept per site fit_poisson2 <- brm( formula = sp_rich ~ habitat + (1|site), family = poisson(), data = ants_df, chains = nchains, - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -937,15 +950,15 @@ pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, How can we objectively compare these models? -# Compare models +## Compare models Based on the PPCs we can already see which model fits the data best. Furthermore, there are some functions that **brms** provides to compare different models. With the function `add_criterion()` you can add model fit criteria to model objects. Type `?add_criterion()` to see which ones are available. See also -## Leave-one-out cross validation +### Leave-one-out cross validation Cross-validation (CV) is a family of techniques that attempts to estimate how well a model would predict unknown data through predictions of the model fitted to the known data. You do not necessarily have to collect new data for this. You can split your own data into a test and training dataset. You fit the model to the training dataset and then use that model to estimate how well it can predict the data in the test dataset. With leave-one-out CV (LOOCV) you leave out one observation each time as test dataset and refit the model based on all other observations (= training dataset). -![Figure LOOCV.](loocv.png) +![Figure LOOCV.](./loocv.png) ```{r compare-loocv} # Add leave-one-out model fit criterium to model objects @@ -977,11 +990,11 @@ comp_loo %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## K-fold cross-validation +### K-fold cross-validation With K-fold cross-validation, the data is split into $K$ groups. We will use $K = 10$ groups (= folds) here. So instead of leaving out a single observation each time, as with leave-one-out CV, we will leave out one $10^{th}$ of the data here. Via the arguments `folds = "stratified"` and `group = "habitat"` we ensure that the relative frequencies of habitat are preserved for each group. This technique will therefore be less precise than the previous one, but will be faster to calculate if you work with a lot of data. -![Figure K-fold CV.](k-foldcv.png) +![Figure K-fold CV.](./k-foldcv.png) ```{r compare-K-fold} # Add K-fold model fit criterium to model objects @@ -1023,7 +1036,7 @@ comp_kfold %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## WAIC +### WAIC The Widely Applicable Information Criterion (WAIC) does not use cross-validation but is a computational way to estimate the ELPD. How this happens exactly is beyond the purpose of this tutorial. It is yet another measure to apply model selection. @@ -1055,10 +1068,12 @@ comp_waic %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## Conclusion +### Conclusion Both based on the PPC and the comparisons with different model selection criteria, we can conclude that the second Poisson model with random intercepts fits the data best. In principle, we could have expected this based on our own intuition and the design of the study, i.e. the use of the Poisson distribution to model numbers and the use of random intercepts to control for a hierarchical design (habitats nested within sites). + + # Final model results When we look at the model fit object, we see results that are similar to results we see when we fit a frequentist model. On the one hand we get an estimate of all parameters with their uncertainty, but on the other hand we see that this is clearly the output of a Bayesian model. We get information about the parameters we used for the MCMC algorithm, we get a 95% credible interval (CI) instead of a confidence interval and we also get the $\hat{R}$ value for each parameter as discussed earlier. @@ -1185,6 +1200,181 @@ t_test_normal1 We see that this indeed produces almost exactly the same results. Our Bayesian model estimates that on average `r round(diff_bog1, 3)` more ant species occur in forests than in bogs (90% credible interval: `r round(ll_diff_bog1, 3)` to `r round(ul_diff_bog1, 3)`). The t-test estimates that on average `r round(diff_bog2, 3)` more ant species occur in forests than in bogs (90% confidence interval: `r round(ll_diff_bog2, 3)` to `r round(ul_diff_bog2, 3)`). +# Deep Dive: `rstan` +## Stan: What? Why?! +The `brms` package is a convenience wrapper for the `rstan` package, which in turn ports `stan` functionality to R. +Stan is a modeling framework written in the `C` programming language, which implements many probabilistic ("Bayesian") modeling tools. +More info can be found on [the Stan website](https://mc-stan.org). + + +The advantage of `brms` is usability: many functions work out-of-the-box, with reasonable default values, and a syntax that is similar to what many statisticians are habituated to. +However, the relative ease-of-use comes at the cost of flexibility, and do some degree, readability. + +In contrast, Stan and `rstan` lean more to the mathematical formulation of models. +Every aspect of the model has to be explicitly set, which can be an advantage (e.g. if you face non-standard use cases), or disadvantage (e.g. if you specify models in non-optimal ways). + + +To briefly give an impression, we will build the same models as above, using the Stan framework. + + +```{r} +library("rstan") +conflicted::conflicts_prefer(rstan::extract) +conflicted::conflicts_prefer(brms::loo) +``` + +## Model Definition +RMarkdown can handle `stan` code chunks, though more general model definition is outsourced to a separate "*.stan" file (e.g. [here](https://github.com/inbo/tutorials/tree/master/static/code/brms_modeling/poisson_model.stan)). +Alternatively, you can define your model in a big text block, as shown below. +The simple poisson model resembles [one of the `stan`-dard examples](https://mc-stan.org/docs/stan-users-guide/posterior-prediction.html#posterior-prediction-for-regressions), which you can refer to for all further details and more. + + +```{r stan_poisson_simple, eval=TRUE, output.var="", class.source='fold-show'} +stan_poisson_model_code = ' +data { + int N; // sample size + vector[N] habitat_obs; // habitat + array[N] int sp_rich; // outcome variable: number of ants +} + +parameters { + real intercept_rich; // intercept + real habitat_effect; // slope +} + +model { + // priors + intercept_rich ~ normal(0, 1); + habitat_effect ~ normal(0, 1); + + // "posterior", "likelihood", ... give it a name! + sp_rich ~ poisson_log(intercept_rich + habitat_effect * habitat_obs); + +} + +' +``` + + +Take this as a "look behind the scenes"! +You have to explicitly define the model structure, priors, even the data types of input variables. +Yet note how even Stan is not without convenience functions: we can use the `poisson_log` posterior to model log rates. + + +When working outside RStudio/RMarkdown, you might prefer loading the model from a file: + +```{r stan_load_model, eval=TRUE, class.source='fold-show'} +stan_poisson_model <- stan_model( + # file = here("code", "brms_modeling", "poisson_model.stan"), + model_code = stan_poisson_model_code, + model_name = "stan poisson model" + ) +``` + + +## Sampling + +Sampling does pretty much the same as above, since at the core, `brms` is just `stan`. + +```{r stan_sampling} +stan_poisson_fit <- sampling( + stan_poisson_model, + list( + N = nrow(ants_df), + habitat_obs = as.integer(ants_df$habitat)-1, + sp_rich = ants_df$sp_rich + ), + iter = niter, + chains = nchains, + cores = nparallel + ) + +stan_poisson_fit +``` + + +A convenient way of quickly inspecting model fits is `shinystan`. +It opens a browser with shiny plots. + +```{r eval=FALSE} +library("shinystan") +launch_shinystan(stan_poisson_fit) +``` + + +Behold: model outcome is exactly as above. +In the present case, the benefit from turning to `stan` is very limited. +In other cases, it might pay off. + +Know that Stan is there for you, do not hesitate to turn to its extensive documentation, and do not fear to give it a try! + + +## Homework: Hierarchical Model +To take your modeling skills even further, you may implement and sample the "random intercept" model. +In "Bayesian" terms, the [general terminology is "hierarchical" model](https://mc-stan.org/docs/stan-users-guide/regression.html#hierarchical-regression). + + +Below is the code which considers site-specific intercepts. +It works slightly different from the `+(1|site)` approach above: the upper one is an "offset" parametrization, whereas this time we use a hyperprior. +These are two common strategies often worth interchanging, with marginal or substantial effects on sampler performance. + + +```{stan, eval=FALSE, output.var="stan_hierarchical_poisson"} +data { + int N; // sample size + int L; // number of sites + vector[N] habitat_obs; // habitat + array[N] int site_obs; // site + array[N] int sp_rich; // outcome variable: number of ants +} + +parameters { + real intercept_rich; // "global" intercept + vector[L] site_effect; // the sr intercept at each site + real site_variation; // variation on site level + real habitat_effect; // habitat slope +} + +model { + // priors + intercept_rich ~ normal(0, 1); + site_variation ~ cauchy(0, 1); + habitat_effect ~ normal(0, 1); + + // site-wise intercept + for (i in 1:L) { + site_effect[i] ~ normal(intercept_rich, site_variation); + } + + // "posterior", "likelihood", ... give it a name! + sp_rich ~ poisson_log(site_effect[site_obs] + habitat_effect * habitat_obs); + +} +``` + + +```{r, eval=FALSE, stan_hierarchical_sampling} +stan_poisson_fit <- sampling( + stan_poisson_model, + list( + N = nrow(ants_df), + L = length(levels(ants_df$site)), + habitat_obs = as.integer(ants_df$habitat)-1, + site_obs = as.integer(ants_df$site), + sp_rich = ants_df$sp_rich + ), + iter = niter, + chains = nchains, + cores = nparallel + ) + +stan_poisson_fit +``` + +With Stan, po(i)ssibilities are almost endless - don't get lost in model building! + + + # References {-} - Stan User Guide [(link)](https://mc-stan.org/docs/2_18/stan-users-guide/index.html). diff --git a/content/tutorials/r_brms/brms_nl/bayesian_statistics_1_script.R b/content/tutorials/r_brms/brms_nl/bayesian_statistics_1_script.R index 14f84b77b..717271bf7 100644 --- a/content/tutorials/r_brms/brms_nl/bayesian_statistics_1_script.R +++ b/content/tutorials/r_brms/brms_nl/bayesian_statistics_1_script.R @@ -48,7 +48,7 @@ confint(lm1, level = 0.9) ## ------------------------------------------------------------------------------------------------------------------------------------- -source(file = "./source/mcmc_functions.R") +source(file = here::here("code", "brms_modeling", "mcmc_functions.R")) ## ------------------------------------------------------------------------------------------------------------------------------------- @@ -232,7 +232,7 @@ fit_normal1 <- brm( family = gaussian(), # we gebruiken de Normaal verdeling data = ants_df, # ingeven data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -298,7 +298,7 @@ as_draws_df(fit_normal1, variable = parameters) %>% ## ----simpel-model-posterior-density2-------------------------------------------------------------------------------------------------- -# Visualisatie density plot van de posterior voor ieder van de chains apart in +# Visualisatie density plot van de posterior voor ieder van de chains apart in # overlay. via Bayesplot package mcmc_dens_overlay(fit_normal1, pars = parameters) @@ -337,7 +337,7 @@ mcmc_neff(ratios_fit_normal1) + yaxis_text(hjust = 1) ## ----simpel-model-fit1---------------------------------------------------------------------------------------------------------------- # Visualiseer model fit via Bayesplot package -pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -348,7 +348,7 @@ fit_poisson1 <- brm( family = poisson(), # we gebruiken de Poisson verdeling data = ants_df, # ingeven data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -368,7 +368,7 @@ mcmc_rhat(rhats_fit_poisson1) + yaxis_text(hjust = 1) ## ----poisson-model-fit-vis------------------------------------------------------------------------------------------------------------ # Visualiseer model fit via Bayesplot package -pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -379,7 +379,7 @@ fit_poisson2 <- brm( family = poisson(), data = ants_df, chains = nchains, - warmup = burnin, + warmup = burnin, iter = niter, cores = nparallel, thin = thinning, @@ -401,7 +401,7 @@ mcmc_rhat(rhats_fit_poisson2) + yaxis_text(hjust = 1) ## ----rand-intercept-model-fit-vis----------------------------------------------------------------------------------------------------- # Visualiseer model fit van het Poisson model met random intercept via # Bayesplot package -pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, +pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, group = "habitat") @@ -520,7 +520,7 @@ fit_poisson2 %>% # bereken gemiddelde aantallen en zet om naar lang formaat voor visualisatie mutate(bog = exp(b_Intercept), forest = exp(b_Intercept + b_habitatForest)) %>% - pivot_longer(cols = c("bog", "forest"), names_to = "habitat", + pivot_longer(cols = c("bog", "forest"), names_to = "habitat", values_to = "sp_rich") %>% # visualiseer via ggplot() ggplot(aes(y = sp_rich, x = habitat)) + diff --git a/content/tutorials/r_brms/brms_nl/install_packages.R b/content/tutorials/r_brms/brms_nl/install_packages.R deleted file mode 100644 index 7b8139503..000000000 --- a/content/tutorials/r_brms/brms_nl/install_packages.R +++ /dev/null @@ -1,7 +0,0 @@ -#https://learnb4ss.github.io/learnB4SS/articles/install-brms.html -#first install rstan (more detailed instructions here https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started ) -install.packages("rstan", repos = "https://cloud.r-project.org/", dependencies = TRUE) -#verify the installation with: -example(stan_model, package = "rstan", run.dontrun = TRUE) #this might take a while the first time -#if the example works, install brms: -install.packages("brms") diff --git a/content/tutorials/r_brms/brms_nl/mcmc_functions.R b/content/tutorials/r_brms/brms_nl/mcmc_functions.R deleted file mode 100644 index 27d509ac6..000000000 --- a/content/tutorials/r_brms/brms_nl/mcmc_functions.R +++ /dev/null @@ -1,93 +0,0 @@ -# Likelihood function for the data (log) -ll.fun <- function(beta_0, beta_1, tau = 0.2){ - ll <- sum(dnorm(x = y, mean = beta_0 + beta_1 * x, - sd = sqrt(1 / tau), log = TRUE)) - return(ll) -} - -# prior function (uninformative) (log) -prior.fun <- function(beta_0, beta_1, tau = 0.2){ - lp <- dnorm(x = beta_0, mean = 0, sd = sqrt(1/10^-6), log = TRUE) + - dnorm(x = beta_1, mean = 0, sd = sqrt(1/10^-6), log = TRUE) + - dgamma(x = tau, shape = 0.001, rate = 0.001, log = TRUE) - return(lp) -} - -# MCMC - metropolis algorithm - -# N = lengte van de MCMC -# sigma = grootte van de stappen in de markov chain -# init_b0 en init_b1 = initiële waarden voor beta_0 en beta_1 - -MCMC_metro <- function(x, y, N = 20, sigma = 1, init_b0, init_b1){ - - set.seed(56) # Optioneel, om bij iedereen hetzelfde resultaat te bekomen. - # We maken een lege dataframe om de output in te bewaren - mcmc <- data.frame(iter = 1:N, # positie in de markov chain - beta_0 = NA, # intercept - beta_1 = NA, # helling - tau = NA, # 1/standard deviation - post = NA, # posterior - accept = NA) # stap wel (1) of niet (0) geaccepteerd. - - mcmc$beta_0[1] <- init_b0 - mcmc$beta_1[1] <- init_b1 - mcmc$tau[1] <- 0.2 # init_tau - - # Berekening van de posterior voor de startwaarde - mcmc$post[1] <- ll.fun(beta_0 = mcmc$beta_0[1], - beta_1 = mcmc$beta_1[1], - tau = mcmc$tau[1]) + - prior.fun(beta_0 = mcmc$beta_0[1], - beta_1 = mcmc$beta_1[1], - tau = mcmc$tau[1]) - mcmc$accept[1] <- 1 # De initiële parameterwaarde altijd accepteren - - for(i in 2:N){ - - # Kies een nieuwe waarde uit een normale verdeling met - # mean = vorige parameterwaarde en sd = sigma (stapgrootte) - new_beta_0 <- rnorm(1, mean = mcmc$beta_0[i - 1], sd = sigma) - new_beta_1 <- rnorm(1, mean = mcmc$beta_1[i - 1], sd = sigma) - #new_tau <- rnorm(1, mean = mcmc$tau[i - 1], sd = 0.1) - new_tau <- mcmc$tau[i -1] #we houden tau fixed voor de eenvoud - - # Bereken de posterior voor de nieuwe parameterwaarden - new_post <- ll.fun(beta_0 = new_beta_0, - beta_1 = new_beta_1, - tau = new_tau) + - prior.fun(beta_0 = new_beta_0, - beta_1 = new_beta_1, - tau = new_tau) - - # Bereken de ratio van de posterios (in de log schaal is dit het verschil) - log.ratio <- new_post - mcmc$post[i-1] - - # Als de ratio > 0 => nieuwe posterior beter = > accepteren - if (log.ratio > 0) { - mcmc$accept[i] <- 1 - # Als log.ratio <= 0 is de kans op accepteren afhankelijk van de log.ratio - }else{ - if (rbinom(n = 1, size = 1, prob = exp(log.ratio))) { - mcmc$accept[i] <- 1 - }else{ - mcmc$accept[i] <- 0 - } - } - - # De nieuwe waarde wordt geaccepteerd - if(mcmc$accept[i]) { - mcmc$beta_0[i] <- new_beta_0 - mcmc$beta_1[i] <- new_beta_1 - mcmc$tau[i] <- new_tau - mcmc$post[i] <- new_post - }else{ - # De nieuwe waarde wordt niet geaccepteerd (de oude behouden) - mcmc$beta_0[i] <- mcmc$beta_0[i-1] - mcmc$beta_1[i] <- mcmc$beta_1[i-1] - mcmc$tau[i] <- mcmc$tau[i-1] - mcmc$post[i] <- mcmc$post[i-1] - } - } - return(mcmc = mcmc) -} diff --git a/content/tutorials/r_brms/brms_nl/workshop_1_mcmc_en_brms.Rmd b/content/tutorials/r_brms/brms_nl/workshop_1_mcmc_en_brms.Rmd index 05f228381..9591d8bdb 100644 --- a/content/tutorials/r_brms/brms_nl/workshop_1_mcmc_en_brms.Rmd +++ b/content/tutorials/r_brms/brms_nl/workshop_1_mcmc_en_brms.Rmd @@ -35,12 +35,18 @@ library(tidybayes) # nabewerking en visualisatie van Bayesiaanse modellen # Conflicten tussen packages conflicted::conflicts_prefer(dplyr::filter) conflicted::conflicts_prefer(dplyr::lag) +conflicted::conflicts_prefer(tidyr::extract) conflicted::conflicts_prefer(brms::ar) conflicted::conflicts_prefer(brms::dstudent_t) conflicted::conflicts_prefer(brms::pstudent_t) conflicted::conflicts_prefer(brms::qstudent_t) conflicted::conflicts_prefer(brms::rstudent_t) conflicted::conflict_prefer("rhat", "brms") + +# html bestandsgrootte verminderen door grafiekparameters te specifiëren +knitr::opts_chunk$set(dev = "jpeg", dpi = 56, + fig.width = 4, fig.height = 3, + out.width = "60%") ``` # Theoretische achtergrond @@ -234,7 +240,7 @@ Bij een stap bergop gaat de MCMC altijd door. Is de stap bergaf, dan gaat de MCM De absolute hoogte van de berg (posterior likelihood) doet er op zich niet toe. Enkel de relatieve verschillen. -![Metropolis regels](../media/MCMC4b.png) +![Metropolis regels](./mcmc4b.png) Een eenvoudig MCMC-algoritme is het Metropolis algoritme. Het wordt beschreven door deze formele regels: @@ -293,7 +299,8 @@ confint(lm1, level = 0.9) We lezen enkele korte functies in om de (log) likelihood en de prior te berekenen en het MCMC metropolis algoritme uit te voeren. ```{r} -source(file = "./source/mcmc_functions.R") + +source(file = here("static", "code", "brms_modeling", "mcmc_functions.R")) ``` Voor dit eenvoudig model met een kleine dataset kunnen we de posterior voor een groot aantal combinaties voor `beta_0` en `beta_1` uitrekenen en plotten. @@ -340,7 +347,7 @@ ggplot(df, aes(x = beta_0, y = beta_1, z = post)) + theme(legend.position="none") ``` -De startwaarde ligt een eind van de maximale likelihood. Het duurt even voor de MCMC stabiliseert in de omgeving van de top van de berg. Daarom wordt het eerste deel van de MCMC nooit gebruikt. Dit is de 'burn-in' (of warmp-up in brms). +De startwaarde ligt een eind van de maximale likelihood. Het duurt even voor de MCMC stabiliseert in de omgeving van de top van de berg. Daarom wordt het eerste deel van de MCMC nooit gebruikt. Dit is de 'warmup' (engels voor opwarmfase; ook 'burn-in' of 'tuning' in andere bibliotheken). We runnen nu hetzelde model, maar met een veel langere mcmc @@ -387,7 +394,7 @@ Voor elke parameter kan de MCMC worden weergegeven in een 'trace plot'. ann_text <- tibble(param = c("beta_0", "beta_1"), x = 210, y = c(max(mcmc_l$beta_0), max(mcmc_l$beta_1)), - lab = "burn-in") + lab = "warmup") mcmc_l %>% dplyr::select(iter, beta_0, beta_1) %>% pivot_longer(cols = c("beta_0", "beta_1"), names_to = "param") %>% @@ -506,12 +513,14 @@ Er bestaan ook enkele alternatieven: - [RStan](https://cran.r-project.org/web/packages/rstan/vignettes/rstan.html) is de basis van `brms`. Je hebt `RStan`nodig om `brms`te kunnen gebruiken omdat het rechtstreeks met Stan communiceert. Een model definiëren in "Pure" Stan code is echter redelijk omslachtig. - [rstanarm](https://mc-stan.org/rstanarm/articles/index.html) is nog makkelijker dan `brms`qua syntax. Bovendien kan je met `rstanarm` gebruik maken van enkele vooral gecompileerde Stan modellen. Hierdoor duurt het fitten van je model minder lang. Een groot deel van de procestijd van `brms`gaat immers naar het compileren van het Stan programma. -![Overzicht van verschillende Stan software (bron: https://jtimonen.github.io/posts/post-01/)](software.png) +![Overzicht van verschillende Stan software (bron: https://jtimonen.github.io/posts/post-01/)](./stan_software.png) +Hoewel het pakket **brms** afhankelijk is van **rstan** als standaard backend, is er nog een andere softwaretool die het vermelden waard is: **cmdstanr**. +[Vergeleken met **rstan** profiteert **cmdstanr**]( https://mc-stan.org/cmdstanr/articles/cmdstanr.html#comparison-with-rstan) van frequentere updates en enkele performance-voordelen, maar vereist het [extra installatiestappen](https://mc-stan.org/cmdstanr/articles/cmdstanr.html#installing-cmdstan). +In de hieronder geïntroduceerde functie `brm(...)` kunt u deze twee omwisselen met het trefwoord `backend = ...`. -# Een model fitten met brms -## Dataset laden en data exploratie +# Dataset laden en data exploratie We laden een dataset in over het aantal mierensoorten in New England (USA). Typ `?ants` in de console voor meer info. @@ -583,6 +592,9 @@ ants_df %>% Als oefening zullen we een model maken om het aantal soorten te vergelijken tussen beide habitats. Uit de data exploratie zagen we al dat het aantal hoger lijkt te liggen in bossen en dat sites met een hoger aantal in moerassen vaak ook een hoger aantal in bossen hebben. + +# Een model fitten met `brms` + ## Specificatie van een lineaire regressie ### Model specificatie @@ -607,8 +619,8 @@ Eerst en vooral besluiten we welke MCMC parameters we zullen gebruiken. Typ `?br ```{r simpel-model-mcmc-par} # Instellen MCMC parameters nchains <- 3 # aantal chains -niter <- 2000 # aantal iteraties (incl. burn-in, zie volgende) -burnin <- niter / 4 # aantal initiële samples om te verwijderen (= burn-in) +niter <- 2000 # aantal iteraties (incl. warmup, zie volgende) +warmup <- niter / 4 # aantal initiële samples om te verwijderen (= warmup) nparallel <- nchains # aantal cores voor parallel computing thinning <- 1 # verdunningsfactor (hier 1 = geen verdunning) ``` @@ -620,14 +632,14 @@ Het model wordt gefit a.d.h.v. de `brm()` functie. De syntax is zeer gelijkaardi - `file` en `file_refit` om het model object op te slaan nadat het gefit is. Als je de code opnieuw runt en het model is al eens opgeslaan, dan zal `brm()` dit model gewoon inladen in plaats van het opnieuw te fitten. -```{r simpel-model-fit-poisson} +```{r simpel-model-fit-poisson, class.source = 'fold-show'} # Fit Normaal model fit_normal1 <- brm( formula = sp_rich ~ habitat, # beschrijving van het model family = gaussian(), # we gebruiken de Normaal verdeling data = ants_df, # ingeven data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -638,7 +650,7 @@ Voor we de resultaten bekijken, controleren we eerst of het model goed convergee ### MCMC convergentie -Er zijn verschillende manieren om convergentie van het MCMC algoritme voor elke parameter te controleren. Hierbij worden de burn-in samples niet in rekening genomen. Eerst en vooral heb je *visuele controles*. +Er zijn verschillende manieren om convergentie van het MCMC algoritme voor elke parameter te controleren. Hierbij worden de warmup samples niet in rekening genomen. Eerst en vooral heb je *visuele controles*. We kunnen de MCMC samples met de `as_draws()` functies verkrijgen ofwel ineens visualiseren via de [bayesplot](https://mc-stan.org/bayesplot/) package die compatibel is met brmsfit objecten. @@ -703,7 +715,7 @@ as_draws_df(fit_normal1, variable = parameters) %>% pivot_longer(cols = starts_with(parameters), names_to = "name", values_to = "value") %>% # splits naam kolom op bij laatste underscore - extract(name, into = c("parameter", "statistiek"), "(.*)_([^_]+)$") %>% + tidyr::extract(name, into = c("parameter", "statistiek"), "(.*)_([^_]+)$") %>% ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) + geom_line(aes(linetype = statistiek), linewidth = 0.8) + facet_wrap(~parameter, nrow = 3, scales = "free") @@ -825,14 +837,14 @@ $$ We moeten dus twee parameters schatten: $\beta_0$ en $\beta_1$ We gebruiken dezelfde MCMC parameters als voordien. Het enige wat we moeten aanpassen is de keuze `family = poisson()`. -```{r poisson-model-fit} +```{r poisson-model-fit, class.source = 'fold-show'} # Fit Poisson model fit_poisson1 <- brm( formula = sp_rich ~ habitat, # beschrijving van het model family = poisson(), # we gebruiken de Poisson verdeling data = ants_df, # ingeven data chains = nchains, # MCMC parameters - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -888,14 +900,14 @@ $$ b_0 \sim N(0, \sigma_b) $$ -```{r rand-intercept-model-fit} +```{r rand-intercept-model-fit, class.source = 'fold-show'} # Fit Poisson model met random intercepten per site fit_poisson2 <- brm( formula = sp_rich ~ habitat + (1|site), family = poisson(), data = ants_df, chains = nchains, - warmup = burnin, + warmup = warmup, iter = niter, cores = nparallel, thin = thinning, @@ -929,15 +941,15 @@ Hoe kunnen we deze modellen nu objectief gaan vergelijken? -# Vergelijken van modellen +## Vergelijken van modellen Op basis van de PPCs kunnen we reeds zien welk model het best past bij de data. Verder zijn er nog enkele functies die **brms** voorziet om verschillende modellen te vergelijken. Met de functie `add_criterion()` kan je model fit criteria toevoegen aan model objecten. Typ `?add_criterion()` om te zien welke beschikbaar zijn. Zie ook -## Leave-one-out cross-validation +### Leave-one-out cross-validation Cross-validation (CV) is een familie van technieken die probeert in te schatten hoe goed een model onbekende data zou voorspellen via predicties van het model gefit op de bekende data. Hiervoor moet je niet per se nieuwe data gaan inzamelen. Je kan jouw eigen data opsplitsen in een test en training dataset. Je fit het model op de training dataset en je gebruikt dat model dan om te schatten hoe goed het de data in de test dataset kan voorspellen. Bij leave-one-out CV (LOOCV) ga je telkens één observatie weglaten en het model opnieuw fitten o.b.v. alle andere observaties. -![Illustratie LOOCV.](loocv.png) +![Illustratie LOOCV.](./loocv.png) ```{r vergelijken-loocv} # Voeg leave-one-out model fit criterium toe aan model objecten @@ -969,11 +981,11 @@ comp_loo %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## K-fold cross-validation +### K-fold cross-validation Bij K-fold cross-validation worden de data in $K$ groepen opgesplitst. Wij zullen hier $K = 10$ groepen (= folds) gebruiken. In plaats van dus telkens één enkele observatie weg te laten zoals bij leave-one-out CV gaan we hier $1/10$e van de data weglaten. Via de argumenten `folds = "stratified"` en `group = "habitat"` zorgen we ervoor dat voor elke groep de relatieve frequenties van habitat bewaard blijven. Deze techniek zal dus minder precies zijn dan de vorige, maar zal sneller zijn om te berekenen indien je met heel veel data werkt. -![Illustratie K-fold CV.](k-foldcv.png) +![Illustratie K-fold CV.](./k-foldcv.png) ```{r vergelijken-K-fold} # Voeg K-fold model fit criterium toe aan model objecten @@ -1014,7 +1026,7 @@ comp_kfold %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## WAIC +### WAIC Het Widely Applicable Information Criterion (WAIC) maakt geen gebruik van cross-validation maar is een computationele manier om de ELPD te schatten. Hoe dit precies gebeurt valt buiten het doel van deze workshop. Het is nog een andere maat om model selectie toe te passen. @@ -1046,10 +1058,12 @@ comp_waic %>% ul_diff = elpd_diff + qnorm(0.95) * se_diff) ``` -## Conclusie +### Conclusie Zowel o.b.v. de PPC als vergelijkingen met verschillende model selectie criteria, kunnen we besluiten dat het tweede Poisson model met random intercepts het best past bij de data. In principe konden we dit ook verwachten op basis van onze eigen intuïtie en het design van de studie, nl. het gebruik van de Poisson distributie om aantallen te modelleren en het gebruik van random intercepts om te controleren voor een hiërarchisch design (habitats genest in sites). + + # Resultaten finale model Als we het model fit object bekijken, zien we resultaten die vergelijkbaar zijn met resultaten zoals we die zien als we een frequentist model gefit hebben. We krijgen enerzijds een schatting van alle parameters met hun onzekerheid maar anderzijds zien we dat dit duidelijk de output is van een Bayesiaans model. Zo krijgen we info over de parameters die we gebruikt hebben voor het MCMC algoritme, we krijgen een 95 % credible interval (CI) in plaats van een confidence interval en we krijgen bij elke parameter ook de $\hat{R}$ waarde die we eerder besproken hebben. @@ -1174,6 +1188,181 @@ t_test_normal1 We zien dat dit inderdaad bijna exact dezelfde resultaten oplevert. Ons Bayesiaans model schat dat er gemiddeld `r round(verschil_moeras1, 3)` meer mierensoorten in bossen voorkomen dan in moerassen (90 % credible interval: `r round(ll_verschil_moeras1, 3)` tot `r round(ul_verschil_moeras1, 3)`). De t-test schat dat er gemiddeld `r round(verschil_moeras2, 3)` meer mierensoorten in bossen voorkomen dan in moerassen (90 % confidence interval: `r round(ll_verschil_moeras2, 3)` tot `r round(ul_verschil_moeras2, 3)`). + +# Deep Dive: `rstan` +## Stan: Wat? Waarom?! +De `brms` package is een gemaksverpakking van `rstan`, dat zelf de functionaliteit van `stan` naar R brengt. +Stan is een modeleringsomgeving, geschreven in de programmeertaal `C`, die de benodigdheden voor probabilistische ("Bayesiaanse") modellen implementeerd. + +Meer info vind je [op het internet](https://mc-stan.org). + +Het voordeel van `brms` is bruikbaarheid: veel functies werken "out-of-the-box", met redelijke standaardwaarden en een syntaxis die lijkt op wat veel statistici gewend zijn. +Het relatieve gebruiksgemak kan echter ten koste gaan van de flexibiliteit en in zekere mate van de leesbaarheid. + +Stan en `rstan` daarentegen leunen meer op de wiskundige formulering van modellen. +Elk aspect van het model moet expliciet worden ingesteld, wat een voordeel kan zijn (bijvoorbeeld als u te maken hebt met niet-standaard usecases) of een nadeel (bijvoorbeeld als u modellen op niet-optimale manieren specificeert). + +Om kort een indruk te geven, bouwen we dezelfde modellen als hierboven, met behulp van het Stan-framework. + + +```{r} +library("rstan") +conflicted::conflicts_prefer(rstan::extract) +conflicted::conflicts_prefer(brms::loo) +``` + +## Model Definition +RMarkdown kan `stan`-chunks verwerken, hoewel de algemenere modeldefinitie is uitbesteed aan een apart "*.stan"-bestand (bv. [hier](https://github.com/inbo/tutorials/tree/master/static/code/brms_modeling/poisson_model.stan)). +Alternatief kan je je model ook als een groot tekstblok definiëren, zoals hieronder te zien. +Het eenvoudige poissonmodel lijkt op [een van de `stan`-daardvoorbeelden](https://mc-stan.org/docs/stan-users-guide/posterior-prediction.html#posterior-prediction-for-regressions), waarnaar u terecht kunt voor verdere details en meer. + + +```{r stan_poisson_simple, eval=TRUE, output.var="", class.source='fold-show'} +stan_poisson_model_code = ' +data { + int N; // sample size + vector[N] habitat_obs; // habitat + array[N] int sp_rich; // outcome variable: number of ants +} + +parameters { + real intercept_rich; // intercept + real habitat_effect; // slope +} + +model { + // priors + intercept_rich ~ normal(0, 1); + habitat_effect ~ normal(0, 1); + + // "posterior", "likelihood", ... give it a name! + sp_rich ~ poisson_log(intercept_rich + habitat_effect * habitat_obs); + +} + +' +``` + + + +Zie dit als een "kijkje achter de schermen"! +Je moet expliciet de modelstructuur, priors en zelfs de gegevenstypen van invoervariabelen definiëren. +Maar let op hoe zelfs Stan niet zonder handige functies is: we kunnen de `poisson_log` posterior gebruiken om de $log(\lambda )$ te modelleren. + +Als je buiten RStudio/RMarkdown werkt, kun je het model het beste laden vanuit een bestand: + +```{r stan_load_model, eval=TRUE, class.source='fold-show'} +stan_poisson_model <- stan_model( + # file = here("code", "brms_modeling", "poisson_model.stan"), + model_code = stan_poisson_model_code, + model_name = "stan poisson model" + ) +``` + + +## Sampling + +Sampling doet vrijwel hetzelfde als hierboven, want in de kern is `brms` gewoon `stan`. + +```{r stan_sampling} +stan_poisson_fit <- sampling( + stan_poisson_model, + list( + N = nrow(ants_df), + habitat_obs = as.integer(ants_df$habitat)-1, + sp_rich = ants_df$sp_rich + ), + iter = niter, + chains = nchains, + cores = nparallel + ) + +stan_poisson_fit +``` + +Een handige manier om modelfits snel te inspecteren is `shinystan`. +Het opent een browser met glanzende plots. + +```{r eval=FALSE} +library("shinystan") +launch_shinystan(stan_poisson_fit) +``` + + +En kijk: de modeluitkomst is precies zoals boven. +In het huidige geval is het voordeel van het overstappen op `stan` zeer beperkt. +In andere gevallen kan het lonen. + +Weet dat Stan er voor u is, aarzel niet om de uitgebreide documentatie te raadplegen en wees niet bang om het te proberen! + + +## Huiswerk: Hierarchical Model +Om je modelleringsvaardigheden nog verder te ontwikkelen, kun je het "random intercept"-model implementeren en sampelen. + +In "Bayesiaanse" termen is de[algemene terminologie een "hiërarchisch" model](https://mc-stan.org/docs/stan-users-guide/regression.html#hierarchical-regression). + +Hieronder staat de code die sitespecifieke intercepts overweegt. +Het werkt iets anders dan de `+(1|site)`-benadering hierboven: de bovenste is een "offset"-parametrisering, terwijl we deze keer een hyperprior gebruiken. +Dit zijn twee veelvoorkomende strategieën die vaak allebei de moeite waard zijn om te proberen, met marginale of substantiële effecten op de samplerprestaties. + + +```{stan, eval=FALSE, output.var="stan_hierarchical_poisson"} +data { + int N; // sample size + int L; // number of sites + vector[N] habitat_obs; // habitat + array[N] int site_obs; // site + array[N] int sp_rich; // outcome variable: number of ants +} + +parameters { + real intercept_rich; // "global" intercept + vector[L] site_effect; // the sr intercept at each site + real site_variation; // variation on site level + real habitat_effect; // habitat slope +} + +model { + // priors + intercept_rich ~ normal(0, 1); + site_variation ~ cauchy(0, 1); + habitat_effect ~ normal(0, 1); + + // site-wise intercept + for (i in 1:L) { + site_effect[i] ~ normal(intercept_rich, site_variation); + } + + // "posterior", "likelihood", ... give it a name! + sp_rich ~ poisson_log(site_effect[site_obs] + habitat_effect * habitat_obs); + +} +``` + + +```{r, eval=FALSE, stan_hierarchical_sampling} +stan_poisson_fit <- sampling( + stan_poisson_model, + list( + N = nrow(ants_df), + L = length(levels(ants_df$site)), + habitat_obs = as.integer(ants_df$habitat)-1, + site_obs = as.integer(ants_df$site), + sp_rich = ants_df$sp_rich + ), + iter = niter, + chains = nchains, + cores = nparallel + ) + +stan_poisson_fit +``` + +Met Stan zijn de mo(g)e(i)lijkheden bijna eindeloos. +Raak niet verstrikt in het bouwen van modellen! + + + # Referenties {-} - Stan User Guide [(link)](https://mc-stan.org/docs/2_18/stan-users-guide/index.html). diff --git a/content/tutorials/r_brms/index.md b/content/tutorials/r_brms/index.md index 2851711d0..5333996f8 100644 --- a/content/tutorials/r_brms/index.md +++ b/content/tutorials/r_brms/index.md @@ -1,15 +1,18 @@ --- title: "Tutorial Bayesian statistics with brms" description: "This tutorial starts with a theoretical overview of bayesian statistics and the MCMC algorithm. Next, we fit, check, and analyse Bayesian models with the brms package." -author: [raisacarmen, wardlangeraert, toonvandaele] +author: [raisacarmen, wardlangeraert, toonvandaele, falkmielke] date: 2024-01-02 categories: ["r", "statistics"] -tags: ["generalized linear regression", "brms", "r", "mixed models"] +tags: ["generalized linear regression", "brms", "r", "mixed models", "stan"] --- In the fall of 2023, a tutorial on Bayesian statistics with the [**brms**](https://paul-buerkner.github.io/brms/) packages was organised at INBO. -Before you start the tutorial, please follow the instructions in [this R script](https://github.com/inbo/tutorials/blob/master/content/tutorials/r_brms/brms_eng/install_packages.R) to install brms properly. +Before you start the tutorial, please follow the instructions in [this R script](https://github.com/inbo/tutorials/blob/master/static/code/brms_modeling/install_packages.R) to install brms properly. +For more advanced use cases, we added an example of how to directly use [**stan**](https://mc-stan.org) via the `rstan` package. All course material can be found in Dutch and English on the following pages: - The [English tutorial](../../html/workshop_1_mcmc_en_brms_eng.html) with all [course material](https://github.com/inbo/tutorials/tree/master/content/tutorials/r_brms/brms_eng) to test the examples yourself. - The [Dutch tutorial](../../html/workshop_1_mcmc_en_brms.html) with all [course material](https://github.com/inbo/tutorials/tree/master/content/tutorials/r_brms/brms_nl) to test the examples yourself. +- [General files](https://github.com/inbo/tutorials/tree/master/static/code/brms_modeling) useful for both tutorials. + diff --git a/data/authors/falkmielke.toml b/data/authors/falkmielke.toml new file mode 100644 index 000000000..bf4aec20f --- /dev/null +++ b/data/authors/falkmielke.toml @@ -0,0 +1,13 @@ +id = "falkmielke" + +[email] +username = "" +host = "" + +[name] +display = "Falk Mielke" + +[social] +email = "" +facebook = "" +twitter = "" diff --git a/content/tutorials/r_brms/brms_eng/install_packages.R b/static/code/brms_modeling/install_packages.R similarity index 100% rename from content/tutorials/r_brms/brms_eng/install_packages.R rename to static/code/brms_modeling/install_packages.R diff --git a/content/tutorials/r_brms/brms_eng/mcmc_functions.R b/static/code/brms_modeling/mcmc_functions.R similarity index 100% rename from content/tutorials/r_brms/brms_eng/mcmc_functions.R rename to static/code/brms_modeling/mcmc_functions.R diff --git a/static/code/brms_modeling/poisson_model.stan b/static/code/brms_modeling/poisson_model.stan new file mode 100644 index 000000000..34b047e1a --- /dev/null +++ b/static/code/brms_modeling/poisson_model.stan @@ -0,0 +1,21 @@ + +data { + int N; // sample size + vector[N] habitat_obs; // habitat + array[N] int sp_rich; // outcome variable: number of ants +} + +parameters { + real intercept_rich; // intercept + real habitat_effect; // slope +} + +model { + // priors + intercept_rich ~ normal(0, 1); + habitat_effect ~ normal(0, 1); + + // "posterior", "likelihood", ... give it a name! + sp_rich ~ poisson_log(intercept_rich + habitat_effect * habitat_obs); + +} diff --git a/static/html/workshop_1_mcmc_en_brms.html b/static/html/workshop_1_mcmc_en_brms.html index 1bfff1a40..e64871a17 100644 --- a/static/html/workshop_1_mcmc_en_brms.html +++ b/static/html/workshop_1_mcmc_en_brms.html @@ -11,7 +11,7 @@ - + Workshop Bayesiaanse statistiek 2023 @@ -1557,7 +1557,7 @@

Workshop Bayesiaanse statistiek 2023

1. MCMC, model specificatie met brms en interpretatie van de output

Raïsa Carmen, Ward Langeraert & Toon Van Daele

-

2023-11-29

+

2024-11-26

@@ -1572,12 +1572,18 @@

2023-11-29

# Conflicten tussen packages conflicted::conflicts_prefer(dplyr::filter) conflicted::conflicts_prefer(dplyr::lag) +conflicted::conflicts_prefer(tidyr::extract) conflicted::conflicts_prefer(brms::ar) conflicted::conflicts_prefer(brms::dstudent_t) conflicted::conflicts_prefer(brms::pstudent_t) conflicted::conflicts_prefer(brms::qstudent_t) conflicted::conflicts_prefer(brms::rstudent_t) -conflicted::conflict_prefer("rhat", "brms") +conflicted::conflict_prefer("rhat", "brms") + +# html bestandsgrootte verminderen door grafiekparameters te specifiëren +knitr::opts_chunk$set(dev = "jpeg", dpi = 56, + fig.width = 4, fig.height = 3, + out.width = "60%")

1 Theoretische achtergrond

@@ -1682,7 +1688,7 @@

1.1.1 Statistische inferentie

-

+

Stel dat we op één bepaalde site 8 verschillende soorten mieren vinden. Onze likelihood leert ons dat er een oneindig aantal waarden zijn voor \(\lambda\) waarmee we de uitkomst 8 krijgen. De figuur hieronder toont dat de kans op \(X = 8\) verschillend is in ieder van de drie getoonde Poisson verdelingen.

@@ -1695,7 +1701,7 @@

1.1.1 Statistische inferentie

-

+

Ieder van deze modellen, met steeds een andere waarde van \(\lambda\), kan ons de waarde van \(X=8\) geven.

De likelihood functie geeft ons, voor iedere pontentiële waarde van de parameter, de kans op de data. De figuur hieronder toont dit voor onze data met slechts één observatie \(X=8\).

@@ -1711,7 +1717,7 @@

1.1.1 Statistische inferentie

-

+

Bij statistische inferentie willen we de likelihood eigenlijk inverteren \(p(X|\lambda) \rightarrow p(\lambda|X)\). Op deze manier hopen we te weten te komen welk model van alle potentiële modellen het meeste zinvol/waarschijnlijk is. Zowel frequentist als bayesiaanse methoden inverteren in essentie \(p(X|\lambda) \rightarrow p(\lambda|X)\) maar de methode verschilt.

@@ -1760,7 +1766,7 @@

1.2 Parameterschatting in Bayesia p(\lambda|X) = \frac{p(X|\lambda) * p(\lambda)}{p(X)} \]

We vermeldden eerder reeds dat de noemer hier heel moeilijk is om te berekenen. -Er zijn twee mogelijke oplossingen:

+Er zijn drie mogelijke oplossingen:

  • Conjugate priors: Voor sommige verdelingen en combinaties van priors en likelihoods is het mogelijk om toch een analytische oplossing te krijgen (zie deze lijst). Dit beperkt echter onze mogelijkheden:
      @@ -1778,9 +1784,9 @@

      1.2.1 Wat is MCMC?

      Een MCMC algoritme neemt hiervoor random samples van de posterior en genereert een reeks die de hele range van de posterior distributie exploreert. De regels zorgen ervoor dat de sampler de hele tijd in de buurt van de top van de berg blijft. Het is een ‘random walk’ in de parameterruimte in de omgeving waar de posterior distributie het hoogst is (de top van de berg). De grootte en de richting van de stap varieert.

      Bij een stap bergop gaat de MCMC altijd door. Is de stap bergaf, dan gaat de MCMC niet altijd door. De kans om bergaf te gaan hangt af van hoe diep de stap naar beneden is. Kleine stapjes naar beneden worden vaak genomen. Bij diepe stappen naar beneden blijft de MCMC meestal staan en kiest de MCMC vervolgens een nieuwe stap.

      De absolute hoogte van de berg (posterior likelihood) doet er op zich niet toe. Enkel de relatieve verschillen.

      -
      - -

      Metropolis regels

      +
      +Metropolis regels +
      Metropolis regels

      Een eenvoudig MCMC-algoritme is het Metropolis algoritme. Het wordt beschreven door deze formele regels:

        @@ -1811,7 +1817,7 @@

        1.2.1 Wat is MCMC?

        ggplot(df_data, aes(x = x, y = y)) +
   geom_point() +
   geom_smooth(method = "lm", se = TRUE)
        -

        +

        # grootte van de plot beperken

        Een lineaire regressie met LM geeft:

        lm1 <- lm(formula = y ~ x, data = df_data)
@@ -1824,7 +1830,7 @@ 
        1.2.1 Wat is MCMC?
        
 ## (Intercept) 5.173664 9.558508
 ## x           1.032228 2.688598

        We lezen enkele korte functies in om de (log) likelihood en de prior te berekenen en het MCMC metropolis algoritme uit te voeren.

        -
        source(file = "./source/mcmc_functions.R")
        +
        source(file = here("static", "code", "brms_modeling", "mcmc_functions.R"))

        Voor dit eenvoudig model met een kleine dataset kunnen we de posterior voor een groot aantal combinaties voor beta_0 en beta_1 uitrekenen en plotten.

        df <- expand.grid(beta_0 = seq(-1, 15, 0.5), beta_1 = seq(-3, 8, 0.5))
 df$post <- NA
@@ -1837,7 +1843,7 @@ 
        1.2.1 Wat is MCMC?
        
 
        ggplot(df, aes(x = beta_0, y = beta_1, z = post)) +
   geom_raster(aes(fill = post)) + geom_contour(colour = "black") +
   scale_fill_gradientn(colours = rainbow(n = 4))
        
-
        
+
        
 
        Om de MCMC te starten moeten we startwaarden kiezen voor \(\beta_0\) en
 \(\beta_1\)
        
 
        mcmc_small <- MCMC_metro(x = x,         # vector x waarden
@@ -1856,8 +1862,8 @@ 
        1.2.1 Wat is MCMC?
        
              aes(x = beta_0, y = beta_1, z = NA), colour = "red") +
   coord_cartesian(xlim = c(-1, 14.5), ylim = c(-3, 7)) +
   theme(legend.position="none")
        
-
        
-
        De startwaarde ligt een eind van de maximale likelihood. Het duurt even voor de MCMC stabiliseert in de omgeving van de top van de berg. Daarom wordt het eerste deel van de MCMC nooit gebruikt. Dit is de ‘burn-in’ (of warmp-up in brms).
        
+
        
+
        De startwaarde ligt een eind van de maximale likelihood. Het duurt even voor de MCMC stabiliseert in de omgeving van de top van de berg. Daarom wordt het eerste deel van de MCMC nooit gebruikt. Dit is de ‘warmup’ (engels voor opwarmfase; ook ‘burn-in’ of ‘tuning’ in andere bibliotheken).
        
 
        We runnen nu hetzelde model, maar met een veel langere mcmc
        
 
        mcmc_l <- MCMC_metro(x = x,
                      y = y,
@@ -1889,12 +1895,12 @@ 
        1.2.1 Wat is MCMC?
        
 
 gridExtra::grid.arrange(p2, ggplot(), p1, p3, ncol = 2, nrow = 2,
   widths = c(4, 1), heights = c(1, 4))
        
-
        
+
        
 
        Voor elke parameter kan de MCMC worden weergegeven in een ‘trace plot’.
        
 
        ann_text <- tibble(param = c("beta_0", "beta_1"),
                    x = 210,
                    y = c(max(mcmc_l$beta_0), max(mcmc_l$beta_1)),
-                   lab = "burn-in")
+                   lab = "warmup")
 mcmc_l %>%
   dplyr::select(iter, beta_0, beta_1) %>%
   pivot_longer(cols = c("beta_0", "beta_1"), names_to = "param") %>%
@@ -1905,7 +1911,7 @@ 
        1.2.1 Wat is MCMC?
        
   geom_label(data = ann_text, aes(label = lab, x = x, y = y),
              hjust = "left", vjust = "top",
              fill = "gold")
        
-
        
+
        
 
        Schatting voor beta_0 (intercept)
        
 
        quantile(mcmc_l$beta_0[200:5000], probs = c(0.05, 0.5, 0.95))
        
 
        ##       5%      50%      95% 
@@ -1981,11 +1987,11 @@ 
        1.2.1 Wat is MCMC?
        
   ggtitle("equal tail credible interval for a bimodal distribution") +
   geom_line(aes(x = x, y = prob))
 p1 
        
-
        
+
        
 
        p2
        
-
        
+
        
 
        p3
        
-
        
+

@@ -1999,16 +2005,17 @@

1.3 MCMC software

  • RStan is de basis van brms. Je hebt RStannodig om brmste kunnen gebruiken omdat het rechtstreeks met Stan communiceert. Een model definiëren in “Pure” Stan code is echter redelijk omslachtig.
  • rstanarm is nog makkelijker dan brmsqua syntax. Bovendien kan je met rstanarm gebruik maken van enkele vooral gecompileerde Stan modellen. Hierdoor duurt het fitten van je model minder lang. Een groot deel van de procestijd van brmsgaat immers naar het compileren van het Stan programma.
    • -
    - -

    Overzicht van verschillende Stan software (bron: https://jtimonen.github.io/posts/post-01/)

    +
    +Overzicht van verschillende Stan software (bron: https://jtimonen.github.io/posts/post-01/) +
    Overzicht van verschillende Stan software (bron: https://jtimonen.github.io/posts/post-01/)
    +

    Hoewel het pakket brms afhankelijk is van rstan als standaard backend, is er nog een andere softwaretool die het vermelden waard is: cmdstanr. +Vergeleken met rstan profiteert cmdstanr van frequentere updates en enkele performance-voordelen, maar vereist het extra installatiestappen. +In de hieronder geïntroduceerde functie brm(...) kunt u deze twee omwisselen met het trefwoord backend = ....

    -
    -

    2 Een model fitten met brms

    -
    -

    2.1 Dataset laden en data exploratie

    +
    +

    2 Dataset laden en data exploratie

    We laden een dataset in over het aantal mierensoorten in New England (USA). Typ ?ants in de console voor meer info.

    # Laad dataset in
 data(ants)
@@ -2045,7 +2052,7 @@ 
    2.1 Dataset laden en data explora
     geom_bar() +
     scale_x_continuous(limits = c(0, NA)) +
     facet_wrap(~habitat)
    -

    +

    # Boxplots aantal soorten per habitat
 ants_df %>%
   ggplot(aes(y = sp_rich, x = habitat)) +
@@ -2053,7 +2060,7 @@ 
    2.1 Dataset laden en data explora
     geom_line(aes(colour = site, group = site), alpha = 0.5) +
     geom_point(aes(colour = site, group = site)) +
     scale_y_continuous(limits = c(0, NA))
    -

    +

    # Scatter plot aantal soorten en latitude per habitat
 
 ants_df %>%
@@ -2062,7 +2069,7 @@ 
    2.1 Dataset laden en data explora
     geom_smooth(method = "loess", formula = "y ~ x", colour = "firebrick",
                 level = 0.9) +
     facet_wrap(~habitat)
    -

    +

    # Scatter plot aantal soorten en hoogte per habitat
 
 ants_df %>%
@@ -2071,14 +2078,16 @@ 
    2.1 Dataset laden en data explora
     geom_smooth(method = "loess", formula = "y ~ x", colour = "firebrick",
                 level = 0.9) +
     facet_wrap(~habitat)
    -

    +

    Als oefening zullen we een model maken om het aantal soorten te vergelijken tussen beide habitats. Uit de data exploratie zagen we al dat het aantal hoger lijkt te liggen in bossen en dat sites met een hoger aantal in moerassen vaak ook een hoger aantal in bossen hebben.

    -
    -

    2.2 Specificatie van een lineaire regressie

    -
    -

    2.2.1 Model specificatie

    -

    We willen een model waarmee we het verschil kunnen onderzoeken in aantal soorten mieren tussen moeras- en boshabitats. Beschouw de responsvariabele \(Y\), het aantal mieren, en \(X_{habitat}\), een dummy variabele die gelijk is aan 0 voor moerassen en 1 voor bossen. We veronderstellen dat \(Y\) een Normaal verdeling volgt (equivalent aan een lineaire regressie met categorische variabele (= ANOVA))

    +
    +

    3 Een model fitten met brms

    +
    +

    3.1 Specificatie van een lineaire regressie

    +
    +

    3.1.1 Model specificatie

    +

    We willen een model waarmee we het verschil kunnen onderzoeken in aantal soorten mieren tussen moeras- en boshabitats. Beschouw de responsvariabele \(Y\), het aantal mieren, en \(X_{habitat}\), een dummy variabele die gelijk is aan 0 voor moerassen en 1 voor bossen. We veronderstellen dat \(Y\) een Normaal verdeling volgt (equivalent aan een lineaire regressie met categorische variabele (= ANOVA)).

    \[ Y \sim N(\mu, \sigma^2) \]

    @@ -2090,8 +2099,8 @@

    2.2.1 Model specificatie

    Eerst en vooral besluiten we welke MCMC parameters we zullen gebruiken. Typ ?brm om te zien wat de standaard instellingen zijn voor deze parameters. Aangeraden is om meerdere chains te gebruiken (nchains) en deze parallel te laten lopen op verschillende cores van je computer (nparallel). Aangezien dat dit een relatief simpel model is, hebben we niet veel iteraties nodig en geen verdunning.

    # Instellen MCMC parameters
 nchains <- 3           # aantal chains
-niter <- 2000          # aantal iteraties (incl. burn-in, zie volgende)
-burnin <- niter / 4    # aantal initiële samples om te verwijderen (= burn-in)
+niter <- 2000          # aantal iteraties (incl. warmup, zie volgende)
+warmup <- niter / 4    # aantal initiële samples om te verwijderen (= warmup)
 nparallel <- nchains   # aantal cores voor parallel computing
 thinning <- 1          # verdunningsfactor (hier 1 = geen verdunning)

    Het model wordt gefit a.d.h.v. de brm() functie. De syntax is zeer gelijkaardig aan functies die in frequentist statistics worden gebruikt zoals lm() en glm(). De brm() functie bevat heel wat argumenten (zie ?brm). Handige argumenten die we hier niet gebruiken zijn bijvoorbeeld:

    @@ -2100,22 +2109,22 @@

    2.2.1 Model specificatie

  • prior om prior distributies te voorzien. Indien deze niet gespecificeerd zijn, gebruikt de functie default (niet informatieve) priors. Het argument sample_prior kan gebruikt worden om van de prior te samplen zodat je kan visualiseren of de gebruikte prior voldoet aan je verwachtingen (ook handig voor prior predictive checks).
  • file en file_refit om het model object op te slaan nadat het gefit is. Als je de code opnieuw runt en het model is al eens opgeslaan, dan zal brm() dit model gewoon inladen in plaats van het opnieuw te fitten.
    • -
    # Fit Normaal model
+
    # Fit Normaal model
 fit_normal1 <- brm(
   formula = sp_rich ~ habitat, # beschrijving van het model
   family = gaussian(),         # we gebruiken de Normaal verdeling
   data = ants_df,              # ingeven data
   chains = nchains,            # MCMC parameters
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
   seed = 123)                  # seed voor reproduceerbare uitkomst
    
 
    Voor we de resultaten bekijken, controleren we eerst of het model goed convergeert.
    -
    -

    2.2.2 MCMC convergentie

    -

    Er zijn verschillende manieren om convergentie van het MCMC algoritme voor elke parameter te controleren. Hierbij worden de burn-in samples niet in rekening genomen. Eerst en vooral heb je visuele controles.

    +
    +

    3.1.2 MCMC convergentie

    +

    Er zijn verschillende manieren om convergentie van het MCMC algoritme voor elke parameter te controleren. Hierbij worden de warmup samples niet in rekening genomen. Eerst en vooral heb je visuele controles.

    We kunnen de MCMC samples met de as_draws() functies verkrijgen ofwel ineens visualiseren via de bayesplot package die compatibel is met brmsfit objecten.

    # Zet kleurenpallet voor duidelijke visualisatie in bayesplot
 color_scheme_set("mix-blue-red")
    @@ -2132,10 +2141,10 @@

    2.2.2 MCMC convergentie

    ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) + geom_line() + facet_wrap(~parameter, nrow = 3, scales = "free")     -

    +

    # Visualisatie via Bayesplot package
 mcmc_trace(fit_normal1, pars = parameters)
    -

    +

    running mean/quantile/… plot

    Naast een trace plot, kunnen we ook kijken hoe bepaalde statistieken zoals het gemiddelde of de kwantielen variëren over de tijd (= met toename aantal iteraties). Je wilt natuurlijk dat deze statistieken stabiliseren na een aantal iteraties want dan weet je dat je genoeg iteraties gebruikt hebt.

    # code om cumulatieve kwantielen te berekenen
@@ -2159,21 +2168,21 @@ 
    2.2.2 MCMC convergentie
    
   pivot_longer(cols = starts_with(parameters), names_to = "name",
                values_to = "value") %>%
   # splits naam kolom op bij laatste underscore
-  extract(name, into = c("parameter", "statistiek"), "(.*)_([^_]+)$") %>%
+  tidyr::extract(name, into = c("parameter", "statistiek"), "(.*)_([^_]+)$") %>%
   ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) +
     geom_line(aes(linetype = statistiek), linewidth = 0.8) +
     facet_wrap(~parameter, nrow = 3, scales = "free")
    -

    +

    posterior density plot

    De vorm van de posterior distributions kan ook informatief zijn om na te kijken of de chains geconvergeerd zijn naar dezelfde distributies. Als we een bimodale verdeling zien (een kamelenrug met twee pieken), dan zou er mogelijks ook iets mis kunnen zijn met de specificatie van ons model.

    # Visualisatie density plot van de posterior voor ieder van de chains apart in 
 # overlay. via Bayesplot package
 mcmc_dens_overlay(fit_normal1, pars = parameters)
    -

    +

    Nog eenvoudiger is de plot() functie gebruiken op het brmsfit object. Deze toont direct de trace plots en de posterior density plots naast elkaar voor elke parameter.

    # trace plots en posterior density voor iedere parameter
 plot(fit_normal1)
    -

    +

    autocorrelation plot

    Als twee variabelen gecorreleerd zijn, wil dat zeggen dat er een verband is tussen beiden. Dit kan een positief of negative verband zijn. Een voorbeeld van positief gecorreleerde variabelen is de grootte en het gewicht van mensen; hoe groter iemand is, hoe meer hij typisch zal wegen. Een voorbeeld van negatief gecorreleerde variabelen is de hoeveelheid zonneschijn en het vochtgehalte in de toplaag van de bodem.

    @@ -2190,7 +2199,7 @@

    2.2.2 MCMC convergentie

    Dit zal de autocorrelatie reduceren maar ook het aantal posterior samples die we overhouden.

    # Visualisatie autocorrelation plots via Bayesplot package
 mcmc_acf(fit_normal1, pars = parameters)
    -

    +

    Naast visuele checks zijn er ook diagnostische controles.

    Gelman-Rubin diagnostic

    Ook wel \(\hat{R}\) (R-hat) of potential scale reduction factor genoemd. Een manier om te controleren of een chain is geconvergeerd, is door zijn gedrag te vergelijken met andere willekeurig geïnitialiseerde ketens. De \(\hat{R}\) statistiek neemt de verhouding van de gemiddelde variantie van samples binnen elke chain tot de variantie van de gepoolde samples over alle chains heen. Als alle ketens convergeren naar een gemeenschappelijke distributie, zullen deze verhoudingen gelijk aan 1 zijn. Als de ketens niet zijn geconvergeerd naar een gemeenschappelijke verdeling, zal de \(\hat{R}\) statistiek groter zijn dan 1. We kunnen de \(\hat{R}\) statistieken extraheren via de rhat() functie en dan eventueel visualiseren met de mcmc_rhat() functie van de bayesplot package.

    @@ -2198,42 +2207,42 @@

    2.2.2 MCMC convergentie

    rhats_fit_normal1 <- rhat(fit_normal1)[parameters] print(rhats_fit_normal1) 
    ##     b_Intercept b_habitatForest           sigma 
-##       1.0003861       1.0004918       0.9998149
    +## 0.999769 1.000314 1.000759     
    # Plot R-hats via Bayesplot package
 mcmc_rhat(rhats_fit_normal1) + yaxis_text(hjust = 1)
    -

    +

    effective sample size

    De effective sample size (\(n_{eff}\)) is een schatting van het aantal onafhankelijke trekkingen van de posterior distribution van elke parameter. Omdat de trekkingen binnen een Markov chain niet onafhankelijk zijn als er sprake is van autocorrelatie, is de effective sample size gewoonlijk kleiner dan de totale steekproefomvang, namelijk het aantal iteraties (\(N\)). Hoe groter de verhouding \(n_{eff}/N\) hoe beter. De bayesplot package biedt een neff_ratio() extractorfunctie. De mcmc_neff() functie kan vervolgens worden gebruikt om de verhoudingen te plotten.

    # Krijg verhoudingen effective sample size
 ratios_fit_normal1 <- neff_ratio(fit_normal1)[parameters]
 print(ratios_fit_normal1)
    ##     b_Intercept b_habitatForest           sigma 
-##       0.7436946       0.6869533       0.7381539
    +## 0.7153819 0.6962991 0.6844347 
    # Plot verhoudingen via Bayesplot package
 mcmc_neff(ratios_fit_normal1) + yaxis_text(hjust = 1)
    -

    +

    Andere packages die van pas kunnen komen bij het controleren van MCMC convergentie zijn mcmcplots en coda. Voor deze packages moet je de MCMC samples uit je brms model wel eerst omzetten naar een mcmc object.

    # voorbeeld
 install.packages("mcmcplots")
 mcmcplots::mcmcplot(as.mcmc(fit_normal1, pars = parameters))

    We kunnen over het algemeen besluiten dat MCMC convergentie goed is voor alle parameters. Indien dit niet het geval zou zijn, kan het zijn dat de MCMC parameters anders moeten gekozen worden, bijvoorbeeld door het aantal iteraties te verhogen. Nu we weten dat de parameters correct geschat zijn, kunnen we controleren of het model goed bij de data past.

    -
    -

    2.2.3 Model fit

    +
    +

    3.1.3 Model fit

    De posterior predictive check (PPC) is een goede manier om te controleren of het model goed past bij de data. Het idee achter de PPC is dat, als een model goed past, de voorspelde waardes o.b.v. het model veel lijken op de data die we hebben gebruikt om het model te fitten. Om de voorspellingen te genereren die worden gebruikt voor PPCs, simuleren we meerdere malen vanuit de posterior predictive distribution (de posterior distribution van de responsvariabele). Visualisatie kan met de bayesplot package. De dikke lijn toont de data en de dunne lijntjes tonen verschillende simulaties o.b.v. het model.

    # Visualiseer model fit via Bayesplot package
 pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    -

    +

    We zien dat de geobserveerde data en de voorspellingen o.b.v. ons model niet helemaal overeenkomen. Mogelijks was onze keuze voor de Normaal verdeling niet goed. Aantallen zijn discrete, positieve waarden, terwijl de Normaal verdeling continue en zowel positieve als negatieve waarden kan voorspellen. De Poisson verdeling is een discrete verdeling die steeds positief is. Ze wordt vaak gebruikt voor het modelleren van geregistreerde aantallen gedurende een gegeven tijdsinterval, afstand, oppervlakte, volume …

    Opmerking:

    Merk op dat MCMC convergentie en model fit twee verschillende dingen zijn. Het is niet omdat MCMC convergentie goed is dat model fit ook goed is. MCMC convergentie is om na te gaan of het MCMC algoritme de distributies van de parameters goed heeft kunnen benaderen. Model fit controleert of het model dat door ons aangenomen is (assumptie normaliteit, lineair verband gemiddelde …) goed past bij de data.

    -
    -

    2.3 Specificatie van een Poisson model

    -
    -

    2.3.1 Model specificatie

    +
    +

    3.2 Specificatie van een Poisson model

    +
    +

    3.2.1 Model specificatie

    We veronderstellen dat \(Y\) nu een Poisson verdeling volgt

    \[ Y \sim Pois(\lambda) @@ -2243,39 +2252,39 @@

    2.3.1 Model specificatie

    \ln(\lambda) = \beta_0 + \beta_1X_{habitat} \]

    We moeten dus twee parameters schatten: \(\beta_0\) en \(\beta_1\) We gebruiken dezelfde MCMC parameters als voordien. Het enige wat we moeten aanpassen is de keuze family = poisson().

    -
    # Fit Poisson model
+
    # Fit Poisson model
 fit_poisson1 <- brm(
   formula = sp_rich ~ habitat, # beschrijving van het model
   family = poisson(),          # we gebruiken de Poisson verdeling
   data = ants_df,              # ingeven data
   chains = nchains,            # MCMC parameters
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
   seed = 123)                  # seed voor reproduceerbare uitkomst
    -
    -

    2.3.2 MCMC convergentie

    +
    +

    3.2.2 MCMC convergentie

    Convergentie ziet er opnieuw goed uit.

    # Visualiseer MCMC convergentie van het Poisson model
 plot(fit_poisson1)
    -

    +

    # Extraheer en visualiseer R-hat waarden
 rhats_fit_poisson1 <- rhat(fit_poisson1)[c("b_Intercept", "b_habitatForest")]
 mcmc_rhat(rhats_fit_poisson1) + yaxis_text(hjust = 1)
    -

    +

    -
    -

    2.3.3 Model fit

    +
    +

    3.2.3 Model fit

    We zien dat alle predicties nu wel positief zijn, maar de fit is nog altijd niet ideaal.

    # Visualiseer model fit via Bayesplot package
 pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    -

    +

    -
    -

    2.3.4 Discussie

    +
    +

    3.2.4 Discussie

    Hoe kunnen we model fit nog verbeteren?

    We weten dat elke site 2x bezocht is. Eenmaal in moeras en eenmaal in bos. Het is natuurlijk mogelijk dat het aantal mieren in het moeras en het bos van dezelfde site gecorreleerd zullen zijn (site effect). Dit zagen we inderdaad ook in de data exploratie. Sites met een hoger aantal soorten in moerassen hebben vaak ook een hoger aantal in bossen en vice versa. Hiervoor kunnen we corrigeren door een random intercept voor elke site toe te voegen ... + (1|site).

    We veronderstellen dat het aantal soorten \(Y\) per site \(j\) (\(j = 1, ..., J\)) een Poisson verdeling volgt

    @@ -2290,13 +2299,13 @@

    2.3.4 Discussie

    \[ b_0 \sim N(0, \sigma_b) \]

    -
    # Fit Poisson model met random intercepten per site
+
    # Fit Poisson model met random intercepten per site
 fit_poisson2 <- brm(
   formula = sp_rich ~ habitat + (1|site),
   family = poisson(),
   data = ants_df,
   chains = nchains,
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
@@ -2304,31 +2313,30 @@ 
    2.3.4 Discussie
    
 
    Convergentie is goed.
    
 
    # Visualiseer MCMC convergentie
 plot(fit_poisson2)
    
-
    
+
    
 
    # Extraheer en visualiseer R-hat waarden
 parameters2 <- c("b_Intercept", "b_habitatForest", "sd_site__Intercept")
 rhats_fit_poisson2 <- rhat(fit_poisson2)[parameters2]
 mcmc_rhat(rhats_fit_poisson2) + yaxis_text(hjust = 1)
    
-
    
+
    
 
    Model fit lijkt nu nog beter.
    
 
    # Visualiseer model fit van het Poisson model met random intercept via
 # Bayesplot package
 pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    
-
    
+
    
 
    Hoe kunnen we deze modellen nu objectief gaan vergelijken?
    -
    -
    -

    3 Vergelijken van modellen

    +
    +

    3.3 Vergelijken van modellen

    Op basis van de PPCs kunnen we reeds zien welk model het best past bij de data. Verder zijn er nog enkele functies die brms voorziet om verschillende modellen te vergelijken. Met de functie add_criterion() kan je model fit criteria toevoegen aan model objecten. Typ ?add_criterion() om te zien welke beschikbaar zijn. Zie ook https://mc-stan.org/loo/articles/online-only/faq.html

    -
    -

    3.1 Leave-one-out cross-validation

    +
    +

    3.3.1 Leave-one-out cross-validation

    Cross-validation (CV) is een familie van technieken die probeert in te schatten hoe goed een model onbekende data zou voorspellen via predicties van het model gefit op de bekende data. Hiervoor moet je niet per se nieuwe data gaan inzamelen. Je kan jouw eigen data opsplitsen in een test en training dataset. Je fit het model op de training dataset en je gebruikt dat model dan om te schatten hoe goed het de data in de test dataset kan voorspellen. Bij leave-one-out CV (LOOCV) ga je telkens één observatie weglaten en het model opnieuw fitten o.b.v. alle andere observaties.

    -
    - -

    Illustratie LOOCV.

    +
    +Illustratie LOOCV. +
    Illustratie LOOCV.
    # Voeg leave-one-out model fit criterium toe aan model objecten
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2343,13 +2351,13 @@ 
    3.1 Leave-one-out cross-validatio
                         criterion = "loo")
 print(comp_loo, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_loo se_elpd_loo p_loo    se_p_loo looic   
-## fit_poisson2    0.000     0.000 -106.276    3.849      11.552    1.522  212.552
-## fit_poisson1  -13.403     5.802 -119.679    8.319       3.714    0.888  239.358
-## fit_normal1   -15.647     3.629 -121.923    5.301       3.013    0.842  243.846
+## fit_poisson2    0.000     0.000 -106.077    3.896      11.399    1.569  212.154
+## fit_poisson1  -13.685     5.728 -119.762    8.309       3.813    0.903  239.524
+## fit_normal1   -15.874     3.597 -121.951    5.329       3.007    0.841  243.902
 ##              se_looic
-## fit_poisson2    7.698
-## fit_poisson1   16.637
-## fit_normal1    10.601
    +## fit_poisson2 7.792 +## fit_poisson1 16.618 +## fit_normal1 10.659

    Als je ervan uitgaat dat dit verschil normaal verdeeld is, kan je met de onzekerheid een betrouwbaarheidsinterval berekenen. We zien dat het tweede Poisson model het best scoort en dat de andere modellen significant lager scoren.

    # Bereken betrouwbaarheidsinterval om modellen te vergelijken met loocv
 comp_loo %>%
@@ -2359,15 +2367,15 @@ 
    3.1 Leave-one-out cross-validatio
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    ##              elpd_diff  se_diff   ll_diff   ul_diff
 ## fit_poisson2   0.00000 0.000000   0.00000  0.000000
-## fit_poisson1 -13.40293 5.802467 -22.94714 -3.858719
-## fit_normal1  -15.64659 3.628604 -21.61511 -9.678068
    +## fit_poisson1 -13.68487 5.727656 -23.10603 -4.263713 +## fit_normal1 -15.87376 3.596654 -21.78973 -9.957793
    -
    -

    3.2 K-fold cross-validation

    +
    +

    3.3.2 K-fold cross-validation

    Bij K-fold cross-validation worden de data in \(K\) groepen opgesplitst. Wij zullen hier \(K = 10\) groepen (= folds) gebruiken. In plaats van dus telkens één enkele observatie weg te laten zoals bij leave-one-out CV gaan we hier \(1/10\)e van de data weglaten. Via de argumenten folds = "stratified" en group = "habitat" zorgen we ervoor dat voor elke groep de relatieve frequenties van habitat bewaard blijven. Deze techniek zal dus minder precies zijn dan de vorige, maar zal sneller zijn om te berekenen indien je met heel veel data werkt.

    -
    - -

    Illustratie K-fold CV.

    +
    +Illustratie K-fold CV. +
    Illustratie K-fold CV.
    # Voeg K-fold model fit criterium toe aan model objecten
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2391,9 +2399,9 @@ 
    3.2 K-fold cross-validation
    
                           criterion = "kfold")
 print(comp_kfold, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_kfold se_elpd_kfold p_kfold  se_p_kfold
-## fit_poisson2    0.000     0.000 -106.773      4.145        12.048    1.944  
-## fit_poisson1  -12.290     5.195 -119.063      8.156         3.097    0.939  
-## fit_normal1   -14.619     3.258 -121.392      5.172         2.482    0.796
    +## fit_poisson2 0.000 0.000 -106.516 3.890 11.837 1.561 +## fit_poisson1 -12.646 5.808 -119.162 8.227 3.212 1.025 +## fit_normal1 -15.101 3.616 -121.617 5.148 2.672 0.717

    De resultaten zijn hetzelfde, maar de verschillen zijn iets kleiner als voordien. Het verschil tussen de twee Poisson modellen is niet meer significant.

    # Bereken betrouwbaarheidsinterval om modellen te vergelijken met K-fold CV
 comp_kfold %>%
@@ -2403,11 +2411,11 @@ 
    3.2 K-fold cross-validation
    
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    ##              elpd_diff  se_diff   ll_diff   ul_diff
 ## fit_poisson2   0.00000 0.000000   0.00000  0.000000
-## fit_poisson1 -12.28990 5.195211 -20.83526 -3.744537
-## fit_normal1  -14.61938 3.258315 -19.97883 -9.259928
    +## fit_poisson1 -12.64595 5.807731 -22.19882 -3.093082 +## fit_normal1 -15.10095 3.616047 -21.04882 -9.153085
    -
    -

    3.3 WAIC

    +
    +

    3.3.3 WAIC

    Het Widely Applicable Information Criterion (WAIC) maakt geen gebruik van cross-validation maar is een computationele manier om de ELPD te schatten. Hoe dit precies gebeurt valt buiten het doel van deze workshop. Het is nog een andere maat om model selectie toe te passen.

    # Voeg WAIC model fit criterium toe aan model objecten
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2422,29 +2430,30 @@ 
    3.3 WAIC
    
                          criterion = "waic")
 print(comp_waic, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_waic se_elpd_waic p_waic   se_p_waic
-## fit_poisson2    0.000     0.000 -105.062     3.683       10.338    1.348 
-## fit_poisson1  -14.594     5.858 -119.656     8.312        3.691    0.882 
-## fit_normal1   -16.833     3.624 -121.895     5.289        2.985    0.830 
+## fit_poisson2    0.000     0.000 -104.809     3.678       10.131    1.344 
+## fit_poisson1  -14.934     5.851 -119.743     8.305        3.794    0.899 
+## fit_normal1   -17.104     3.620 -121.913     5.313        2.969    0.824 
 ##              waic     se_waic 
-## fit_poisson2  210.124    7.367
-## fit_poisson1  239.313   16.624
-## fit_normal1   243.790   10.578
    +## fit_poisson2 209.618 7.355 +## fit_poisson1 239.486 16.609 +## fit_normal1 243.827 10.626 
    # Bereken betrouwbaarheidsinterval om modellen te vergelijken met WAIC
 comp_waic %>%
   as.data.frame() %>%
   select(waic, elpd_diff, se_diff) %>%
   mutate(ll_diff = elpd_diff  + qnorm(0.05) * se_diff,
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    -
    ##                  waic elpd_diff  se_diff   ll_diff    ul_diff
-## fit_poisson2 210.1239   0.00000 0.000000   0.00000   0.000000
-## fit_poisson1 239.3128 -14.59444 5.857609 -24.22935  -4.959532
-## fit_normal1  243.7901 -16.83311 3.624302 -22.79456 -10.871665
    +
    ##                  waic elpd_diff  se_diff   ll_diff   ul_diff
+## fit_poisson2 209.6176   0.00000 0.000000   0.00000   0.00000
+## fit_poisson1 239.4857 -14.93407 5.850998 -24.55810  -5.31003
+## fit_normal1  243.8265 -17.10447 3.620256 -23.05926 -11.14968
    -
    -

    3.4 Conclusie

    +
    +

    3.3.4 Conclusie

    Zowel o.b.v. de PPC als vergelijkingen met verschillende model selectie criteria, kunnen we besluiten dat het tweede Poisson model met random intercepts het best past bij de data. In principe konden we dit ook verwachten op basis van onze eigen intuïtie en het design van de studie, nl. het gebruik van de Poisson distributie om aantallen te modelleren en het gebruik van random intercepts om te controleren voor een hiërarchisch design (habitats genest in sites).

    +

    4 Resultaten finale model

    Als we het model fit object bekijken, zien we resultaten die vergelijkbaar zijn met resultaten zoals we die zien als we een frequentist model gefit hebben. We krijgen enerzijds een schatting van alle parameters met hun onzekerheid maar anderzijds zien we dat dit duidelijk de output is van een Bayesiaans model. Zo krijgen we info over de parameters die we gebruikt hebben voor het MCMC algoritme, we krijgen een 95 % credible interval (CI) in plaats van een confidence interval en we krijgen bij elke parameter ook de \(\hat{R}\) waarde die we eerder besproken hebben.

    @@ -2457,15 +2466,15 @@

    4 Resultaten finale model

    ## Draws: 3 chains, each with iter = 2000; warmup = 500; thin = 1; ## total post-warmup draws = 4500 ## -## Group-Level Effects: +## Multilevel Hyperparameters: ## ~site (Number of levels: 22) ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS -## sd(Intercept) 0.38 0.10 0.22 0.60 1.00 1748 2738 +## sd(Intercept) 0.38 0.09 0.23 0.60 1.00 1733 2561 ## -## Population-Level Effects: +## Regression Coefficients: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS -## Intercept 1.51 0.13 1.25 1.77 1.00 2873 2723 -## habitatForest 0.64 0.12 0.41 0.88 1.00 7102 3251 +## Intercept 1.51 0.13 1.24 1.75 1.00 2452 2883 +## habitatForest 0.64 0.12 0.41 0.87 1.00 5062 3193 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential @@ -2482,14 +2491,14 @@

    4 Resultaten finale model

    q_20 = quantile(.value, probs = 0.20), gemiddelde = mean(.value), mediaan = median(.value), - q_20 = quantile(.value, probs = 0.80), + q_80 = quantile(.value, probs = 0.80), q_95 = quantile(.value, probs = 0.95), max = max(.value)) -
    ## # A tibble: 2 × 8
-##   .variable         min  q_05  q_20 gemiddelde mediaan  q_95   max
-##   <chr>           <dbl> <dbl> <dbl>      <dbl>   <dbl> <dbl> <dbl>
-## 1 b_Intercept     1.04  1.29  1.63       1.51    1.52  1.74   1.97
-## 2 b_habitatForest 0.268 0.432 0.742      0.638   0.639 0.831  1.00
    +
    ## # A tibble: 2 × 9
+##   .variable         min  q_05  q_20 gemiddelde mediaan  q_80  q_95   max
+##   <chr>           <dbl> <dbl> <dbl>      <dbl>   <dbl> <dbl> <dbl> <dbl>
+## 1 b_Intercept     1.06  1.29  1.41       1.52    1.52  1.62  1.72   1.95
+## 2 b_habitatForest 0.275 0.454 0.544      0.641   0.639 0.742 0.827  1.07

    Handige functies van de tidybayes package zijn ook median_qi(), mean_qi() … na gather_draws() die je kan gebruiken in plaats van group_by() en summarise().

    We zouden nu graag het geschatte aantal soorten visualiseren per habitattype met bijbehorende onzekerheid. Met de functie spread_draws() kan je een bepaald aantal samples van de posterior distributie nemen en omzetten in een tabel in wijd formaat. Het gemiddeld aantal soorten in moerassen volgens ons model is \(\exp(\beta_0)\) en in bossen \(\exp(\beta_0+\beta_1)\). We tonen de posterior distributions met de posterior mediaan en 60 en 90 % credible intervals.

    fit_poisson2 %>%
@@ -2504,7 +2513,7 @@ 
    4 Resultaten finale model
    
   ggplot(aes(y = sp_rich, x = habitat)) +
     stat_eye(point_interval = "median_qi", .width = c(0.6, 0.9)) +
     scale_y_continuous(limits = c(0, NA))
    -

    +

    Naast stat_eye() vind je hier nog leuke manieren om posterior distributions te visualiseren.

    We zien een duidelijk verschil in aantal soorten tussen beide habitats. Is er een significant verschil tussen het aantal soorten in moerassen en bossen? We testen de hypothese dat het aantal in moerassen en bossen gelijk is

    \[ @@ -2520,7 +2529,7 @@

    4 Resultaten finale model

    hyp 
    ## Hypothesis Tests for class b:
 ##            Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
-## 1 (habitatForest) = 0     0.64      0.12     0.44     0.84         NA        NA
+## 1 (habitatForest) = 0     0.64      0.12     0.45     0.83         NA        NA
 ##   Star
 ## 1    *
 ## ---
@@ -2530,7 +2539,7 @@ 
    4 Resultaten finale model
    
 ## Posterior probabilities of point hypotheses assume equal prior probabilities.
    # Plot posterior distribution hypothese
 plot(hyp)
    -

    +

    Finaal visualiseren we de random effecten van de sites. We sorteren ze van hogere soortenrijkdom naar lager.

    # Neem het gemiddelde van sd van random effects
 # om verderop aan figuur toe te voegen
@@ -2550,7 +2559,7 @@ 
    4 Resultaten finale model
    
                color = "darkgrey", linetype = 2) +
     stat_halfeye(point_interval = "median_qi", .width = 0.9, size = 2/3,
                  fill = "cornflowerblue")
    -

    +

    5 Vergelijking met frequentist statistics

    @@ -2569,14 +2578,14 @@

    5 Vergelijking met frequentist st ## Draws: 3 chains, each with iter = 2000; warmup = 500; thin = 1; ## total post-warmup draws = 4500 ## -## Population-Level Effects: +## Regression Coefficients: ## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS -## Intercept 4.84 0.80 3.53 6.16 1.00 4194 3347 -## habitatForest 4.33 1.14 2.44 6.20 1.00 4113 3091 +## Intercept 4.80 0.81 3.46 6.11 1.00 4196 3219 +## habitatForest 4.33 1.14 2.42 6.15 1.00 4502 3133 ## -## Family Specific Parameters: +## Further Distributional Parameters: ## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS -## sigma 3.73 0.41 3.13 4.44 1.00 3925 3322 +## sigma 3.73 0.41 3.12 4.46 1.00 3557 3080 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential @@ -2599,7 +2608,142 @@

    5 Vergelijking met frequentist st ## sample estimates: ## mean in group Bog mean in group Forest ## 4.863636 9.181818 -

    We zien dat dit inderdaad bijna exact dezelfde resultaten oplevert. Ons Bayesiaans model schat dat er gemiddeld 4.326 meer mierensoorten in bossen voorkomen dan in moerassen (90 % credible interval: 2.439 tot 6.203). De t-test schat dat er gemiddeld 4.318 meer mierensoorten in bossen voorkomen dan in moerassen (90 % confidence interval: 2.445 tot 6.192).

    +

    We zien dat dit inderdaad bijna exact dezelfde resultaten oplevert. Ons Bayesiaans model schat dat er gemiddeld 4.329 meer mierensoorten in bossen voorkomen dan in moerassen (90 % credible interval: 2.423 tot 6.149). De t-test schat dat er gemiddeld 4.318 meer mierensoorten in bossen voorkomen dan in moerassen (90 % confidence interval: 2.445 tot 6.192).

    +

    +
    +

    6 Deep Dive: rstan

    +
    +

    6.1 Stan: Wat? Waarom?!

    +

    De brms package is een gemaksverpakking van rstan, dat zelf de functionaliteit van stan naar R brengt. +Stan is een modeleringsomgeving, geschreven in de programmeertaal C, die de benodigdheden voor probabilistische (“Bayesiaanse”) modellen implementeerd.

    +

    Meer info vind je op het internet.

    +

    Het voordeel van brms is bruikbaarheid: veel functies werken “out-of-the-box”, met redelijke standaardwaarden en een syntaxis die lijkt op wat veel statistici gewend zijn. +Het relatieve gebruiksgemak kan echter ten koste gaan van de flexibiliteit en in zekere mate van de leesbaarheid.

    +

    Stan en rstan daarentegen leunen meer op de wiskundige formulering van modellen. +Elk aspect van het model moet expliciet worden ingesteld, wat een voordeel kan zijn (bijvoorbeeld als u te maken hebt met niet-standaard usecases) of een nadeel (bijvoorbeeld als u modellen op niet-optimale manieren specificeert).

    +

    Om kort een indruk te geven, bouwen we dezelfde modellen als hierboven, met behulp van het Stan-framework.

    +
    library("rstan")
+conflicted::conflicts_prefer(rstan::extract)
+conflicted::conflicts_prefer(brms::loo)
    +
    +
    +

    6.2 Model Definition

    +

    RMarkdown kan stan-chunks verwerken, hoewel de algemenere modeldefinitie is uitbesteed aan een apart “*.stan”-bestand (bv. hier). +Alternatief kan je je model ook als een groot tekstblok definiëren, zoals hieronder te zien. +Het eenvoudige poissonmodel lijkt op een van de stan-daardvoorbeelden, waarnaar u terecht kunt voor verdere details en meer.

    +
    stan_poisson_model_code = '
+data {
+  int<lower=1> N;                   // sample size
+  vector[N] habitat_obs;            // habitat
+  array[N] int<lower=0> sp_rich;    // outcome variable: number of ants
+}
+
+parameters {
+  real intercept_rich;              // intercept
+  real habitat_effect;              // slope
+}
+
+model {
+  // priors
+  intercept_rich ~ normal(0, 1);
+  habitat_effect ~ normal(0, 1);
+
+  // "posterior", "likelihood", ... give it a name!
+  sp_rich ~ poisson_log(intercept_rich + habitat_effect * habitat_obs);
+
+}
+
+'
    +

    Zie dit als een “kijkje achter de schermen”! +Je moet expliciet de modelstructuur, priors en zelfs de gegevenstypen van invoervariabelen definiëren. +Maar let op hoe zelfs Stan niet zonder handige functies is: we kunnen de poisson_log posterior gebruiken om de \(log(\lambda )\) te modelleren.

    +

    Als je buiten RStudio/RMarkdown werkt, kun je het model het beste laden vanuit een bestand:

    +
    stan_poisson_model <- stan_model(
+  # file = here("code", "brms_modeling", "poisson_model.stan"),
+  model_code = stan_poisson_model_code, 
+  model_name = "stan poisson model"
+  )
    +
    +
    +

    6.3 Sampling

    +

    Sampling doet vrijwel hetzelfde als hierboven, want in de kern is brms gewoon stan.

    +
    stan_poisson_fit <- sampling(
+  stan_poisson_model,
+  list(
+    N = nrow(ants_df),
+    habitat_obs = as.integer(ants_df$habitat)-1,
+    sp_rich = ants_df$sp_rich
+    ),
+  iter = niter,
+  chains = nchains,
+  cores = nparallel
+  )
+
+stan_poisson_fit
    +

    Een handige manier om modelfits snel te inspecteren is shinystan. +Het opent een browser met glanzende plots.

    +
    library("shinystan")
+launch_shinystan(stan_poisson_fit)
    +

    En kijk: de modeluitkomst is precies zoals boven. +In het huidige geval is het voordeel van het overstappen op stan zeer beperkt. +In andere gevallen kan het lonen.

    +

    Weet dat Stan er voor u is, aarzel niet om de uitgebreide documentatie te raadplegen en wees niet bang om het te proberen!

    +
    +
    +

    6.4 Huiswerk: Hierarchical Model

    +

    Om je modelleringsvaardigheden nog verder te ontwikkelen, kun je het “random intercept”-model implementeren en sampelen.

    +

    In “Bayesiaanse” termen is dealgemene terminologie een “hiërarchisch” model.

    +

    Hieronder staat de code die sitespecifieke intercepts overweegt. +Het werkt iets anders dan de +(1|site)-benadering hierboven: de bovenste is een “offset”-parametrisering, terwijl we deze keer een hyperprior gebruiken. +Dit zijn twee veelvoorkomende strategieën die vaak allebei de moeite waard zijn om te proberen, met marginale of substantiële effecten op de samplerprestaties.

    +
    data {
+  int<lower=1> N;                   // sample size
+  int<lower=1> L;                   // number of sites
+  vector[N] habitat_obs;            // habitat
+  array[N] int<lower=1, upper=L> site_obs; // site
+  array[N] int<lower=0> sp_rich;    // outcome variable: number of ants
+}
+
+parameters {
+  real intercept_rich;              // "global" intercept
+  vector[L] site_effect;            // the sr intercept at each site
+  real<lower=0> site_variation;     // variation on site level
+  real habitat_effect;              // habitat slope
+}
+
+model {
+  // priors
+  intercept_rich ~ normal(0, 1);
+  site_variation ~ cauchy(0, 1);
+  habitat_effect ~ normal(0, 1);
+
+  // site-wise intercept
+  for (i in 1:L) {
+    site_effect[i] ~ normal(intercept_rich, site_variation);
+  }
+
+  // "posterior", "likelihood", ... give it a name!
+  sp_rich ~ poisson_log(site_effect[site_obs] + habitat_effect * habitat_obs);
+
+}
    +
    stan_poisson_fit <- sampling(
+  stan_poisson_model,
+  list(
+    N = nrow(ants_df),
+    L = length(levels(ants_df$site)),
+    habitat_obs = as.integer(ants_df$habitat)-1,
+    site_obs = as.integer(ants_df$site),
+    sp_rich = ants_df$sp_rich
+    ),
+  iter = niter,
+  chains = nchains,
+  cores = nparallel
+  )
+
+stan_poisson_fit
    +

    Met Stan zijn de mo(g)e(i)lijkheden bijna eindeloos. +Raak niet verstrikt in het bouwen van modellen!

    +

    Referenties

    diff --git a/static/html/workshop_1_mcmc_en_brms_eng.html b/static/html/workshop_1_mcmc_en_brms_eng.html index d44879cd9..7c1831040 100644 --- a/static/html/workshop_1_mcmc_en_brms_eng.html +++ b/static/html/workshop_1_mcmc_en_brms_eng.html @@ -11,7 +11,7 @@ - + Tutorial Bayesian statistics @@ -1557,7 +1557,7 @@

    Tutorial Bayesian statistics

    1. MCMC, model specification with brms and output interpretation

    Raïsa Carmen, Ward Langeraert & Toon Van Daele

    -

    2024-01-02

    +

    2024-11-26

    @@ -1572,12 +1572,18 @@

    2024-01-02

    # Conflicten tussen packages conflicted::conflicts_prefer(dplyr::filter) conflicted::conflicts_prefer(dplyr::lag) +conflicted::conflicts_prefer(tidyr::extract) conflicted::conflicts_prefer(brms::ar) conflicted::conflicts_prefer(brms::dstudent_t) conflicted::conflicts_prefer(brms::pstudent_t) conflicted::conflicts_prefer(brms::qstudent_t) conflicted::conflicts_prefer(brms::rstudent_t) -conflicted::conflict_prefer("rhat", "brms") +conflicted::conflict_prefer("rhat", "brms") + +# reduce html export file size by specifying graphics size +knitr::opts_chunk$set(dev = "jpeg", dpi = 56, + fig.width = 4, fig.height = 3, + out.width = "60%")

    1 Theoretical background

    @@ -1682,7 +1688,7 @@

    1.1.1 Statistical inference

    geom_errorbar(aes(ymin = 0, ymax = prob), width = 0) + xlab("number of ant species at a site") + ylab("probability mass") -

    +

    Suppose we find 8 different species of ants on one particular site. Our likelihood tells us that there are an infinite number of values for \(\lambda\), which gives us the result 8. The figure below shows that the probability of \(X = 8\) is different in each of the three shown Poisson distributions.

    @@ -1695,7 +1701,7 @@

    1.1.1 Statistical inference

    xlab("Number of ant species at a site") + ylab("probability mass") + facet_grid(cols = vars(lambda)) -

    +

    Each of these models, with a different value of \(\lambda\), can give us the value of \(X=8\) with a different probability.

    The likelihood function gives us, for each potential value of the parameter, the probability of the data. The figure below shows this for our data with only one observation \(x=8\).

    @@ -1711,7 +1717,7 @@

    1.1.1 Statistical inference

    geom_vline(aes(xintercept = 8), color = "red") + geom_label(aes(label = "ML", x = 8, y = 0.16), fill = "red", alpha = 0.5) -

    +

    With statistical inference we actually want to invert the likelihood \(p(X|\lambda) \rightarrow p(\lambda|X)\). In this way, we hope to find out which model of all potential models is the most meaningful or probable. Both frequentist and Bayesian methods essentially invert \(p(X|\lambda) \rightarrow p(\lambda|X)\), but the method differs.

    @@ -1781,9 +1787,9 @@

    1.2.1 What is MCMC?

    When taking a step uphill, the MCMC always continues. If the step is downhill, the MCMC does not always continue. The probability of going downhill depends on how deep the step down is. Small steps downwards are often taken. If the step downwards is very steep, the MCMC usually will not go into that direction and sample a new step.

    The absolute height of the mountain (posterior likelihood) in itself does not matter. Only the relative differences or shape of the mountain matters. This means that we do not need to calculate the denominator of Bayes’ rule, which is the reason why the algorithm is frequently used in Bayesian statistics.

    -
    - -

    Metropolis rules

    +
    +Metropolis rules +
    Metropolis rules

    A simple MCMC algorithm is the Metropolis algorithm. It is described by these formal rules:

      @@ -1814,7 +1820,7 @@

      1.2.1 What is MCMC?

      ggplot(df_data, aes(x = x, y = y)) +
   geom_point() +
   geom_smooth(method = "lm", se = TRUE)
      -

      +

      # grootte van de plot beperken

      A linear regression with LM yields:

      lm1 <- lm(formula = y ~ x, data = df_data)
@@ -1827,7 +1833,7 @@ 
      1.2.1 What is MCMC?
      
 ## (Intercept) 5.173664 9.558508
 ## x           1.032228 2.688598

      We source some short functions to calculate the (log) likelihood and the prior and to execute the MCMC metropolis algorithm.

      -
      source(file = "./source/mcmc_functions.R")
      +
      source(file = here("static", "code", "brms_modeling", "mcmc_functions.R"))

      For this simple model with a small data set, we can calculate and plot the posterior for a large number of combinations of ‘beta_0’ and ‘beta_1’.

      df <- expand.grid(beta_0 = seq(-1, 15, 0.5), beta_1 = seq(-3, 8, 0.5))
 df$post <- NA
@@ -1840,7 +1846,7 @@ 
      1.2.1 What is MCMC?
      
 
      ggplot(df, aes(x = beta_0, y = beta_1, z = post)) +
   geom_raster(aes(fill = post)) + geom_contour(colour = "black") +
   scale_fill_gradientn(colours = rainbow(n = 4))
      
-
      
+
      
 
      To start the MCMC, we need to choose starting values for \(\beta_0\) and \(\beta_1\):
      
 
      mcmc_small <- MCMC_metro(x = x,         # vector x waarden
                          y = y,         # vector y waarden
@@ -1858,10 +1864,10 @@ 
      1.2.1 What is MCMC?
      
              aes(x = beta_0, y = beta_1, z = NA), colour = "red") +
   coord_cartesian(xlim = c(-1, 14.5), ylim = c(-3, 7)) +
   theme(legend.position = "none")
      
-
      
+
      
 
      The starting value is quite far from the maximum likelihood.
 It takes a while for the MCMC to stabilize in the vicinity of the top of the mountain.
-Therefore, the first part of the MCMC is never used. This is called the ‘burn-in’ or ‘warm-up’.
      
+Therefore, the first part of the MCMC is never used. This is called the ‘warmup’, ‘burn-in’ or ‘tuning’.
      
 
      We now run the same model, but with a much longer MCMC (more iterations):
      
 
      mcmc_l <- MCMC_metro(x = x,
                      y = y,
@@ -1893,12 +1899,12 @@ 
      1.2.1 What is MCMC?
      
 
 gridExtra::grid.arrange(p2, ggplot(), p1, p3, ncol = 2, nrow = 2,
   widths = c(4, 1), heights = c(1, 4))
      
-
      
+
      
 
      The MCMC for each parameter can be displayed in a so called ‘trace plot’ where we visualise the sampled parameter values over time.
      
 
      ann_text <- tibble(param = c("beta_0", "beta_1"),
                    x = 210,
                    y = c(max(mcmc_l$beta_0), max(mcmc_l$beta_1)),
-                   lab = "burn-in")
+                   lab = "warmup")
 mcmc_l %>%
   dplyr::select(iter, beta_0, beta_1) %>%
   pivot_longer(cols = c("beta_0", "beta_1"), names_to = "param") %>%
@@ -1909,7 +1915,7 @@ 
      1.2.1 What is MCMC?
      
   geom_label(data = ann_text, aes(label = lab, x = x, y = y),
              hjust = "left", vjust = "top",
              fill = "gold")
      
-
      
+
      
 
      Summary statistics of beta_0 (intercept)
      
 
      quantile(mcmc_l$beta_0[200:5000], probs = c(0.05, 0.5, 0.95))
      
 
      ##       5%      50%      95% 
@@ -1992,11 +1998,11 @@ 
      1.2.1 What is MCMC?
      
   ggtitle("equal tail credible interval for a bimodal distribution") +
   geom_line(aes(x = x, y = prob))
 p1 
      
-
      
+
      
 
      p2
      
-
      
+
      
 
      p3
      
-
      
+
    @@ -2010,16 +2016,17 @@

    1.3 MCMC software

  • rstan is the basis of brms. You need rstan to use brms because it communicates directly with Stan. A model composed in “pure” Stan code can be fitted with the rstan package in R (so without using brms), but can be quite complex.
  • rstanarm is similar to brms in terms of complexity of the syntax. Moreover, with the rstanarm package you can use some pre-compiled Stan models. This means that fitting your model takes less time. A large part of the process time of brms goes into compiling the Stan program.
    • -
    - -

    Overview of various Stan software (source: https://jtimonen.github.io/posts/post-01/)

    +
    +Overview of various Stan software (source: https://jtimonen.github.io/posts/post-01/) +
    Overview of various Stan software (source: https://jtimonen.github.io/posts/post-01/)
    +

    Although the package brms depends on rstan as a default backend, there is another software tool worth mentioning: cmdstanr. +Compared to rstan, cmdstanr benefits from more frequent updates and some performance advantages, yet it requires extra installation steps. +In the brm(...) function introduced below, you can switch these two with the backend = ... keyword.

    -
    -

    2 Fitting a model with brms

    -
    -

    2.1 Loading the dataset and data exploration

    +
    +

    2 Loading the dataset and data exploration

    We load a dataset on the number of ant species in New England (USA). Type ?ants into the console for more info.

    # Load data
 data(ants)
@@ -2056,7 +2063,7 @@ 
    2.1 Loading the dataset and data
     geom_bar() +
     scale_x_continuous(limits = c(0, NA)) +
     facet_wrap(~habitat)
    -

    +

    # Boxplots of the number of species per habitat
 ants_df %>%
   ggplot(aes(y = sp_rich, x = habitat)) +
@@ -2064,7 +2071,7 @@ 
    2.1 Loading the dataset and data
     geom_line(aes(colour = site, group = site), alpha = 0.5) +
     geom_point(aes(colour = site, group = site)) +
     scale_y_continuous(limits = c(0, NA))
    -

    +

    # Scatter plot of the number of species, related to the latitude per habitat
 ants_df %>%
   ggplot(aes(y = sp_rich, x = latitude)) +
@@ -2072,7 +2079,7 @@ 
    2.1 Loading the dataset and data
     geom_smooth(method = "loess", formula = "y ~ x", colour = "firebrick",
                 level = 0.9) +
     facet_wrap(~habitat)
    -

    +

    # Scatter plot of the number of species, related to the height per habitat
 ants_df %>%
   ggplot(aes(y = sp_rich, x = elevation)) +
@@ -2080,14 +2087,16 @@ 
    2.1 Loading the dataset and data
     geom_smooth(method = "loess", formula = "y ~ x", colour = "firebrick",
                 level = 0.9) +
     facet_wrap(~habitat)
    -

    +

    As an exercise, we will create a model to compare the number of species between both habitats. From the data exploration, we already saw that the number of species seems to be higher in forests and that sites with a higher number in bogs often also have a higher number in forests.

    -
    -

    2.2 Specification of a linear regression

    -
    -

    2.2.1 Model specification

    +
    +

    3 Fitting a model with brms

    +
    +

    3.1 Specification of a linear regression

    +
    +

    3.1.1 Model specification

    We want a model that allows us to investigate the difference in number of ant species between bog and forest habitats. Consider the response variable \(Y\), the number of ants, and \(X_{habitat}\), a dummy variable equal to 0 for bogs and 1 for forests. We assume that \(Y\) follows a Normal distribution (equivalent to a linear regression with categorical variable (= ANOVA)).

    \[ Y \sim N(\mu, \sigma^2) @@ -2100,8 +2109,8 @@

    2.2.1 Model specification

    First of all we decide which MCMC parameters we will use. Type ?brm to see what the default settings are for these parameters. It is recommended to use multiple chains (nchains) and run them in parallel on different cores of your computer (nparallel). Since this is a relatively simple model, we don’t need many iterations and no thinning.

    # Set MCMC parameters
 nchains <- 3 # number of chains
-niter <- 2000 # number of iterations (incl. burn-in, see next)
-burnin <- niter / 4 # number of initial samples to remove (= burn-in)
+niter <- 2000 # number of iterations (incl. warmup, see next)
+warmup <- niter / 4 # number of initial samples to remove (= warmup)
 nparallel <- nchains # number of cores for parallel computing
 thinning <- 1 # thinning factor (here 1 = no thinning)

    The model is fitted using the brm() function. The syntax is very similar to functions used in frequentist statistics such as lm() and glm(). The brm() function contains many arguments (see ?brm). Useful arguments that we do not use here are, for example:

    @@ -2110,22 +2119,22 @@

    2.2.1 Model specification

  • prior to provide prior distributions. If these are not specified, the function uses default (non-informative) priors. The argument sample_prior can be used to sample from the prior so that you can visualise whether the used prior meets your expectations (also useful for prior predictive checks) .
  • file and file_refit to save the model object after it has been fitted. If you run the code again and the model has already been saved, brm() will simply load this model instead of refitting it.
    • -
    # Fit Normal model
+
    # Fit Normal model
 fit_normal1 <- brm(
   formula = sp_rich ~ habitat, # specify the model
   family = gaussian(),         # we use the Normal distribution
   data = ants_df,              # specify data
   chains = nchains,            # MCMC parameters
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
   seed = 123)                  # seed for reproducibility
    
 
    Before we look at the results, we first check whether the model converges well.
    -
    -

    2.2.2 MCMC convergence

    -

    There are several ways to check convergence of the MCMC algorithm for each parameter. The burn-in samples are not taken into account. First and foremost, you have visual controls.

    +
    +

    3.1.2 MCMC convergence

    +

    There are several ways to check convergence of the MCMC algorithm for each parameter. The warmup samples are not taken into account. First and foremost, you have visual controls.

    We can obtain the MCMC samples with the as_draws() functions or visualise them at once via the bayesplot package that is compatible with brmsfit objects.

    # Set colour palette for clear visualisation in bayesplot
 color_scheme_set("mix-blue-red")
    @@ -2142,10 +2151,10 @@

    2.2.2 MCMC convergence

    ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) + geom_line() + facet_wrap(~parameter, nrow = 3, scales = "free")     -

    +

    # Visualise with bayesplot package
 mcmc_trace(fit_normal1, pars = parameters)
    -

    +

    running mean/quantile/… plot

    In addition to a trace plot, we can also look at how certain statistics, such as the average or the quantiles, vary over time (= with increasing number of iterations). You want these statistics to stabilize after a number of iterations because then you know that you have used enough iterations.

    # Code to calculate cumulative quantiles
@@ -2169,21 +2178,21 @@ 
    2.2.2 MCMC convergence
    
   pivot_longer(cols = starts_with(parameters), names_to = "name",
                values_to = "value") %>%
   # split columns names at the last underscore
-  extract(name, into = c("parameter", "statistic"), "(.*)_([^_]+)$") %>%
+  tidyr::extract(name, into = c("parameter", "statistic"), "(.*)_([^_]+)$") %>%
   ggplot(aes(y = value, x = .iteration, colour = factor(.chain))) +
     geom_line(aes(linetype = statistic), linewidth = 0.8) +
     facet_wrap(~parameter, nrow = 3, scales = "free")
    -

    +

    posterior density plot

    The shape of the posterior distributions can also be informative to check whether the chains have converged to the same distributions. If we see a bimodal distribution (a camel’s ridge with two peaks) for example, then there could also be something wrong with the specification of our model.

    # Visualise density plot of the posterior for each of the chains separately
 # via bayesplot package.
 mcmc_dens_overlay(fit_normal1, pars = parameters)
    -

    +

    It is even easier to use the plot() function on the brmsfit object. This immediately shows the trace plots and the posterior density plots side by side for each parameter.

    # Trace plots and posterior density for each parameter.
 plot(fit_normal1)
    -

    +

    autocorrelation plot

    If two variables are correlated, it means that there is a relationship between them. This can be a positive or negative relationship. An example of positively correlated variables is people’s size and weight; the taller a person is, the more they will typically weigh. An example of negatively correlated variables is the amount of sunshine and the moisture content in the top layer of the soil.

    @@ -2200,7 +2209,7 @@

    2.2.2 MCMC convergence

    This will reduce the autocorrelation but also the number of posterior samples we have left.

    # Visualisation of autocorrelation plots via bayesplot package
 mcmc_acf(fit_normal1, pars = parameters)
    -

    +

    In addition to visual checks, there are also diagnostic checks.

    Gelman-Rubin diagnostics

    Also called \(\hat{R}\) (R-hat) or potential scale reduction factor. One way to check whether a chain has converged is to compare its behaviour with other randomly initialized chains. The \(\hat{R}\) statistic takes the ratio of the average variance of samples within each chain to the variance of the pooled samples across all chains. If all chains converge to a common distribution, these ratios will equal 1. If the chains have not converged to a common distribution, the \(\hat{R}\) statistic will be greater than 1. We can extract the \(\hat{R}\) statistics via the rhat() function and then visualise them if desired with the mcmc_rhat() function of the bayesplot package.

    @@ -2208,42 +2217,42 @@

    2.2.2 MCMC convergence

    rhats_fit_normal1 <- rhat(fit_normal1)[parameters] print(rhats_fit_normal1) 
    ##     b_Intercept b_habitatForest           sigma 
-##       1.0003861       1.0004918       0.9998149
    +## 0.999769 1.000314 1.000759     
    # Plot R-hats via bayesplot package
 mcmc_rhat(rhats_fit_normal1) + yaxis_text(hjust = 1)
    -

    +

    effective sample size

    The effective sample size (\(n_{eff}\)) is an estimate of the number of independent draws from the posterior distribution of each parameter. Since the draws within a Markov chain are not independent when there is autocorrelation, the effective sample size is usually smaller than the total sample size, namely the number of iterations (\(N\)). The larger the ratio \(n_{eff}/N\) the better. The bayesplot package provides a neff_ratio() extractor function. The mcmc_neff() function can then be used to plot the ratios.

    # Get ratio of effective sample size
 ratios_fit_normal1 <- neff_ratio(fit_normal1)[parameters]
 print(ratios_fit_normal1)
    ##     b_Intercept b_habitatForest           sigma 
-##       0.7436946       0.6869533       0.7381539
    +## 0.7153819 0.6962991 0.6844347 
    # Plot ratio via bayesplot package
 mcmc_neff(ratios_fit_normal1) + yaxis_text(hjust = 1)
    -

    +

    Other packages that may be useful for checking MCMC convergence are mcmcplots and coda. For these packages you must first convert the MCMC samples from your brms model to an mcmc object.

    # example
 install.packages("mcmcplots")
 mcmcplots::mcmcplot(as.mcmc(fit_normal1, pars = parameters))

    In general, we can conclude that MCMC convergence is good for all parameters in our Normal model. If this were not the case, the MCMC parameters may have to be chosen differently, for example by increasing the number of iterations. Now that we know that the parameters have been estimated correctly, we can check whether the model fits the data well.

    -
    -

    2.2.3 Model fit

    +
    +

    3.1.3 Model fit

    The posterior predictive check (PPC) is a good way to check whether the model fits the data well. The idea behind the PPC is that, if a model fits well, the predicted values based on the model are very similar to the data we used to fit the model. To generate the predictions used for PPCs, we simulate several times from the posterior predictive distribution (the posterior distribution of the response variable). visualisation is possible with the bayesplot package. The thick line shows the data and the thin lines show different simulations based on the model.

    # Visualise model fit via bayesplot package
 pp_check(fit_normal1, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    -

    +

    We see that the observed data and the predictions based on our model do not quite match. Our choice for the Normal distribution may have been wrong. Numbers are discrete, positive values, while the Normal distribution can predict continuous and both positive and negative values. The Poisson distribution is a discrete distribution that is always positive. It is often used to model numbers recorded over a given time interval, distance, area, volume…

    Remark:

    Note that MCMC convergence and model fit are two different things. If MCMC convergence is good, this does not necessarily mean that model fit is also good. MCMC convergence is to check whether the MCMC algorithm has been able to correctly approximate the distributions of the parameters. Model fit checks whether the model we have adopted (normality assumption, linear relationship average…) fits the data well.

    -
    -

    2.3 Specification of a Poisson model

    -
    -

    2.3.1 Model specification

    +
    +

    3.2 Specification of a Poisson model

    +
    +

    3.2.1 Model specification

    We now assume that \(Y\) now follows a Poisson distribution

    \[ Y \sim Pois(\lambda) @@ -2253,39 +2262,39 @@

    2.3.1 Model specification

    \ln(\lambda) = \beta_0 + \beta_1X_{habitat} \]

    So we need to estimate two parameters: \(\beta_0\) and \(\beta_1\) We use the same MCMC parameters as before. The only thing we need to adjust is the choice family = poisson().

    -
    # Fit Poisson model
+
    # Fit Poisson model
 fit_poisson1 <- brm(
   formula = sp_rich ~ habitat, # specify the model
   family = poisson(),          # we use the Poisson distribution
   data = ants_df,              # specify the data
   chains = nchains,            # MCMC parameters
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
   seed = 123)                  # seed for reproducibility
    -
    -

    2.3.2 MCMC convergence

    +
    +

    3.2.2 MCMC convergence

    Convergence looks good.

    # Visualise MCMC convergence of the Poisson model
 plot(fit_poisson1)
    -

    +

    # Extract and visualise R-hat values
 rhats_fit_poisson1 <- rhat(fit_poisson1)[c("b_Intercept", "b_habitatForest")]
 mcmc_rhat(rhats_fit_poisson1) + yaxis_text(hjust = 1)
    -

    +

    -
    -

    2.3.3 Model fit

    +
    +

    3.2.3 Model fit

    We see that all predictions are now positive, but the fit is still not perfect.

    # Visualise model fit via bayesplot package
 pp_check(fit_poisson1, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    -

    +

    -
    -

    2.3.4 Discussion

    +
    +

    3.2.4 Discussion

    How can we further improve model fit?

    We know that each site has been visited twice. Once in bog and once in forest. It is of course possible that the number of ants in the bog and forest of the same site will be correlated (site effect). Indeed, this is something we noticed during data exploration as well. Sites with higher numbers of species in bogs often also have higher numbers in forests and vice versa. We can correct this by adding a random intercept for each site ... + (1|site).

    We assume that the number of species \(Y\) per site \(j\) (\(j = 1, ..., J\)) follows a Poisson distribution

    @@ -2300,13 +2309,13 @@

    2.3.4 Discussion

    \[ b_0 \sim N(0, \sigma_b) \]

    -
    # Fit Poisson model with random intercept per site
+
    # Fit Poisson model with random intercept per site
 fit_poisson2 <- brm(
   formula = sp_rich ~ habitat + (1|site),
   family = poisson(),
   data = ants_df,
   chains = nchains,
-  warmup = burnin, 
+  warmup = warmup, 
   iter = niter,
   cores = nparallel,
   thin = thinning,
@@ -2314,31 +2323,30 @@ 
    2.3.4 Discussion
    
 
    Convergence looks good.
    
 
    # Visualise MCMC convergence
 plot(fit_poisson2)
    
-
    
+
    
 
    # Extract and visualise R-hat values
 parameters2 <- c("b_Intercept", "b_habitatForest", "sd_site__Intercept")
 rhats_fit_poisson2 <- rhat(fit_poisson2)[parameters2]
 mcmc_rhat(rhats_fit_poisson2) + yaxis_text(hjust = 1)
    
-
    
+
    
 
    Model fit looks even better than before.
    
 
    # Visualise model fit of the Poisson model with random intercept via
 # bayesplot package
 pp_check(fit_poisson2, type = "dens_overlay_grouped", ndraws = 100, 
          group = "habitat")
    
-
    
+
    
 
    How can we objectively compare these models?
    -
    -
    -

    3 Compare models

    +
    +

    3.3 Compare models

    Based on the PPCs we can already see which model fits the data best. Furthermore, there are some functions that brms provides to compare different models. With the function add_criterion() you can add model fit criteria to model objects. Type ?add_criterion() to see which ones are available. See also https://mc-stan.org/loo/articles/online-only/faq.html

    -
    -

    3.1 Leave-one-out cross validation

    +
    +

    3.3.1 Leave-one-out cross validation

    Cross-validation (CV) is a family of techniques that attempts to estimate how well a model would predict unknown data through predictions of the model fitted to the known data. You do not necessarily have to collect new data for this. You can split your own data into a test and training dataset. You fit the model to the training dataset and then use that model to estimate how well it can predict the data in the test dataset. With leave-one-out CV (LOOCV) you leave out one observation each time as test dataset and refit the model based on all other observations (= training dataset).

    -
    - -

    Figure LOOCV.

    +
    +Figure LOOCV. +
    Figure LOOCV.
    # Add leave-one-out model fit criterium to model objects
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2353,13 +2361,13 @@ 
    3.1 Leave-one-out cross validatio
                         criterion = "loo")
 print(comp_loo, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_loo se_elpd_loo p_loo    se_p_loo looic   
-## fit_poisson2    0.000     0.000 -106.276    3.849      11.552    1.522  212.552
-## fit_poisson1  -13.403     5.802 -119.679    8.319       3.714    0.888  239.358
-## fit_normal1   -15.647     3.629 -121.923    5.301       3.013    0.842  243.846
+## fit_poisson2    0.000     0.000 -106.077    3.896      11.399    1.569  212.154
+## fit_poisson1  -13.685     5.728 -119.762    8.309       3.813    0.903  239.524
+## fit_normal1   -15.874     3.597 -121.951    5.329       3.007    0.841  243.902
 ##              se_looic
-## fit_poisson2    7.698
-## fit_poisson1   16.637
-## fit_normal1    10.601
    +## fit_poisson2 7.792 +## fit_poisson1 16.618 +## fit_normal1 10.659

    If you assume that this difference is normally distributed, you can calculate a confidence interval based on the uncertainty. We see that the second Poisson model scores best and that the other models score significantly lower.

    # Calculate confidence interval to compare models with loocv
 comp_loo %>%
@@ -2369,15 +2377,15 @@ 
    3.1 Leave-one-out cross validatio
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    ##              elpd_diff  se_diff   ll_diff   ul_diff
 ## fit_poisson2   0.00000 0.000000   0.00000  0.000000
-## fit_poisson1 -13.40293 5.802467 -22.94714 -3.858719
-## fit_normal1  -15.64659 3.628604 -21.61511 -9.678068
    +## fit_poisson1 -13.68487 5.727656 -23.10603 -4.263713 +## fit_normal1 -15.87376 3.596654 -21.78973 -9.957793
    -
    -

    3.2 K-fold cross-validation

    +
    +

    3.3.2 K-fold cross-validation

    With K-fold cross-validation, the data is split into \(K\) groups. We will use \(K = 10\) groups (= folds) here. So instead of leaving out a single observation each time, as with leave-one-out CV, we will leave out one \(10^{th}\) of the data here. Via the arguments folds = "stratified" and group = "habitat" we ensure that the relative frequencies of habitat are preserved for each group. This technique will therefore be less precise than the previous one, but will be faster to calculate if you work with a lot of data.

    -
    - -

    Figure K-fold CV.

    +
    +Figure K-fold CV. +
    Figure K-fold CV.
    # Add K-fold model fit criterium to model objects
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2401,9 +2409,9 @@ 
    3.2 K-fold cross-validation
    
                           criterion = "kfold")
 print(comp_kfold, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_kfold se_elpd_kfold p_kfold  se_p_kfold
-## fit_poisson2    0.000     0.000 -106.773      4.145        12.048    1.944  
-## fit_poisson1  -12.290     5.195 -119.063      8.156         3.097    0.939  
-## fit_normal1   -14.619     3.258 -121.392      5.172         2.482    0.796
    +## fit_poisson2 0.000 0.000 -106.516 3.890 11.837 1.561 +## fit_poisson1 -12.646 5.808 -119.162 8.227 3.212 1.025 +## fit_normal1 -15.101 3.616 -121.617 5.148 2.672 0.717

    The results are the same, but the differences are slightly smaller than before.

    # Calculate confidence interval to compare models with K-fold CV
 comp_kfold %>%
@@ -2413,11 +2421,11 @@ 
    3.2 K-fold cross-validation
    
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    ##              elpd_diff  se_diff   ll_diff   ul_diff
 ## fit_poisson2   0.00000 0.000000   0.00000  0.000000
-## fit_poisson1 -12.28990 5.195211 -20.83526 -3.744537
-## fit_normal1  -14.61938 3.258315 -19.97883 -9.259928
    +## fit_poisson1 -12.64595 5.807731 -22.19882 -3.093082 +## fit_normal1 -15.10095 3.616047 -21.04882 -9.153085
    -
    -

    3.3 WAIC

    +
    +

    3.3.3 WAIC

    The Widely Applicable Information Criterion (WAIC) does not use cross-validation but is a computational way to estimate the ELPD. How this happens exactly is beyond the purpose of this tutorial. It is yet another measure to apply model selection.

    # Add WAIC model fit criterion to model objects
 fit_normal1 <- add_criterion(fit_normal1,
@@ -2432,29 +2440,30 @@ 
    3.3 WAIC
    
                          criterion = "waic")
 print(comp_waic, simplify = FALSE, digits = 3)
    ##              elpd_diff se_diff  elpd_waic se_elpd_waic p_waic   se_p_waic
-## fit_poisson2    0.000     0.000 -105.062     3.683       10.338    1.348 
-## fit_poisson1  -14.594     5.858 -119.656     8.312        3.691    0.882 
-## fit_normal1   -16.833     3.624 -121.895     5.289        2.985    0.830 
+## fit_poisson2    0.000     0.000 -104.809     3.678       10.131    1.344 
+## fit_poisson1  -14.934     5.851 -119.743     8.305        3.794    0.899 
+## fit_normal1   -17.104     3.620 -121.913     5.313        2.969    0.824 
 ##              waic     se_waic 
-## fit_poisson2  210.124    7.367
-## fit_poisson1  239.313   16.624
-## fit_normal1   243.790   10.578
    +## fit_poisson2 209.618 7.355 +## fit_poisson1 239.486 16.609 +## fit_normal1 243.827 10.626 
    # Calculate confidence interval to compare models with WAIC
 comp_waic %>%
   as.data.frame() %>%
   select(waic, elpd_diff, se_diff) %>%
   mutate(ll_diff = elpd_diff  + qnorm(0.05) * se_diff,
          ul_diff = elpd_diff  + qnorm(0.95) * se_diff)
    -
    ##                  waic elpd_diff  se_diff   ll_diff    ul_diff
-## fit_poisson2 210.1239   0.00000 0.000000   0.00000   0.000000
-## fit_poisson1 239.3128 -14.59444 5.857609 -24.22935  -4.959532
-## fit_normal1  243.7901 -16.83311 3.624302 -22.79456 -10.871665
    +
    ##                  waic elpd_diff  se_diff   ll_diff   ul_diff
+## fit_poisson2 209.6176   0.00000 0.000000   0.00000   0.00000
+## fit_poisson1 239.4857 -14.93407 5.850998 -24.55810  -5.31003
+## fit_normal1  243.8265 -17.10447 3.620256 -23.05926 -11.14968
    -
    -

    3.4 Conclusion

    +
    +

    3.3.4 Conclusion

    Both based on the PPC and the comparisons with different model selection criteria, we can conclude that the second Poisson model with random intercepts fits the data best. In principle, we could have expected this based on our own intuition and the design of the study, i.e. the use of the Poisson distribution to model numbers and the use of random intercepts to control for a hierarchical design (habitats nested within sites).

    +

    4 Final model results

    When we look at the model fit object, we see results that are similar to results we see when we fit a frequentist model. On the one hand we get an estimate of all parameters with their uncertainty, but on the other hand we see that this is clearly the output of a Bayesian model. We get information about the parameters we used for the MCMC algorithm, we get a 95% credible interval (CI) instead of a confidence interval and we also get the \(\hat{R}\) value for each parameter as discussed earlier.

    @@ -2467,15 +2476,15 @@

    4 Final model results

    ## Draws: 3 chains, each with iter = 2000; warmup = 500; thin = 1; ## total post-warmup draws = 4500 ## -## Group-Level Effects: +## Multilevel Hyperparameters: ## ~site (Number of levels: 22) ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS -## sd(Intercept) 0.38 0.10 0.22 0.60 1.00 1748 2738 +## sd(Intercept) 0.38 0.09 0.23 0.60 1.00 1733 2561 ## -## Population-Level Effects: +## Regression Coefficients: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS -## Intercept 1.51 0.13 1.25 1.77 1.00 2873 2723 -## habitatForest 0.64 0.12 0.41 0.88 1.00 7102 3251 +## Intercept 1.51 0.13 1.24 1.75 1.00 2452 2883 +## habitatForest 0.64 0.12 0.41 0.87 1.00 5062 3193 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential @@ -2498,8 +2507,8 @@

    4 Final model results

    ## # A tibble: 2 × 9
 ##   .variable         min  q_05  q_20  mean median  q_80  q_95   max
 ##   <chr>           <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
-## 1 b_Intercept     1.04  1.29  1.40  1.51   1.52  1.63  1.74   1.97
-## 2 b_habitatForest 0.268 0.432 0.537 0.638  0.639 0.742 0.831  1.00
    +## 1 b_Intercept 1.06 1.29 1.41 1.52 1.52 1.62 1.72 1.95 +## 2 b_habitatForest 0.275 0.454 0.544 0.641 0.639 0.742 0.827 1.07

    Useful functions of the tidybayes package are also median_qi(), mean_qi() … after gather_draws() which you can use instead of group_by() and summarise() .

    We would now like to visualise the estimated number of species per habitat type with associated uncertainty. With the function spread_draws() you can take a certain number of samples from the posterior distribution and convert them into a wide format table. The average number of species in bogs according to our model is \(\exp(\beta_0)\) and in forests \(\exp(\beta_0+\beta_1)\). We show the posterior distributions with the posterior median and 60 and 90% credible intervals.

    fit_poisson2 %>%
@@ -2514,7 +2523,7 @@ 
    4 Final model results
    
   ggplot(aes(y = sp_rich, x = habitat)) +
     stat_eye(point_interval = "median_qi", .width = c(0.6, 0.9)) +
     scale_y_continuous(limits = c(0, NA))
    -

    +

    In addition to stat_eye() you will find here some nice ways to visualise posterior distributions .

    We see a clear difference in the number of species between the two habitats. Is there a significant difference between the number of species in bogs and forests? We test the hypothesis that numbers are equal in bogs and forests.

    \[ @@ -2530,7 +2539,7 @@

    4 Final model results

    hyp 
    ## Hypothesis Tests for class b:
 ##            Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
-## 1 (habitatForest) = 0     0.64      0.12     0.44     0.84         NA        NA
+## 1 (habitatForest) = 0     0.64      0.12     0.45     0.83         NA        NA
 ##   Star
 ## 1    *
 ## ---
@@ -2540,7 +2549,7 @@ 
    4 Final model results
    
 ## Posterior probabilities of point hypotheses assume equal prior probabilities.
    # Plot posterior distribution hypothesis
 plot(hyp)
    -

    +

    We can conclude that there is a significant difference since 0 is not included in the 90% credible interval.

    Finally, we visualise the random effects of the sites. We sort them from high to low species richness.

    # Take the mean of SD of random effects
@@ -2561,7 +2570,7 @@ 
    4 Final model results
    
                color = "darkgrey", linetype = 2) +
     stat_halfeye(point_interval = "median_qi", .width = 0.9, size = 2/3,
                  fill = "cornflowerblue")
    -

    +

    5 Comparison with frequentist statistics

    @@ -2580,14 +2589,14 @@

    5 Comparison with frequentist sta ## Draws: 3 chains, each with iter = 2000; warmup = 500; thin = 1; ## total post-warmup draws = 4500 ## -## Population-Level Effects: +## Regression Coefficients: ## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS -## Intercept 4.84 0.80 3.53 6.16 1.00 4194 3347 -## habitatForest 4.33 1.14 2.44 6.20 1.00 4113 3091 +## Intercept 4.80 0.81 3.46 6.11 1.00 4196 3219 +## habitatForest 4.33 1.14 2.42 6.15 1.00 4502 3133 ## -## Family Specific Parameters: +## Further Distributional Parameters: ## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS -## sigma 3.73 0.41 3.13 4.44 1.00 3925 3322 +## sigma 3.73 0.41 3.12 4.46 1.00 3557 3080 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential @@ -2610,7 +2619,141 @@

    5 Comparison with frequentist sta ## sample estimates: ## mean in group Bog mean in group Forest ## 4.863636 9.181818 -

    We see that this indeed produces almost exactly the same results. Our Bayesian model estimates that on average 4.326 more ant species occur in forests than in bogs (90% credible interval: 2.439 to 6.203). The t-test estimates that on average 4.318 more ant species occur in forests than in bogs (90% confidence interval: 2.445 to 6.192).

    +

    We see that this indeed produces almost exactly the same results. Our Bayesian model estimates that on average 4.329 more ant species occur in forests than in bogs (90% credible interval: 2.423 to 6.149). The t-test estimates that on average 4.318 more ant species occur in forests than in bogs (90% confidence interval: 2.445 to 6.192).

    +

    +
    +

    6 Deep Dive: rstan

    +
    +

    6.1 Stan: What? Why?!

    +

    The brms package is a convenience wrapper for the rstan package, which in turn ports stan functionality to R. +Stan is a modeling framework written in the C programming language, which implements many probabilistic (“Bayesian”) modeling tools. +More info can be found on the Stan website.

    +

    The advantage of brms is usability: many functions work out-of-the-box, with reasonable default values, and a syntax that is similar to what many statisticians are habituated to. +However, the relative ease-of-use comes at the cost of flexibility, and do some degree, readability.

    +

    In contrast, Stan and rstan lean more to the mathematical formulation of models. +Every aspect of the model has to be explicitly set, which can be an advantage (e.g. if you face non-standard use cases), or disadvantage (e.g. if you specify models in non-optimal ways).

    +

    To briefly give an impression, we will build the same models as above, using the Stan framework.

    +
    library("rstan")
+conflicted::conflicts_prefer(rstan::extract)
+conflicted::conflicts_prefer(brms::loo)
    +
    +
    +

    6.2 Model Definition

    +

    RMarkdown can handle stan code chunks, though more general model definition is outsourced to a separate “*.stan” file (e.g. here). +Alternatively, you can define your model in a big text block, as shown below. +The simple poisson model resembles one of the stan-dard examples, which you can refer to for all further details and more.

    +
    stan_poisson_model_code = '
+data {
+  int<lower=1> N;                   // sample size
+  vector[N] habitat_obs;            // habitat
+  array[N] int<lower=0> sp_rich;    // outcome variable: number of ants
+}
+
+parameters {
+  real intercept_rich;              // intercept
+  real habitat_effect;              // slope
+}
+
+model {
+  // priors
+  intercept_rich ~ normal(0, 1);
+  habitat_effect ~ normal(0, 1);
+
+  // "posterior", "likelihood", ... give it a name!
+  sp_rich ~ poisson_log(intercept_rich + habitat_effect * habitat_obs);
+
+}
+
+'
    +

    Take this as a “look behind the scenes”! +You have to explicitly define the model structure, priors, even the data types of input variables. +Yet note how even Stan is not without convenience functions: we can use the poisson_log posterior to model log rates.

    +

    When working outside RStudio/RMarkdown, you might prefer loading the model from a file:

    +
    stan_poisson_model <- stan_model(
+  # file = here("code", "brms_modeling", "poisson_model.stan"),
+  model_code = stan_poisson_model_code, 
+  model_name = "stan poisson model"
+  )
    +
    +
    +

    6.3 Sampling

    +

    Sampling does pretty much the same as above, since at the core, brms is just stan.

    +
    stan_poisson_fit <- sampling(
+  stan_poisson_model,
+  list(
+    N = nrow(ants_df),
+    habitat_obs = as.integer(ants_df$habitat)-1,
+    sp_rich = ants_df$sp_rich
+    ),
+  iter = niter,
+  chains = nchains,
+  cores = nparallel
+  )
+
+stan_poisson_fit
    +

    A convenient way of quickly inspecting model fits is shinystan. +It opens a browser with shiny plots.

    +
    library("shinystan")
+launch_shinystan(stan_poisson_fit)
    +

    Behold: model outcome is exactly as above. +In the present case, the benefit from turning to stan is very limited. +In other cases, it might pay off.

    +

    Know that Stan is there for you, do not hesitate to turn to its extensive documentation, and do not fear to give it a try!

    +
    +
    +

    6.4 Homework: Hierarchical Model

    +

    To take your modeling skills even further, you may implement and sample the “random intercept” model. +In “Bayesian” terms, the general terminology is “hierarchical” model.

    +

    Below is the code which considers site-specific intercepts. +It works slightly different from the +(1|site) approach above: the upper one is an “offset” parametrization, whereas this time we use a hyperprior. +These are two common strategies often worth interchanging, with marginal or substantial effects on sampler performance.

    +
    data {
+  int<lower=1> N;                   // sample size
+  int<lower=1> L;                   // number of sites
+  vector[N] habitat_obs;            // habitat
+  array[N] int<lower=1, upper=L> site_obs; // site
+  array[N] int<lower=0> sp_rich;    // outcome variable: number of ants
+}
+
+parameters {
+  real intercept_rich;              // "global" intercept
+  vector[L] site_effect;            // the sr intercept at each site
+  real<lower=0> site_variation;     // variation on site level
+  real habitat_effect;              // habitat slope
+}
+
+model {
+  // priors
+  intercept_rich ~ normal(0, 1);
+  site_variation ~ cauchy(0, 1);
+  habitat_effect ~ normal(0, 1);
+
+  // site-wise intercept
+  for (i in 1:L) {
+    site_effect[i] ~ normal(intercept_rich, site_variation);
+  }
+
+  // "posterior", "likelihood", ... give it a name!
+  sp_rich ~ poisson_log(site_effect[site_obs] + habitat_effect * habitat_obs);
+
+}
    +
    stan_poisson_fit <- sampling(
+  stan_poisson_model,
+  list(
+    N = nrow(ants_df),
+    L = length(levels(ants_df$site)),
+    habitat_obs = as.integer(ants_df$habitat)-1,
+    site_obs = as.integer(ants_df$site),
+    sp_rich = ants_df$sp_rich
+    ),
+  iter = niter,
+  chains = nchains,
+  cores = nparallel
+  )
+
+stan_poisson_fit
    +

    With Stan, po(i)ssibilities are almost endless - don’t get lost in model building!

    +

    References