Merge pull request #346 from inbo/rstan

BRMS tutorial: add `rstan` example

Author: falkmielke

.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

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)) +

content/tutorials/r_brms/brms_eng/workshop_1_mcmc_en_brms_eng.Rmd

@@ -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)) +