load("./calculate_fuel_load.rda")
source("./scripts/get_fuel_data.r")
transectid <- c("site", "treatment", "corner", "azi")
library(glmmTMB)
library(ggdist)
library(bayesplot)
library(brms)
library(tidybayes)
library(DHARMa)
library(dplyr)
library(ggplot2)
library(marginaleffects)
library(tidyr)
library(stringr)
library(purrr)Exploration of frequentist and Bayesian approaches can lead to more robust conclusions.
The often recommended Pillai’s Trace Test is robust to the normality assumption. Follow up with linear discriminant analysis, or multiple one-way anovas dependidng on research question. Using a Bonferroni correction for rejecting the null of alpha / m, for m hypothesis, we get an alpha of 0.008 for an alpha of 0.05 and 6 tests.
This suggests that it is unlikely that all treatemnts are equal.
load_vars <- c("onehr", "tenhr", "hundhr", "dufflitter", "thoushr", "veg")
tl <- load2("wide", "pre", treatment, all_of(load_vars))
myexpr <- expr(cbind(!!!syms(load_vars)) ~ treatment)
test1 <- do.call("manova", list(myexpr, data = quote(tl)))
summary(test1) Df Pillai approx F num Df den Df Pr(>F)
treatment 3 0.35947 2.7454 18 363 0.0001889 ***
Residuals 124
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1One way anova (using the welch test) can either assume constant variance or not. A levene test (using median) indicates onehr, tenhr, and veg may all have different variances between groups.
The one-way anova test results are the same though between equal and unequal variance assumptions. These tests support the notion that we can’t assume that the mean vegetatvie and onehr fuel loading are equal across all treatments, but there isn’t such evidence for the other fuel loading classes.
d <- load2("long", "pre", treatment, all_of(load_vars)) |> group_by(class)
d |>
nest() |>
rowwise() |>
transmute(
levene = car::leveneTest(load ~ factor(treatment), data)[[3]][1],
welch_uneq_var = oneway.test(load ~ treatment, data)$p.value,
welch_eq_var = oneway.test(
load ~ treatment,
var.equal = TRUE,
data = data
)$p.value,
) |>
knitr::kable(digits = 3)| class | levene | welch_uneq_var | welch_eq_var |
|---|---|---|---|
| onehr | 0.043 | 0.001 | 0.000 |
| tenhr | 0.020 | 0.491 | 0.596 |
| hundhr | 0.937 | 0.648 | 0.636 |
| dufflitter | 0.118 | 0.152 | 0.136 |
| thoushr | 0.955 | 0.919 | 0.955 |
| veg | 0.006 | 0.010 | 0.001 |
We can use the Games Howell test for pairwise comparisons to follow up on the welches test for differences between means when there is unequal variance among groups. These p-values provide evidence that for onehr fuels, the mean value of ha is greater than gs and ld, and the mean value for hd is also greater than gs and ld. Also, for vegetation, gs is greater than ha only. While this test is robust to the assumptions of normality, some of our data is highly skewed. Also, because of the nesting of our data, observations are not independent, so our effective sample size is not what is assumed by this test.
gh_test <- d |>
rstatix::games_howell_test(load ~ treatment) |>
filter(p.adj.signif != "ns") |>
rstatix::add_y_position(scales = "free", step.increase = 0.5)
ggpubr::ggboxplot(d, x = "treatment", y = "load", facet.by = "class") +
facet_wrap(~class, scales = "free") +
ggpubr::stat_pvalue_manual(gh_test, label = "p.adj") +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
We have transects nested within plot corners, corners nested within plots, and plots nested within sites. We would like to detect a treatment effect, while accounting for the non-independence of this nested data structure. The following model, I believe, captures this grouping structure.
form <- load ~ treatment + (1 | site / treatment / corner)This will estimate a group-wise intercept adjustments for each site, plot, and corner, based on modeled variances for each of these grouping levels.
d <- load2("long", "pre", site, treatment, corner, all_of(load_vars)) |>
group_by(class)
m1 <- d |>
nest() |>
rowwise() |>
transmute(
mod = list(lme4::lmer(form, data = data)),
emmeans = list(emmeans::emmeans(mod, "treatment")),
pairs = list(as_tibble(pairs(emmeans, infer = TRUE)))
)boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')Pairwise comparisons with Tukey adjustment for each of 6 multilevel models representing different fuel loading classes reveals that the only evidence for differences in means among treatments is with vegetation between the gs and ha treatments. Another sizeable difference in means is between gs and ha for the onehr fuels (Figure 6.2).
select(m1, pairs) |>
unnest(pairs) |>
filter(p.value <= 0.05) |>
knitr::kable(digits = 3)Adding missing grouping variables: `class`| class | contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|---|
| veg | gs - ha | 25.31 | 6.879 | 9 | 3.835 | 46.786 | 3.679 | 0.022 |
group_map(
m1,
~ plot(.x$emmeans[[1]], comparisons = TRUE) + ggtitle(.x$class)
) |>
patchwork::wrap_plots()
Hypothesis testing with multi-level models is not as straight forward with multi-level models. The problem, explained here is two fold. For GLMMs and unbalanced experimental designs, the null distribution for the F-statistic may not be F-distributed.
For us, we have a balanced design (I think) and so the F-statistic should be F distributed and degrees of freedom should be clear from the details of the design. Because of our balanced design, the Kenward-Rogers approach and “inner-outer” design approach (which is used by nlme::lme) give the same result of 9 DF for all of the pairwise tests.
Using the package pbkrtest we can get parametric bootstrap liklihood ratio statistics and test this statistic in a number of different ways. The PBtest should probably be the most reliable, but I’ve included descriptions of the others from the package documentation for reference. I’m also including an F-test in which degrees of freedom are estimated with Kenward-Rogers approach.
d <- load2("long", "pre", all_of(c(transectid, load_vars)))
dd <- ungroup(d) |> split(~class)
if (file.exists("pmod.rda")) {
load("pmod.rda")
} else {
cluster <- parallel::makeCluster(rep("localhost", parallel::detectCores()))
pmod <- imap(dd, function(d, i) {
form <- load ~ treatment + (1 | site / treatment / corner)
amod <- lme4::lmer(form, d, REML = FALSE)
nmod <- update(amod, . ~ . - treatment)
krtest <- pbkrtest::KRmodcomp(amod, nmod) |>
pluck("test", \(x) slice(x, 1)) |>
rename(df = ndf) |>
rownames_to_column("test")
pbkrtest::PBmodcomp(amod, nmod, cl = cluster) |>
pluck(summary, "test") |>
rownames_to_column("test") |>
bind_rows(krtest) |>
mutate(class = i, .before = 1) |>
relocate(c(df, ddf, F.scaling), .after = stat)
}) |>
list_rbind() |>
filter(test %in% c("LRT", "PBtest"))
parallel::stopCluster(cluster)
save(pmod, file = "pmod.rda")
}
pmod |> knitr::kable(digits = c(NA, NA, 2, 1, 2, 1, 4))| class | test | stat | df | ddf | F.scaling | p.value |
|---|---|---|---|---|---|---|
| onehr | LRT | 7.48 | 3 | NA | NA | 0.0580 |
| onehr | PBtest | 7.48 | NA | NA | NA | 0.1678 |
| tenhr | LRT | 0.77 | 3 | NA | NA | 0.8576 |
| tenhr | PBtest | 0.77 | NA | NA | NA | 0.9018 |
| hundhr | LRT | 0.70 | 3 | NA | NA | 0.8725 |
| hundhr | PBtest | 0.70 | NA | NA | NA | 0.9101 |
| dufflitter | LRT | 6.05 | 3 | NA | NA | 0.1093 |
| dufflitter | PBtest | 6.05 | NA | NA | NA | 0.1109 |
| thoushr | LRT | 0.22 | 3 | NA | NA | 0.9739 |
| thoushr | PBtest | 0.22 | NA | NA | NA | 0.9657 |
| veg | LRT | 11.70 | 3 | NA | NA | 0.0085 |
| veg | PBtest | 11.70 | NA | NA | NA | 0.0340 |
Taking a look at residual vs. fitted and qqplots of the model, it looks like our residuals are not normally distributed and there is not constant variance.
# These are functions to plot for each model, residuals vs fitted and normal
# quantiles. The third function is a wrapper to do both.
resid_plot <- function(data) {
data |>
ggplot(aes(fitted, resid)) +
geom_point() +
facet_wrap(~class, scales = "free") +
geom_hline(yintercept = 0)
}
qq_plot <- function(data) {
data |>
ggplot(aes(sample = resid)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~class, scales = "free")
}
resid_qq_plot <- function(data) {
data <- unnest(data, c(resid, fitted))
list(
a = resid_plot(data),
b = qq_plot(data)
)
}
d <- load2("long", "pre", all_of(c(transectid, load_vars))) |>
group_by(class) |>
nest() |>
rowwise()form <- load ~ treatment + (1 | site / treatment / corner)
mod1 <- d |>
mutate(
mod = list(lme4::lmer(form, data)),
fitted = list(fitted(mod)),
resid = list(resid(mod, type = "pearson", scaled = TRUE)),
.keep = "unused"
)boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')resid_qq_plot(mod1) |> patchwork::wrap_plots(ncol = 1)
lme4 The residuals are not homogenous.
I’ll try to control the variance by refitting the model with nlme::lme and using the weights argument. I’ll be using the pearson residuals which are corrected for heteroscedasticity.
I had to use the control argument sigma = 1 for the model to fit. I’m not sure why, I read it in the documentation for nlme::varConstProp. I’m modeling variance as a constant proportion of the fitted values of the model. This seems to have cleaned up the variance, but the the model for 1,000-hr fuels is causing a “singular matrix” error. Perhaps this is because of the excess of zeros?
Additionaly, residuals are still not normally distributed.
mod2 <- d |>
filter(class != "thoushr") |>
mutate(
mod = list(nlme::lme(
fixed = load ~ treatment,
random = ~ 1 | site / treatment / corner,
data = data,
weights = nlme::varConstProp(),
control = nlme::lmeControl(sigma = 1)
)),
fitted = list(fitted(mod)),
resid = list(resid(mod, type = "pearson")),
.keep = "unused"
)
patchwork::wrap_plots(resid_qq_plot(mod2), ncol = 1)
nlme. The (scaled) residuals are more homogenous now.
first, I want to see if the models produced by lme and lmer are equivalent they seem equivalent enough, although the random effects variances estimated by lmer are somewhat smaller.
d |>
mutate(
mod1 = list(lme4::lmer(form, data)),
mod2 = list(nlme::lme(
fixed = load ~ treatment,
random = ~ 1 | site / treatment / corner,
data = data
)),
.keep = "unused"
) |>
pivot_longer(-class, names_to = "model") |>
rowwise() |>
mutate(s = list(broom.mixed::tidy(value, effect = "fixed"))) |>
select(class, model, s) |>
unnest(everything()) |>
arrange(class, term) |>
print(n = 48)boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')# A tibble: 48 × 9
class model effect term estimate std.error statistic df p.value
<fct> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 onehr mod1 fixed (Interc… 0.598 0.191 3.13 NA NA
2 onehr mod2 fixed (Interc… 0.598 0.191 3.13 64 2.65e- 3
3 onehr mod1 fixed treatme… 0.603 0.256 2.35 NA NA
4 onehr mod2 fixed treatme… 0.603 0.256 2.35 9 4.33e- 2
5 onehr mod1 fixed treatme… 0.439 0.256 1.71 NA NA
6 onehr mod2 fixed treatme… 0.439 0.256 1.71 9 1.21e- 1
7 onehr mod1 fixed treatme… 0.0697 0.256 0.272 NA NA
8 onehr mod2 fixed treatme… 0.0697 0.256 0.272 9 7.92e- 1
9 tenhr mod1 fixed (Interc… 3.75 0.784 4.78 NA NA
10 tenhr mod2 fixed (Interc… 3.75 0.784 4.78 64 1.06e- 5
11 tenhr mod1 fixed treatme… -0.550 1.11 -0.496 NA NA
12 tenhr mod2 fixed treatme… -0.550 1.11 -0.496 9 6.32e- 1
13 tenhr mod1 fixed treatme… -0.828 1.11 -0.747 NA NA
14 tenhr mod2 fixed treatme… -0.828 1.11 -0.747 9 4.74e- 1
15 tenhr mod1 fixed treatme… -0.364 1.11 -0.329 NA NA
16 tenhr mod2 fixed treatme… -0.364 1.11 -0.329 9 7.50e- 1
17 hundhr mod1 fixed (Interc… 11.8 2.69 4.39 NA NA
18 hundhr mod2 fixed (Interc… 11.8 2.69 4.39 64 4.34e- 5
19 hundhr mod1 fixed treatme… -2.09 3.73 -0.562 NA NA
20 hundhr mod2 fixed treatme… -2.09 3.73 -0.562 9 5.88e- 1
21 hundhr mod1 fixed treatme… -2.34 3.73 -0.627 NA NA
22 hundhr mod2 fixed treatme… -2.34 3.73 -0.627 9 5.46e- 1
23 hundhr mod1 fixed treatme… -2.27 3.73 -0.608 NA NA
24 hundhr mod2 fixed treatme… -2.27 3.73 -0.608 9 5.58e- 1
25 dufflitter mod1 fixed (Interc… 48.4 5.94 8.14 NA NA
26 dufflitter mod2 fixed (Interc… 48.4 5.94 8.14 64 1.83e-11
27 dufflitter mod1 fixed treatme… -8.15 6.15 -1.33 NA NA
28 dufflitter mod2 fixed treatme… -8.15 6.15 -1.33 9 2.17e- 1
29 dufflitter mod1 fixed treatme… 6.58 6.15 1.07 NA NA
30 dufflitter mod2 fixed treatme… 6.58 6.15 1.07 9 3.12e- 1
31 dufflitter mod1 fixed treatme… -3.50 6.15 -0.570 NA NA
32 dufflitter mod2 fixed treatme… -3.50 6.15 -0.570 9 5.83e- 1
33 thoushr mod1 fixed (Interc… 44.9 17.1 2.63 NA NA
34 thoushr mod2 fixed (Interc… 44.9 17.1 2.63 64 1.08e- 2
35 thoushr mod1 fixed treatme… -3.27 24.2 -0.135 NA NA
36 thoushr mod2 fixed treatme… -3.27 24.2 -0.135 9 8.95e- 1
37 thoushr mod1 fixed treatme… -5.41 24.2 -0.224 NA NA
38 thoushr mod2 fixed treatme… -5.41 24.2 -0.224 9 8.28e- 1
39 thoushr mod1 fixed treatme… -9.65 24.2 -0.399 NA NA
40 thoushr mod2 fixed treatme… -9.65 24.2 -0.399 9 6.99e- 1
41 veg mod1 fixed (Interc… 38.2 6.12 6.24 NA NA
42 veg mod2 fixed (Interc… 38.2 6.12 6.24 64 3.98e- 8
43 veg mod1 fixed treatme… -25.3 6.88 -3.68 NA NA
44 veg mod2 fixed treatme… -25.3 6.88 -3.68 9 5.08e- 3
45 veg mod1 fixed treatme… -20.9 6.88 -3.04 NA NA
46 veg mod2 fixed treatme… -20.9 6.88 -3.04 9 1.40e- 2
47 veg mod1 fixed treatme… -17.0 6.88 -2.47 NA NA
48 veg mod2 fixed treatme… -17.0 6.88 -2.47 9 3.56e- 2Now, lets compare the two lme models using AIC. I’m fitting with REML because I’m not changing the fixed effects structure. I still have to remove 1,000-hr fuels from the mix.
This indicates that the model with modeled variance fits the data better than the model without.
mod_c <- d |>
filter(class != "thoushr") |>
mutate(
unweighted = list(nlme::lme(
fixed = load ~ treatment,
random = ~ 1 | site / treatment / corner,
data = data,
)),
weighted = list(nlme::lme(
fixed = load ~ treatment,
random = ~ 1 | site / treatment / corner,
data = data,
weights = nlme::varConstProp(),
control = nlme::lmeControl(sigma = 1),
)),
.keep = "unused"
)
mod_c |>
mutate(aic = list(across(matches("weight"), ~ AIC(.x)))) |>
select(aic) |>
unnest(aic) |>
knitr::kable()Adding missing grouping variables: `class`| class | unweighted | weighted |
|---|---|---|
| onehr | 251.8197 | 220.0447 |
| tenhr | 589.9042 | 486.1699 |
| hundhr | 892.3760 | 860.7936 |
| dufflitter | 1180.4072 | 1177.7730 |
| veg | 1171.0855 | 1122.1291 |
Does this change our conclusions about the effect of the treatment?
No, it doesn’t seem to make much of a difference.
d <- load2("long", "pre", site, treatment, corner, all_of(load_vars))
m2 <- d |>
group_by(class) |>
nest() |>
rowwise() |>
filter(class != "thoushr") |>
transmute(
mod = list(nlme::lme(
fixed = load ~ treatment,
random = ~ 1 | site / treatment / corner,
data = data,
weights = nlme::varConstProp(),
control = nlme::lmeControl(sigma = 1),
)),
emmeans = list(emmeans::emmeans(mod, "treatment")),
pairs = list(as_tibble(pairs(emmeans, infer = TRUE)))
)select(m2, pairs) |>
unnest(pairs) |>
filter(p.value <= 0.05) |>
knitr::kable(digits = 3)Adding missing grouping variables: `class`| class | contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|---|
| veg | gs - ha | 25.31 | 7.12 | 9 | 3.083 | 47.538 | 3.555 | 0.026 |
group_map(
m2,
~ plot(.x$emmeans[[1]], comparisons = TRUE) + ggtitle(.x$class)
) |>
patchwork::wrap_plots()
I’m not sure I’m using the correct random effects specification. The somewhat confusing thing is that I have a random effects nested above and below my fixed effect. This means that when I specify my random effect using the nesting notation: 1 | site/treatment/corner, I’m estimating a variance for corner:treatment:site, treatment:site, and site. The interaction of treatment and site here is analagous to a plot effect, of which there are 16.
It has been difficult fitting these models zeros, non-constant variance, and non-normal response. A more flexible approach may be to implement my models in a Bayesian context. This has the added benefits of ease of interpretability of the results and quantified estimates of uncertainty.
I’ll use brms with the same formula I used for the lmm above. First we’ll reload our data.
All the data is nested to facilitate modeling each fuel class separately. We’ll look at the average loading to get an idea of the data.
d <- load2("long", "pre", all_of(c(transectid, load_vars))) |>
group_by(class) |>
nest() |>
rowwise()I started out using a Gaussian model for load, with mostly default priors on the random effects. This was for convenience more than for anything. A model with positive support only, and more approriate priors would improve the model. This parameterization also results in the first treatment level (gs) being interpreted as the global intercept. This implies the unfortunate assumption that the other levels entail more variability as their priors result from the combination of the intercept prior as well as the treatment prior, wheres the prior for the first level only contains the variability in the intercept. Thus, parameterizing the model by removing the intercept for the main effect would allow equivalent interpretations of all four treatment levels. I use this approach below in the Gamma model.
We are using mostly uninformative priors for our un-pooled (fixed effect) estimates of treatment intercepts. They are all set as normal distributions, centered at the median of the fuel load with a sd of 2.5 times the sd of the data. While it is not possible to have negative values here, a current limitation of the brms package is that you can’t put bounds on individual coefficients, and we would not want to contstrain our treatment effects to be positive, as they are centered around the mean and will be both positive and negative.
# add calcualted priors to the data
d <- mutate(
d,
priors = list(
brms::set_prior(
str_glue(
"normal( {round(median(data$load))}, {round(2.5 * sd(data$load))} )"
),
class = "b",
coef = "Intercept"
) +
brms::set_prior(
str_glue(
"normal( 0, {round(1.5 * sd(data$load))} )"
),
class = "b"
)
)
)
# Plot the priors
d |>
mutate(
prior_dist = list(
tidybayes::parse_dist(priors) |>
rename(dist_class = class)
)
) |>
unnest(prior_dist) |>
ggplot(aes(ydist = .dist_obj, color = interaction(coef, dist_class))) +
stat_slab(normalize = "panels", fill = NA) +
stat_pointinterval(
position = position_dodge(width = 0.3, preserve = "single")
) +
geom_text(
aes(label = format(.dist_obj), y = mean(.dist_obj), x = 0.97, vjust = 0),
position = position_dodge(width = 0.3),
show.legend = FALSE
) +
coord_flip() +
facet_wrap(~class, scales = "free_x") +
scale_color_hue(name = "prior", labels = c("treatment", "intercept", "sd"))
Now we’ll fit the model.
form <- load ~ 0 + Intercept + treatment + (1 | site / treatment / corner)
bf2 <- mutate(
d,
mod = list(brms::brm(
form,
data,
warmup = 3000,
iter = 4000,
cores = 4,
control = list(adapt_delta = 0.99),
prior = priors,
sample_prior = TRUE,
file = paste0("fits/bf2_", class)
))
)# Plod density of model priors against posterior to get a sense of the
# informativeness of the prior. Average over all fixed effects.
plot_pri_post <- function(mod) {
tidy_draws(mod) |>
select(c(matches("b_"), matches("sd"), matches("sigma"), matches("hu"))) |>
pivot_longer(everything()) |>
mutate(
name = if_else(
str_starts(name, "prior"),
name,
paste0("posterior_", name)
)
) |>
separate_wider_regex(
name,
c(phase = "prior|posterior", "_", name = ".*")
) |>
mutate(name = str_remove(name, "__.*")) |>
ggplot(aes(x = value, color = phase, fill = phase)) +
stat_slab(normalize = "panels", alpha = 3 / 4) +
stat_pointinterval(
position = position_dodge(width = 0.3, preserve = "single")
) +
facet_wrap(~name, scales = "free", labeller = fuel_class_labeller) +
labs(x = "Parameter value", y = "Density") +
scale_fill_manual(values = c("gray60", "gray30")) +
scale_color_manual(values = c("gray60", "gray30"))
}
posterior_predictive_check <- function(models) {
models <- mutate(
models,
pred = list(tidybayes::add_predicted_draws(data, mod, ndraws = 15)),
lims = list(tibble(
xmin = min(
quantile(data$load, .001),
quantile(pred$.prediction, .001)
),
xmax = max(
quantile(data$load, .999),
quantile(pred$.prediction, .999)
)
))
)
ggplot() +
geom_line(
data = unnest(models, data),
aes(x = load),
stat = "density",
size = 1
) +
geom_line(
data = unnest(models, pred),
aes(x = .prediction, group = .draw),
stat = "density",
alpha = 0.15,
size = 1
) +
facet_wrap(~class, scales = "free", labeller = fuel_class_labeller) +
coord_cartesian_panels(
panel_limits = unnest(select(models, lims), lims)
)
}Here is a posterior predictive check.
posterior_predictive_check(bf2)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.Adding missing grouping variables: `class`
The Gaussian distribution is symmetric and doesn’t capture well the peak near zero and the long right tail of our observed values for most of the fuel classes. It also dramatically overpredicts negative values (which are absent from our data, despite the fact that the density smoothing of the observed values seems to suggest there are some.)
The fact that the model predicts negative values suggests that it is not right for our data, and could potentially be biasing comparissons between treatments. A Gamma distribution for the response makes more sense.
The gamma model has support for only postive values, which makes sense for our weight per area data. Initial data exploration also revealed that the gamma appears to be a good fit for our data. This makes sense, as our data is fundamentally transformed count data, and the gamma distribution is the continuous generalization of the negative binomial, which is used for modeling count data.
I’m assuming the outcome is hurdle gamma distributed.
\[ \operatorname{HurdleGamma}(y \mid \mu, \text{shape}, \text{hu}) = \begin{cases} \text{hu} &\quad\text{if } y = 0, \text{ and}\\ (1 - \text{hu}) \operatorname{Gamma}(y \mid \alpha, \beta) &\quad\text{if } y > 0, \end{cases} \]
Where \(\mu = \frac{\alpha}{\beta}\), and \(\text{shape} = \beta\). This is a mixture model where the proportion of zeros are estimated and the rest of the observations are assumed to be gamma distributed. Because the gamma distribution doesn’t incldue support for zero, we don’t have to worry about any “overlapping” zero predictions from the gamma portion of the model.
Our model form is as follows:
\[ \begin{align*} y_i &\sim \operatorname{HurdleGamma}(\mu, \text{shape}, \text{hu})\\ \log(\mu_i) &= \bar{\alpha} + U_{\text{site}[i]}\sigma_{u} + V_{\text{plot}[i]}\sigma_{v} + W_{\text{corner}[i]}\sigma_{w} + \beta_{treatment[i]}\\ \bar{\alpha} &\sim \operatorname{normal}(M, S)\\ U_j &\sim \operatorname{normal}(0, 1) \qquad \text{for } j=1 \dots 4\\ V_j &\sim \operatorname{normal}(0, 1) \qquad \text{for } j=1 \dots 16\\ W_j &\sim \operatorname{normal}(0, 1) \qquad \text{for } j=1 \dots 64\\ \sigma_u &\sim \operatorname{normal^+}(0, 1)\\ \sigma_v &\sim \operatorname{normal^+}(0, 1)\\ \sigma_w &\sim \operatorname{normal^+}(0, 1)\\ \beta_j &\sim \operatorname{normal}(0, 1) \qquad \text{for } j=1 \dots 4\\ \text{shape} &\sim \operatorname{Gamma}(0.01, 0.01)\\ \text{hu} &\sim \operatorname{beta}(1, 1) \end{align*} \tag{6.1}\]
Here is an interpretation for each line of the above model.
brms. It is chosen to be minimally informative.One problem with these independent prior assumptions is that the multiplicative effect of the terms on the response scale (due to exponentiation) can lead to impossibly large predictiontions, particularly when multiple random or fixed effect parameters are sampled simultaneosly in their upper tails. For instance, if the effects of site, plot, and corner where sampled at 2 (on the log scale) and the effect of treatment, at 1.8, this could lead to a predicted \(\mu\) more than 2,400 times greater than the grand mean.
A better choice would be to use a joint prior, like the Dirichlet, so a prior could be set on the total variance and the individual component variances would vary (in either a eqaul, or unequally weighted fashioin) in their repsective proportions of that total.
In this regard, with either independent, or independet priors, further work could be done on establshing the relative importance of each variance component and setting priors that reflect these differences. For our nesting structure, this means defining the length scale over which fuels are expected to vary most/least, from 10 m to 50 m to > 100 m, to thousands of meters.
First I’ll set my prior. This is complicated because we are now working with a log link. According to Solomon Kurz, this is how to transform your mean and sd for a normal distribution on the identity scale, to the equivalent normal distribution on the log scale (a lognormal distribution). This is what we’ll use for our prior on the grand mean (\(\bar{\alpha}\)).
d <- load2("long", "pre", all_of(c(transectid, load_vars))) |>
group_by(class) |>
nest() |>
rowwise()
lnp <- function(data) {
# Desired values
m <- mean(data)
s <- 2.5 * sd(data)
# use the equations
mu <- log(m / sqrt(s^2 / m^2 + 1))
sigma <- sqrt(log(s^2 / m^2 + 1))
# output mu and sigma on lognormals own scale
list2env(list(mu = round(mu, 2), sigma = round(sigma, 2)))
}# data frame to plot the priors
pd <- d |>
mutate(
priors = list(brms::set_prior(
str_glue("normal({mu}, {sigma})", .envir = lnp(data$load))
)),
prior_dist = list(tidybayes::parse_dist(priors)),
lims = list(
tibble(
xmin = quantile(exp(prior_dist$.dist_obj), .01),
xmax = quantile(exp(prior_dist$.dist_obj), .99)
)
)
) |>
rename(fuel_class = class)
ggplot(data = unnest(pd, prior_dist)) +
tidybayes::stat_halfeye(
aes(xdist = exp(.dist_obj)),
normalize = "panels"
) +
geom_text(
aes(
label = format(.dist_obj),
x = mean(exp(.dist_obj)),
y = 0.97,
hjust = 0
)
) +
facet_wrap(~fuel_class, scales = "free_x", labeller = fuel_class_labeller) +
coord_cartesian_panels(panel_limits = unnest(select(pd, lims), lims)) +
labs(x = expression(exp(bar(alpha))))
The below code ran with very many divergent transitions and the results did not really make sense. This may have to do with a comment I found here: independent priors specified on the log link usually don’t sample.
An alternative is refitting the whole model with different priors and comparing the output to infer the prior influence. I’m not going to do that right now.
I will though, sample from the priors while fitting the model in order to plot those along with the posterior for our variables of interest.
Here I actually fit several models. They were each explored to some degree in an iterative process of model fitting and checking. For the sake of brevity, I I have focussed on just one of them (bf4a) for displaying final results. This model is the one described mathematically above.
bf4a <- mutate(
d,
priors = list(
set_prior(
str_glue("normal({mu}, {sigma})", .envir = lnp(data$load)),
nlpar = "a",
coef = "Intercept"
) +
set_prior("normal(0, 1)", nlpar = "a", class = "sd") +
set_prior("normal(0, 1)", nlpar = "b", coef = "treatmentgs") +
set_prior("normal(0, 1)", nlpar = "b", coef = "treatmentha") +
set_prior("normal(0, 1)", nlpar = "b", coef = "treatmenthd") +
set_prior("normal(0, 1)", nlpar = "b", coef = "treatmentld")
),
mod = list(brms::brm(
brms::bf(
load ~ a + b,
a ~ 1 + (1 | site / treatment / corner),
b ~ 0 + treatment,
nl = TRUE
),
data,
warmup = 4000,
iter = 5000,
cores = 4,
control = list(adapt_delta = 0.99),
family = brms::hurdle_gamma(),
prior = priors,
sample_prior = TRUE,
file = paste0("fits/bf4a_", class)
))
)bf3 <- mutate(
d,
priors = list(brms::set_prior(
str_glue("normal({mu}, {sigma})", .envir = lnp(data$load))
)),
mod = list(brms::brm(
form,
data,
warmup = 5000,
iter = 6000,
cores = 4,
control = list(adapt_delta = 0.99),
family = brms::hurdle_gamma(),
prior = priors,
file = paste0("fits/bf3_", class)
))
)
bf4 <- mutate(
d,
priors = list(
set_prior(
str_glue("normal({mu}, {sigma})", .envir = lnp(data$load)),
nlpar = "a",
coef = "Intercept"
) +
set_prior("student_t(3, 0, 0.35)", nlpar = "a", class = "sd") +
set_prior("student_t(3, 0, 0.35)", nlpar = "b", coef = "treatmentgs") +
set_prior("student_t(3, 0, 0.35)", nlpar = "b", coef = "treatmentha") +
set_prior("student_t(3, 0, 0.35)", nlpar = "b", coef = "treatmenthd") +
set_prior("student_t(3, 0, 0.35)", nlpar = "b", coef = "treatmentld")
),
mod = list(brms::brm(
brms::bf(
load ~ a + b,
a ~ 1 + (1 | site / treatment / corner),
b ~ 0 + treatment,
nl = TRUE
),
data,
warmup = 4000,
iter = 5000,
cores = 4,
control = list(adapt_delta = 0.99),
family = brms::hurdle_gamma(),
prior = priors,
sample_prior = TRUE,
file = paste0("fits/bf4_", class)
))
)
bf5 <- mutate(
d,
priors = list(brms::set_prior(
str_glue("normal({mu}, {sigma})", .envir = lnp(data$load))
)),
mod = list(brms::brm(
brms::bf(
load ~ 0 + Intercept + (1 | site / treatment / corner) + (1 | treatment)
),
data,
warmup = 4000,
iter = 5000,
cores = 4,
control = list(adapt_delta = 0.99),
family = brms::hurdle_gamma(),
prior = priors,
file = paste0("fits/bf5_", class)
))
)my_summary <- function(mod) {
as_draws_df(mod) |>
select(starts_with(c("b_", "sd_")), shape, hu) |>
posterior::summarize_draws(
"mean",
~ posterior::quantile2(.x, c(0.05, .5, .95)),
"rhat",
"ess_bulk",
"ess_tail"
)
}
bf4a |>
mutate(
summary = list(my_summary(mod))
) |>
select(summary) |>
unnest(summary) |>
mutate(
variable = str_remove_all(variable, "b_._|_a_")
) |>
knitr::kable(digits = c(NA, NA, 2, 2, 2, 2, 2, 0, 0))| class | variable | mean | q5 | q50 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|
| onehr | Intercept | -0.32 | -1.17 | -0.32 | 0.53 | 1.00 | 1261 | 1666 |
| onehr | treatmentgs | -0.31 | -1.15 | -0.30 | 0.52 | 1.00 | 1278 | 1690 |
| onehr | treatmentha | 0.44 | -0.39 | 0.44 | 1.24 | 1.00 | 1311 | 1640 |
| onehr | treatmenthd | 0.27 | -0.56 | 0.27 | 1.08 | 1.00 | 1363 | 1569 |
| onehr | treatmentld | -0.13 | -0.99 | -0.13 | 0.72 | 1.00 | 1256 | 1500 |
| onehr | sd_site_Intercept | 0.38 | 0.05 | 0.32 | 0.92 | 1.00 | 1131 | 1232 |
| onehr | sd_site:treatment_Intercept | 0.35 | 0.14 | 0.34 | 0.60 | 1.00 | 1221 | 1245 |
| onehr | sd_site:treatment:corner_Intercept | 0.12 | 0.01 | 0.11 | 0.28 | 1.00 | 1200 | 1374 |
| onehr | shape | 2.35 | 1.90 | 2.34 | 2.86 | 1.00 | 3019 | 2789 |
| onehr | hu | 0.01 | 0.00 | 0.01 | 0.02 | 1.00 | 3317 | 1647 |
| tenhr | Intercept | 1.06 | 0.26 | 1.07 | 1.86 | 1.00 | 1384 | 1902 |
| tenhr | treatmentgs | 0.28 | -0.53 | 0.28 | 1.07 | 1.00 | 1363 | 1914 |
| tenhr | treatmentha | 0.07 | -0.78 | 0.07 | 0.88 | 1.00 | 1477 | 2044 |
| tenhr | treatmenthd | 0.00 | -0.82 | -0.01 | 0.81 | 1.00 | 1418 | 2169 |
| tenhr | treatmentld | 0.10 | -0.72 | 0.11 | 0.88 | 1.00 | 1425 | 1911 |
| tenhr | sd_site_Intercept | 0.25 | 0.02 | 0.18 | 0.71 | 1.00 | 1396 | 1959 |
| tenhr | sd_site:treatment_Intercept | 0.38 | 0.17 | 0.37 | 0.63 | 1.00 | 912 | 669 |
| tenhr | sd_site:treatment:corner_Intercept | 0.11 | 0.01 | 0.10 | 0.27 | 1.00 | 1302 | 1596 |
| tenhr | shape | 2.56 | 2.01 | 2.54 | 3.18 | 1.00 | 3337 | 3046 |
| tenhr | hu | 0.06 | 0.03 | 0.06 | 0.10 | 1.00 | 5628 | 2452 |
| hundhr | Intercept | 2.25 | 1.40 | 2.26 | 3.06 | 1.00 | 2087 | 2520 |
| hundhr | treatmentgs | 0.26 | -0.56 | 0.26 | 1.10 | 1.00 | 2094 | 2651 |
| hundhr | treatmentha | 0.03 | -0.78 | 0.03 | 0.88 | 1.00 | 2195 | 2537 |
| hundhr | treatmenthd | 0.02 | -0.79 | 0.02 | 0.84 | 1.00 | 2116 | 2424 |
| hundhr | treatmentld | 0.15 | -0.66 | 0.15 | 1.01 | 1.00 | 2129 | 2714 |
| hundhr | sd_site_Intercept | 0.30 | 0.03 | 0.24 | 0.80 | 1.00 | 1772 | 2044 |
| hundhr | sd_site:treatment_Intercept | 0.25 | 0.05 | 0.24 | 0.47 | 1.00 | 1015 | 1407 |
| hundhr | sd_site:treatment:corner_Intercept | 0.15 | 0.01 | 0.13 | 0.32 | 1.00 | 1408 | 2164 |
| hundhr | shape | 2.80 | 2.18 | 2.77 | 3.52 | 1.00 | 3476 | 2851 |
| hundhr | hu | 0.12 | 0.08 | 0.12 | 0.17 | 1.00 | 9182 | 2704 |
| dufflitter | Intercept | 3.74 | 2.94 | 3.74 | 4.53 | 1.00 | 1027 | 1411 |
| dufflitter | treatmentgs | 0.10 | -0.66 | 0.11 | 0.87 | 1.00 | 1025 | 1469 |
| dufflitter | treatmentha | -0.07 | -0.84 | -0.06 | 0.68 | 1.00 | 1032 | 1524 |
| dufflitter | treatmenthd | 0.23 | -0.54 | 0.24 | 0.99 | 1.00 | 1058 | 1572 |
| dufflitter | treatmentld | 0.05 | -0.73 | 0.05 | 0.82 | 1.00 | 1031 | 1734 |
| dufflitter | sd_site_Intercept | 0.29 | 0.06 | 0.23 | 0.74 | 1.00 | 1237 | 1178 |
| dufflitter | sd_site:treatment_Intercept | 0.10 | 0.01 | 0.08 | 0.24 | 1.00 | 1539 | 1640 |
| dufflitter | sd_site:treatment:corner_Intercept | 0.10 | 0.01 | 0.09 | 0.23 | 1.00 | 993 | 1652 |
| dufflitter | shape | 3.93 | 3.15 | 3.90 | 4.85 | 1.00 | 2983 | 2787 |
| dufflitter | hu | 0.01 | 0.00 | 0.01 | 0.02 | 1.00 | 3668 | 1545 |
| thoushr | Intercept | 3.58 | 2.65 | 3.59 | 4.49 | 1.00 | 1268 | 1860 |
| thoushr | treatmentgs | 0.35 | -0.57 | 0.35 | 1.28 | 1.00 | 1340 | 2428 |
| thoushr | treatmentld | -0.09 | -1.02 | -0.09 | 0.85 | 1.00 | 1468 | 2023 |
| thoushr | treatmentha | -0.05 | -0.99 | -0.04 | 0.87 | 1.00 | 1668 | 2247 |
| thoushr | treatmenthd | 0.29 | -0.63 | 0.29 | 1.24 | 1.00 | 1534 | 2291 |
| thoushr | sd_site_Intercept | 0.40 | 0.03 | 0.33 | 1.03 | 1.00 | 1370 | 2076 |
| thoushr | sd_site:treatment_Intercept | 0.58 | 0.23 | 0.56 | 0.96 | 1.00 | 913 | 721 |
| thoushr | sd_site:treatment:corner_Intercept | 0.38 | 0.07 | 0.38 | 0.67 | 1.00 | 794 | 1242 |
| thoushr | shape | 1.41 | 1.06 | 1.39 | 1.81 | 1.00 | 2063 | 2865 |
| thoushr | hu | 0.25 | 0.19 | 0.25 | 0.32 | 1.00 | 5170 | 2843 |
| veg | Intercept | 2.75 | 1.84 | 2.75 | 3.66 | 1.00 | 1119 | 2021 |
| veg | treatmentgs | 0.53 | -0.32 | 0.52 | 1.39 | 1.01 | 1293 | 2170 |
| veg | treatmentha | -0.39 | -1.28 | -0.37 | 0.48 | 1.00 | 1274 | 1827 |
| veg | treatmenthd | 0.00 | -0.90 | 0.01 | 0.88 | 1.00 | 1071 | 2035 |
| veg | treatmentld | 0.21 | -0.70 | 0.21 | 1.09 | 1.00 | 1147 | 1846 |
| veg | sd_site_Intercept | 0.46 | 0.07 | 0.39 | 1.07 | 1.00 | 1377 | 1222 |
| veg | sd_site:treatment_Intercept | 0.38 | 0.06 | 0.37 | 0.73 | 1.00 | 716 | 801 |
| veg | sd_site:treatment:corner_Intercept | 0.51 | 0.23 | 0.52 | 0.74 | 1.00 | 598 | 552 |
| veg | shape | 1.97 | 1.47 | 1.95 | 2.56 | 1.00 | 1094 | 1691 |
| veg | hu | 0.02 | 0.00 | 0.01 | 0.04 | 1.00 | 4763 | 2070 |
While it was not possible to sample from the prior using the MCMC sampler which was likely due to the extremely long tails implied by our exponentiated, independent priors, another way of describing priors is simply sampling random from random number generators. The results of these random samples are plotted against the posterior distributions for the variables that we set priors for.
#| label: fig-gamma-pri-post
#| fig-cap: >
#| Posterior and prior distributions for variables we set priors for, for one
#| fuel class model.
#| fig-subcap:
#| - ""
#| - ""
#| - ""
#| - ""
#| - ""
#| - ""
plot_pri_post(bf4a$mod[[1]]) + ggtitle("1-hr")
plot_pri_post(bf4a$mod[[2]]) + ggtitle("10-hr")
plot_pri_post(bf4a$mod[[3]]) + ggtitle("100-hr")
plot_pri_post(bf4a$mod[[4]]) + ggtitle("duff & litter")
plot_pri_post(bf4a$mod[[5]]) + ggtitle("1000-hr")
plot_pri_post(bf4a$mod[[6]]) + ggtitle("Vegetation")
Adding missing grouping variables: `class`
The gamma model fits the data better. There are no predictions below zero anymore. It does seem like the gamma distribution tends to predict higher densities of lower values than we observed, as seen in the plots for the tenhr, thoushr, and veg fuel classes. But generally, the predictions appear to agree with the observed data pretty well.
# get posterior predictions for treatments for models for all fuel classes,
# ignoring random effects. These are predictions for the expected value across
# all sites.
predict_posterior_expected <- function(data, plot = TRUE, re_formula = NA) {
newdata <- tidyr::expand(data$data[[1]], nesting(treatment))
data <- mutate(
data,
.keep = "none",
pred = list(
tidybayes::epred_draws(
mod,
newdata,
re_formula = re_formula,
value = "pred"
)
),
lims = list(
tibble(xmin = 0, xmax = quantile(pred$pred, .995))
)
)
data
}
plot_posterior_predicted <- function(data) {
p <- data |>
unnest(pred) |>
ggplot(aes(pred, treatment)) +
tidybayes::stat_halfeye(normalize = "panels") +
facet_wrap(~class, scales = "free_x", labeller = fuel_class_labeller) +
coord_cartesian_panels(panel_limits = unnest(select(data, lims), lims)) +
scale_y_discrete(labels = toupper) +
labs(x = expression(Load ~ (Mg %.% ha^-1)), y = "Treatment")
p
}Adding missing grouping variables: `class`
These are the expected predictions, or predictions for the mean. It only includes the uncertainty in the mean and not the variance in predictions estimated by the model.
There is quite a bit of uncertainty about the mean all around, but there is a notable difference in that uncertainty among treatments for the onehr and veg fuel classes.
Table Table 6.7 shows these data in a tabular format. I’m using the highest density continuous interval because, while its hard to see in Figure 6.10, the highest desity region is actually slightly discontinuous.
expected_predictions |>
mutate(
summary = list(tidybayes::mode_hdci(pred))
) |>
select(summary) |>
unnest(summary) |>
select(treatment, prediction = pred, lower = .lower, upper = .upper) |>
knitr::kable(digits = 1)Adding missing grouping variables: `class`
Adding missing grouping variables: `class`| class | treatment | prediction | lower | upper |
|---|---|---|---|---|
| onehr | gs | 0.5 | 0.2 | 0.9 |
| onehr | ld | 0.6 | 0.3 | 1.1 |
| onehr | ha | 1.1 | 0.5 | 1.9 |
| onehr | hd | 0.9 | 0.4 | 1.6 |
| tenhr | gs | 3.5 | 1.8 | 5.8 |
| tenhr | ld | 2.8 | 1.5 | 4.9 |
| tenhr | ha | 2.7 | 1.4 | 4.7 |
| tenhr | hd | 2.6 | 1.3 | 4.4 |
| hundhr | gs | 10.5 | 5.5 | 16.7 |
| hundhr | ld | 9.4 | 5.2 | 15.6 |
| hundhr | ha | 8.6 | 4.7 | 13.9 |
| hundhr | hd | 8.3 | 4.2 | 13.2 |
| dufflitter | gs | 45.4 | 26.8 | 65.9 |
| dufflitter | ld | 43.3 | 26.0 | 64.4 |
| dufflitter | ha | 38.1 | 23.9 | 57.8 |
| dufflitter | hd | 51.4 | 30.4 | 77.2 |
| thoushr | gs | 33.1 | 11.7 | 78.0 |
| thoushr | ld | 21.2 | 7.2 | 53.3 |
| thoushr | ha | 21.4 | 6.8 | 52.9 |
| thoushr | hd | 31.4 | 13.0 | 78.2 |
| veg | gs | 25.1 | 9.4 | 49.9 |
| veg | ld | 17.9 | 7.1 | 35.8 |
| veg | ha | 9.3 | 4.0 | 20.1 |
| veg | hd | 14.7 | 5.4 | 29.4 |
predict_expected_contrasts <- function(
data,
rope_size,
plot = TRUE,
re_formula = NA
) {
# Assume treatment levels are the same for all models: they are.
newdata <- tidyr::expand(data$data[[1]], tidyr::nesting(treatment))
d <- data |>
mutate(
pred = list(
tidybayes::epred_draws(mod, newdata, re_formula = NA, value = "pred") |>
tidybayes::compare_levels(pred, by = treatment) |>
select(contrast = treatment, pred)
),
rope = rope_size * sd(data$load),
lims = list(
tibble(
xmin = quantile(pred$pred, .001),
xmax = quantile(pred$pred, .999)
)
),
.keep = "none"
)
}
plot_expected_contrasts <- function(data) {
p <- data |>
unnest(c(pred)) |>
ggplot(aes(x = pred, y = contrast)) +
tidybayes::stat_halfeye(normalize = "panels") +
geom_vline(aes(xintercept = rope)) +
geom_vline(aes(xintercept = -rope)) +
facet_wrap(~class, scales = "free_x", labeller = fuel_class_labeller) +
coord_cartesian_panels(
panel_limits = unnest(dplyr::select(data, lims), lims)
) +
scale_y_discrete(labels = toupper) +
labs(x = expression(Load ~ (Mg %.% ha^-1)), y = "Treatment")
p
}Adding missing grouping variables: `class`
Adding missing grouping variables: `class`
Adding missing grouping variables: `class`| class | contrast | prediction | lower | upper | prob |
|---|---|---|---|---|---|
| onehr | ha - gs | 0.53 | 0.05 | 1.27 | 0.99 |
| onehr | ha - ld | 0.47 | -0.07 | 1.22 | 0.96 |
| onehr | hd - gs | 0.40 | -0.06 | 1.02 | 0.97 |
| onehr | hd - ha | -0.18 | -0.91 | 0.53 | 0.72 |
| onehr | hd - ld | 0.24 | -0.20 | 0.90 | 0.91 |
| onehr | ld - gs | 0.07 | -0.27 | 0.55 | 0.74 |
| tenhr | ha - gs | -0.58 | -2.90 | 1.75 | 0.75 |
| tenhr | ha - ld | -0.01 | -2.15 | 1.97 | 0.54 |
| tenhr | hd - gs | -0.81 | -2.97 | 1.32 | 0.82 |
| tenhr | hd - ha | -0.28 | -2.25 | 1.73 | 0.59 |
| tenhr | hd - ld | -0.25 | -2.31 | 1.69 | 0.63 |
| tenhr | ld - gs | -0.42 | -2.94 | 1.68 | 0.73 |
| hundhr | ha - gs | -2.01 | -8.03 | 2.94 | 0.81 |
| hundhr | ha - ld | -1.13 | -6.18 | 4.34 | 0.69 |
| hundhr | hd - gs | -1.81 | -7.59 | 2.95 | 0.83 |
| hundhr | hd - ha | -0.03 | -4.79 | 4.68 | 0.52 |
| hundhr | hd - ld | -1.15 | -6.14 | 3.86 | 0.71 |
| hundhr | ld - gs | -0.92 | -6.91 | 4.51 | 0.66 |
| dufflitter | ha - gs | -5.61 | -23.14 | 6.25 | 0.86 |
| dufflitter | ha - ld | -4.15 | -18.16 | 9.96 | 0.78 |
| dufflitter | hd - gs | 6.63 | -10.40 | 23.21 | 0.80 |
| dufflitter | hd - ha | 13.96 | -2.39 | 30.07 | 0.97 |
| dufflitter | hd - ld | 7.78 | -6.65 | 26.21 | 0.88 |
| dufflitter | ld - gs | -2.38 | -17.13 | 13.10 | 0.63 |
| thoushr | ha - gs | -12.04 | -53.02 | 24.63 | 0.79 |
| thoushr | ha - ld | 1.79 | -29.22 | 32.75 | 0.54 |
| thoushr | hd - gs | -3.42 | -45.30 | 41.34 | 0.55 |
| thoushr | hd - ha | 6.08 | -24.61 | 48.67 | 0.76 |
| thoushr | hd - ld | 6.82 | -24.83 | 50.03 | 0.78 |
| thoushr | ld - gs | -9.37 | -55.97 | 19.62 | 0.82 |
| veg | ha - gs | -14.06 | -36.94 | -0.29 | 0.99 |
| veg | ha - ld | -6.80 | -22.39 | 4.75 | 0.95 |
| veg | hd - gs | -8.37 | -32.01 | 6.94 | 0.92 |
| veg | hd - ha | 4.25 | -5.81 | 17.37 | 0.86 |
| veg | hd - ld | -3.24 | -18.57 | 11.20 | 0.70 |
| veg | ld - gs | -5.12 | -28.88 | 11.62 | 0.81 |
bf4a |>
mutate(
R_squared = list(as_tibble(bayes_R2(mod)))
) |>
select(R_squared) |>
unnest(R_squared) |>
knitr::kable(digits = 2)Adding missing grouping variables: `class`| class | Estimate | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|
| onehr | 0.36 | 0.08 | 0.21 | 0.51 |
| tenhr | 0.23 | 0.08 | 0.09 | 0.39 |
| hundhr | 0.18 | 0.07 | 0.07 | 0.33 |
| dufflitter | 0.19 | 0.07 | 0.08 | 0.34 |
| thoushr | 0.17 | 0.08 | 0.06 | 0.35 |
| veg | 0.44 | 0.10 | 0.23 | 0.63 |
The final model is described in Equation 6.1. The description of the meaning of the priors in the following list should probably be included for explanation.
brms should be included in an appendix.save(
expected_predictions,
expected_contrasts,
fuel_class_labeller,
fuel_class_labels,
fuel_class_labels2,
file = "fuel_data_modeling.rda"
)