Some of these things are somewhat like the others.
Today’s goal is to look at a couple of different model structures that we saw yesterday.
library(tidyverse)
library(cmdstanr)This is cmdstanr version 0.7.1- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr- CmdStan path: /home/andrew/software/cmdstan- CmdStan version: 2.34.1library(tidybayes)
library(palmerpenguins)Are populations of penguins on different islands different in their body mass?
The Palmer penguins are found on three different islands. Let’s look at the distribution of body mass of each species on each island.
penguin_mass_island <- penguins |>
select(species, island, body_mass_g) |>
drop_na(body_mass_g) |>
unite(sp_island, species, island) |>
## center mass and change the units
mutate(mass_kg = (body_mass_g)/1000)penguin_mass_island |>
ggplot(aes(y = sp_island,
x = mass_kg,
colour = sp_island)) +
geom_jitter(alpha = 0.8, height = 0.1, width = 0) +
scale_color_brewer(palette = "Dark2")
Are the sample sizes equal among the species-island combinations?
penguin_mass_island |>
count(sp_island) |>
knitr::kable()| sp_island | n |
|---|---|
| Adelie_Biscoe | 44 |
| Adelie_Dream | 56 |
| Adelie_Torgersen | 51 |
| Chinstrap_Dream | 68 |
| Gentoo_Biscoe | 123 |
We’ll begin by fitting a model that assumes that body size for each of these five groups is completely independent:
\[ \begin{align} \text{Body mass}_i &\sim \text{Normal}(\mu_i, \sigma_{\text{obs}}) \\ \mu_i &= \bar\beta + \beta_{\text{group}[i]} \\ \bar\beta &\sim \text{Normal}(5, 2) \\ \beta_{\text{group}} &\sim \text{Normal}(0, 1) \\ \sigma_{\text{obs}} &\sim \text{Exponential}(.5) \end{align} \]
Here’s a little trick to get group indexes (numbers) from a character vector:
group_names <- unique(penguin_mass_island$sp_island)
group_numbers <- seq_along(group_names)
names(group_numbers) <- group_names
group_numbersAdelie_Torgersen Adelie_Biscoe Adelie_Dream Gentoo_Biscoe
1 2 3 4
Chinstrap_Dream
5 penguin_groupid <- penguin_mass_island |>
mutate(group_id = group_numbers[sp_island])
penguin_groupid# A tibble: 342 × 4
sp_island body_mass_g mass_kg group_id
<chr> <int> <dbl> <int>
1 Adelie_Torgersen 3750 3.75 1
2 Adelie_Torgersen 3800 3.8 1
3 Adelie_Torgersen 3250 3.25 1
4 Adelie_Torgersen 3450 3.45 1
5 Adelie_Torgersen 3650 3.65 1
6 Adelie_Torgersen 3625 3.62 1
7 Adelie_Torgersen 4675 4.68 1
8 Adelie_Torgersen 3475 3.48 1
9 Adelie_Torgersen 4250 4.25 1
10 Adelie_Torgersen 3300 3.3 1
# ℹ 332 more rowsAs you can see, we’re set up now with the names and the indexes we need.
Now we can simulate data and plot it:
ngroup <- length(group_numbers)
overall_mean <- rnorm(1, mean = 5, sd = 2)
group_diffs <- rnorm(n = ngroup, mean = 0, sd = 1)
sigma_obs <- rexp(1, .5)
penguin_pred_obs <- penguin_groupid |>
mutate(fake_mass_avg = overall_mean + group_diffs[group_id],
fake_mass_obs = rnorm(length(fake_mass_avg),
mean = fake_mass_avg,
sd = sigma_obs))
penguin_pred_obs |>
ggplot(aes(y = sp_island,
x = fake_mass_obs,
colour = sp_island)) +
geom_jitter(alpha = 0.8, height = 0.1, width = 0) +
scale_color_brewer(palette = "Dark2")
Run the above code a few times! if you want, try different prior values.
fixed_groups <- cmdstan_model(stan_file = "topics/03_one_random_effect/fixed_groups.stan")
fixed_groupsdata {
int<lower=0> N;
vector[N] y;
int<lower=0> Ngroup;
array[N] int<lower=0, upper=Ngroup> group_id;
}
parameters {
real b_avg;
vector[Ngroup] b_group;
real<lower=0> sigma;
}
model {
y ~ normal(b_avg + b_group[group_id], sigma);
b_group ~ std_normal();
b_avg ~ normal(5, 2);
sigma ~ exponential(.5);
}
generated quantities {
vector[Ngroup] group_averages;
for (k in 1:Ngroup){
group_averages[k] = b_avg + b_group[k];
}
// predict making one new observation per group
vector[Ngroup] one_obs_per_group;
for (k in 1:Ngroup) {
one_obs_per_group[k] = normal_rng(group_averages[k], sigma);
}
}peng_group_list <- with(penguin_groupid,
list(
N = length(mass_kg),
y = mass_kg,
Ngroup = max(group_id),
group_id = group_id
))
fixed_groups_samples <- fixed_groups$sample(
data = peng_group_list,
refresh = 0,
parallel_chains = 4
)Running MCMC with 4 parallel chains...Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpH4tcI4/model-4700c1ddd09d5.stan', line 13, column 2 to column 47)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 1 finished in 0.4 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.4 seconds.Let’s begin by plotting the averages for each group.
fixed_groups_samples |>
gather_rvars(group_averages[group_id]) |>
mutate(sp_island = names(group_numbers)[group_id]) |>
ggplot(aes(y = sp_island, dist = .value)) +
stat_pointinterval() +
geom_jitter(data = penguin_mass_island,
aes(y = sp_island,
x = mass_kg,
colour = sp_island),
pch = 21, inherit.aes = FALSE,
alpha = 0.8, height = 0.1, width = 0) +
scale_colour_brewer(palette = "Dark2")
Some things to notice about the code above:
group_numbers to re-create the column sp_island. This is my technique for making sure I always use the correct label, but you can do this any way you want.tidybayes::stat_pointinterval() to summarize the posterior distribution.penguin_mass_island) with geom_jitter(). We’re adding noise vertically to make the visualization better, but not adding any horizontal noise.Repeat the exercise above using the value of one_obs_per_group. Why are the results different? What additional error is included in these predictions?
fixed_groups_samples |>
tidybayes::gather_rvars(one_obs_per_group[group_id]) |>
mutate(sp_island = group_names[group_id]) |>
ggplot(aes(y = sp_island,
dist = .value,
colour = sp_island)) +
stat_pointinterval(colour = "black") +
geom_jitter(
aes(y = sp_island,
x = mass_kg,
colour = sp_island),
inherit.aes = FALSE,
alpha = .2, data = penguin_groupid, height = .2, width = 0) +
scale_colour_brewer(palette = "Dark2")
\[ \begin{align} \text{Body mass}_i &\sim \text{Normal}(\mu_i, \sigma_{\text{obs}}) \\ \mu_i &= \bar\beta + \beta_{\text{group}[i]} \\ \bar\beta &\sim \text{Normal}(5, 2) \\ \beta_{\text{group}} &\sim \text{Normal}(0, 1) \\ \sigma_{\text{obs}} &\sim \text{Exponential}(.5) \end{align} \]
\[ \begin{align} \text{Body mass}_i &\sim \text{Normal}(\mu_i, \sigma_{\text{obs}}) \\ \mu_i &= \bar\beta + \beta_{\text{group}[i]} \\ \bar\beta &\sim \text{Normal}(5, 2) \\ \beta_{\text{group}} &\sim \text{Normal}(0, \sigma_{\text{sp}}) \\ \sigma_{\text{obs}} &\sim \text{Exponential}(.5) \\ \sigma_{\text{sp}} &\sim \text{Exponential}(1) \end{align} \]
Simulate from the model above. Base your approach on the code for simulation the non-hierarchical version. Remember to simulate one additional number: the standard deviation of group differences
ngroup <- length(group_numbers)
overall_mean <- rnorm(1, mean = 5, sd = 2)
sigma_group <- rexp(1, .1)
group_diffs <- rnorm(n = ngroup, mean = 0, sd = sigma_group)
sigma_obs <- rexp(1, .5)
penguin_pred_obs <- penguin_groupid |>
mutate(fake_mass_avg = overall_mean + group_diffs[group_id],
fake_mass_obs = rnorm(length(fake_mass_avg),
mean = fake_mass_avg,
sd = sigma_obs))
penguin_pred_obs |>
ggplot(aes(y = sp_island,
x = fake_mass_obs,
colour = sp_island)) +
geom_jitter(alpha = 0.8, height = 0.1, width = 0) +
scale_color_brewer(palette = "Dark2")
Below I’m comparing the two Stan programs side-by-side. Compare them to the models above!
hierarchical_groups <- cmdstan_model(stan_file = "topics/03_one_random_effect/hierarchical_groups.stan")fixed_groupsdata {
int<lower=0> N;
vector[N] y;
int<lower=0> Ngroup;
array[N] int<lower=0, upper=Ngroup> group_id;
}
parameters {
real b_avg;
vector[Ngroup] b_group;
real<lower=0> sigma;
}
model {
y ~ normal(b_avg + b_group[group_id], sigma);
b_group ~ std_normal();
b_avg ~ normal(5, 2);
sigma ~ exponential(.5);
}
generated quantities {
vector[Ngroup] group_averages;
for (k in 1:Ngroup){
group_averages[k] = b_avg + b_group[k];
}
// predict making one new observation per group
vector[Ngroup] one_obs_per_group;
for (k in 1:Ngroup) {
one_obs_per_group[k] = normal_rng(group_averages[k], sigma);
}
}hierarchical_groupsdata {
int<lower=0> N;
vector[N] y;
int<lower=0> Ngroup;
array[N] int<lower=0, upper=Ngroup> group_id;
}
parameters {
real b_avg;
vector[Ngroup] b_group;
real<lower=0> sigma_obs;
real<lower=0> sigma_grp;
}
model {
y ~ normal(b_avg + b_group[group_id], sigma_obs);
b_group ~ normal(0, sigma_grp);
b_avg ~ normal(5, 2);
sigma_obs ~ exponential(.5);
sigma_grp ~ exponential(1);
}
generated quantities {
vector[Ngroup] group_averages;
for (k in 1:Ngroup){
group_averages[k] = b_avg + b_group[k];
}
// predict making one new observation per group
vector[Ngroup] one_obs_per_group;
for (k in 1:Ngroup) {
one_obs_per_group[k] = normal_rng(group_averages[k], sigma_obs);
}
// difference for a new group
real new_b_group = normal_rng(0, sigma_grp);
// observations from that new group
real one_obs_new_group = normal_rng(b_avg + new_b_group, sigma_obs);
}hierarchical_groups_samples <- hierarchical_groups$sample(
data = peng_group_list, refresh = 0, parallel_chains = 4)Running MCMC with 4 parallel chains...Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpWhWKuf/model-473001618f1b5.stan', line 15, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 finished in 0.4 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.4 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.4 seconds.
Total execution time: 0.4 seconds.hierarchical_groups_samples variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ 88.39 88.71 2.14 1.99 84.52 91.22 1.00 950 1692
b_avg 4.02 4.01 0.37 0.32 3.45 4.65 1.01 519 643
b_group[1] -0.31 -0.29 0.38 0.33 -0.95 0.27 1.01 539 642
b_group[2] -0.31 -0.29 0.38 0.32 -0.96 0.28 1.01 535 632
b_group[3] -0.33 -0.32 0.38 0.33 -0.97 0.24 1.01 523 652
b_group[4] 1.05 1.06 0.37 0.32 0.43 1.62 1.01 530 638
b_group[5] -0.29 -0.28 0.38 0.32 -0.92 0.29 1.01 529 660
sigma_obs 0.47 0.46 0.02 0.02 0.44 0.50 1.00 1808 1534
sigma_grp 0.79 0.69 0.37 0.25 0.41 1.43 1.00 1167 977
group_averages[1] 3.71 3.71 0.07 0.07 3.60 3.81 1.00 4589 2590
# showing 10 of 21 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)hierarchical_groups_samples |>
tidybayes::gather_rvars(b_group[group_id],
new_b_group) |>
mutate(sp_island = group_names[group_id],
sp_island = if_else(is.na(sp_island),
true = "New Group",
false = sp_island)) |>
ggplot(aes(y = sp_island,
dist = .value,
colour = sp_island)) +
stat_pointinterval()
hierarchical_groups_samples |>
tidybayes::gather_rvars(one_obs_per_group[group_id],
one_obs_new_group) |>
mutate(sp_island = group_names[group_id],
sp_island = if_else(is.na(sp_island),
true = "New Group",
false = sp_island)) |>
ggplot(aes(y = sp_island,
dist = .value,
colour = sp_island)) +
stat_pointinterval() +
geom_point(aes(y = sp_island,
x = mass_kg,
colour = sp_island),
inherit.aes = FALSE,
alpha = .2, data = penguin_groupid)
year as a part of the group-creating factor (i.e. in the call to unite() above). What changes?generated quantities block to simulate a fake observation for EVERY row of the dataset. This opens the possibility of using bayesplot to make predictions. Look back at the code from Day 1 and create a posterior predictive check for both models. (e.g. using ppc_dens_overlay)sex as a grouping factor, but sex has missing values in it! Why is this a problem for this kind of model? What would it take to address that? (Discussion only; missing values are unfortunately outside the scope of the class!)Let’s write a model to answer the question:
How does the total abundance of the mite community change as water content increases?
Here’s a partially complete model for species richness over time
\[ \begin{align} \text{S}_i &\sim \text{Poisson}(e^a) \\ a &= \bar\beta + \beta_{\text{water}} \cdot \text{water}_i \\ \bar\beta &\sim \text{Normal}(?, ?) \\ \beta_{\text{water}} &\sim \text{Normal}(?, ?) \\ \end{align} \]
Simulate from this model, and look at your simulations to decide on a reasonable prior for the data.
n <- 30
water <- seq(from = -5, to = 5, length.out = n)
b0 <- rnorm(1, mean = log(17), sd = .3)
b1 <- rnorm(1, mean = 0, sd = .2)
S <- rpois(n, lambda = exp(b0 + b1*water))
plot(water, S)
First we need to load and prepare the data:
data(mite, package = "vegan")
data("mite.env", package = "vegan")
# combine data and environment
mite_data_long <- mite |>
tibble::rownames_to_column(var = "site_id") |>
bind_cols(mite.env) |>
pivot_longer(Brachy:Trimalc2,
names_to = "spp", values_to = "abd")First let’s transform the mite dataset into a dataframe of total community abundance (N) per site. We’ll also standardize the water content while we’re at it:
mite_community_abd <- mite_data_long |>
group_by(site_id, WatrCont) |>
summarize(N = sum(abd)) |>
ungroup() |>
mutate(water_c = (WatrCont - mean(WatrCont))/100)`summarise()` has grouped output by 'site_id'. You can override using the
`.groups` argument.knitr::kable(head(mite_community_abd))| site_id | WatrCont | N | water_c |
|---|---|---|---|
| 1 | 350.15 | 140 | -0.6048571 |
| 10 | 220.73 | 166 | -1.8990571 |
| 11 | 134.13 | 216 | -2.7650571 |
| 12 | 405.91 | 213 | -0.0472571 |
| 13 | 243.70 | 177 | -1.6693571 |
| 14 | 239.51 | 269 | -1.7112571 |
We get a nice histogram of community abundance, and a clear negative relationship with water volume:
mite_community_abd |>
ggplot(aes(x = N)) +
geom_histogram()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.mite_community_abd |>
ggplot(aes(x = water_c, y = N)) +
geom_point()
poisson_regression <- cmdstan_model(stan_file = "topics/03_one_random_effect/poisson_regression.stan")water_for_pred <- seq(from = -3, to = 4.5, length.out = 15)
abd_data_list <- list(N = length(mite_community_abd$N),
water = mite_community_abd$water_c,
y = mite_community_abd$N,
Npred = 15,
water_pred = water_for_pred)
poisson_regression_sample <- poisson_regression$sample(
data = abd_data_list,
refresh = 0)Running MCMC with 4 sequential chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.0 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.5 seconds.poisson_regression_sample |>
tidybayes::gather_rvars(line_obs[i]) |>
mutate(water = water_for_pred) |>
ggplot(aes(x = water, dist = .value)) +
stat_lineribbon() +
geom_point(aes(x = water_c, y = N),
data = mite_community_abd,
inherit.aes = FALSE) +
scale_fill_brewer(palette = "Greens")
fake_obs_S <- poisson_regression_sample$draws(variables = "fake_obs")
fake_obs_S_matrix <- posterior::as_draws_matrix(fake_obs_S)
bayesplot::ppc_dens_overlay(y = mite_community_abd$N,
yrep = head(fake_obs_S_matrix, 50))
Remember, on the left we are plotting the Prediction interval here: it’s showing the distribution of probable observations according to the model. Notice that the the model predicts much narrower variation than we really find!
On the right hand side we have the posterior predictive check, which once again shows that the model is overconfident and predicts a range of observations that are far too narrow.
Discuss with your neighbours: would you trust this model? Would you publish it? The technical name for this phenomenon is “overdisperson”. Have you checked for this in previous count models you’ve done?
Add a random effect for every individual observation in the model. Begin by writing the mathematical notation for this new model!
fit the model and re-create the two figures above. What do you notice?
\[ \begin{align} \text{S}_i &\sim \text{Poisson}(e^a) \\ a &= \bar\beta + \beta_{\text{water}} \cdot \text{water}_i + \text{site}_i\\ \bar\beta &\sim \text{Normal}(?, ?) \\ \beta_{\text{water}} &\sim \text{Normal}(?, ?) \\ \text{site} &\sim \text{Normal}(?, \sigma) \\ \sigma &\sim \text{Exponential}(?) \end{align} \]
poisson_regression_overdisp <- cmdstan_model(stan_file = "topics/03_one_random_effect/poisson_regression_overdisp.stan")
poisson_regression_overdisp_sample <- poisson_regression_overdisp$sample(
data = abd_data_list, refresh = 0)Running MCMC with 4 sequential chains...Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpIhvutV/model-471231bcd7800.stan', line 18, column 2 to column 42)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpIhvutV/model-471231bcd7800.stan', line 18, column 2 to column 42)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 finished in 0.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 1.3 seconds.fake_obs_S <- poisson_regression_overdisp_sample$draws(variables = "fake_obs")
fake_obs_S_matrix <- posterior::as_draws_matrix(fake_obs_S)
bayesplot::ppc_dens_overlay(y = mite_community_abd$N,
yrep = head(fake_obs_S_matrix, 50))
poisson_regression_overdisp_sample |>
tidybayes::gather_rvars(line_obs[i]) |>
mutate(water = water_for_pred) |>
ggplot(aes(x = water, dist = .value)) +
stat_lineribbon() +
geom_point(aes(x = water_c, y = N),
data = mite_community_abd,
inherit.aes = FALSE)
Another great way to model overdispersion is via the Negative Binomial distribution. Look at the Stan documentation for neg_binomial_2_log and adapt your model to use it (don’t forget to drop the random effect when you do!).