Willow Macroplot Statistical Analyses

Repository:
https://github.com/egage/EVMP

## custom ggplot emm

emm.func.tc.x.fenced <- function(model){
  emmeans(model, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  ggplot() +
  geom_linerange(aes(y=contrast, xmin=lower.HPD, xmax=upper.HPD), color="#babefa", size=2.5) +
  geom_point(aes(x=estimate, y=contrast)) +
  theme_minimal() +
  geom_vline(xintercept = 0, lty="dashed", size=1) +
  facet_wrap(~fenced, ncol=1)
}

# custom gradients and palettes
dbhclass.gradient5 <- c("#26185F", "#0095AF", "#9ADCBB", "#FCFFDD", "grey50")

colfunc1<-colorRampPalette(c("#26185F","#0095AF","#9ADCBB", "#FCFFDD"))
colfunc2<-colorRampPalette(c("#0095AF","#9ADCBB", "#FCFFDD"))
colfunc3<-colorRampPalette(c("#0095AF","#9ADCBB"))
colfunc3_RdBu <- colorRampPalette(c("#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF"))

# plot(rep(1,50),col=(colfunc1(50)), pch=19,cex=2)
# plot(rep(1,5),col=(colfunc1(5)), pch=19,cex=2)
pal.5a <- colfunc1(5)
pal.7a <- colfunc1(7)

# palette3 <- RColorBrewer::brewer.pal(6,"RdBu")
# RColorBrewer::display.brewer.pal(6,"RdBu")

Introduction

This document contains code for running statistical analyses of of Elk Vegetation Management Plan (EVMP) data collected through the 2018 sampling season. Analyses are limited to macroplot data from plots established in riparian areas in areas of the elk winter range and Kawuneeche Valley (i.e., “willow” plots, see Zeigenfuss et al. 2011). Analyses are provided as a supplement to the NPS NRR Report “Monitoring of Vegetation Response to Elk Population and Habitat Management in Rocky Mountain National Park - Analysis of Elk Vegetation Management Plan monitoring Data: 2008–2018.” For information on sampling and the broader analysis and interpretation of the results, refer to the project report, past analyses (Linda C. Zeigenfuss and Johnson 2015), and the original EVMP monitoring plan (Linda C. Zeigenfuss, Johnson, and Wiebe 2011).

Methods

The code and results presented here build on derived data produced in a separate code document focused on ingesting, compiling, and cleaning raw data provided by RMNP staff (Appendix B). Extensive “munging” of raw files was required, so the decision was made to separate code used to pre-process raw data from code used for the analyses presented here. Likewise, separate code documents were developed for other distinct elements of the larger EVMP (e.g., aspen and upland plots). In general, packages included in the “Tidyverse” ecosystem of R packages were used for data transformation and visualization (Wickham et al. 2019), although specialized packages particular to specific tasks (e.g., the “bayesplot” plot for visualization of Bayesian posterior distributions) were also used (Gabry et al. 2019).

Bayesian models view a model parameter \(\theta\) as a [random variable]. In constrast, frequentist models treat model parameters as unknown constant. Rather than estimating an unknown constant, Bayesian modeling focuses on an unknown distribution of parameter values.

Using Bayes’ Law, model parameters are estimated:

\[ \begin{equation} \label{eq:bayeslaw} P(\theta | y, X) = \frac{P(y|\theta,X) \cdot P(\theta)}{P(y)} \propto P(y|\theta,X) \cdot P(\theta). \end{equation} \]

Where:

\(X\): The observed data
\(y\): The outcomes
\(P(\theta|y,X)\): Our view of the possible value of our model parameters after seeing the data
\(P(y|\theta,X)\): The likelihood of seeing the current outcome conditioned on the data we have and a specific parameter value
\(P(\theta)\): Our view of the possible values of our model parameters before observing any actual data, i.e. the prior
\(P(y)\): The unconditional probability of the outcome, i.e. the probability after considering all possible parameter. The objective is to estimate the posterior density of parameter \(P(\theta|y,X)\), typically described in terms of a point estimate of the expected value of posterior density (e.g., the median of the distribution) and a measure of the variance. To solve for \(P(\theta|y,X)\) analytically, the data likelihood \(P(y|\theta,X)\), prior \(P(\theta)\), and marginal likelihood of outcome \(P(y)\) are needed, but because there is usually no closed-form analystical solution, Markov chain Monte Carlo methods are used to numerically solve the posterior density by directly generating random draws of parameters. The general steps in a Bayesian analysis followed included: Specify a joint distribution for the outcome(s) and all the unknowns, which proportional to a posterior distribution of the unknowns conditional on the observed data
Use Markov Chain Monte Carlo (MCMC) to draw from posterior distribution.
*Evaluation of model fit and revise the model as appropriate.

Bayesian repeated measures analyses were fit separately for combined core and non-core winter range, core winter range (WC) and non-core winter range plots (WNC). For willow height, models were fit using weekly informative priors and a gamma distribution with the “stan_glmer” function in the “rstanarm” package (Goodrich et al. 2020)(Brilleman et al. 2018). Gamma regression is commonly used when the response variable is continuous and positive.

For the simplest case a GLM for a continuous outcome is simply a linear model and the likelihood for one observation is a conditionally normal probability density function.
\[\frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{1}{2} \left(\frac{y - \mu}{\sigma}\right)^2},\] where \(\mu = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) is a linear predictor and \(\sigma\) is the standard deviation of the error in predicting the outcome, \(y\). A linear predictor \(\eta = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) can more generally be related to the conditional mean of the outcome via a link function \(g\). This maps the range of values on which the outcome is defined. For the gamma distribution used in the modeling of willow height, the log link function was used. The likelihood is the product of the likelihood contributions of individual observations.

Bayesian analysis requires specifying prior distributions \(f(\alpha)\) and \(f(\boldsymbol{\beta})\) for the intercept and vector of regression coefficients. Prior distributions were represented by normal distributions with a mean zero and a small standard deviation (scale)). The joint posterior distribution for \(\alpha\) and \(\boldsymbol{\beta}\) is proportional to the product of the priors and the \(N\) likelihood contributions:

\[f\left(\boldsymbol{\beta} | \mathbf{y},\mathbf{X}\right) \propto f\left(\alpha\right) \times \prod_{k=1}^K f\left(\beta_k\right) \times \prod_{i=1}^N {f(y_i|\eta_i)},\]

where \(\mathbf{X}\) is the matrix of predictors and \(\eta\) the linear predictor, i.e., the posterior distribution that stan_glmer draws from using MCMC.

Willow cover data were modeled using the beta distribution as the likelihood for the data, \[ f(y_i | a, b) = \frac{y_i^{(a-1)}(1-y_i)^{(b-1)}}{B(a,b)} \] where \(B(\cdot)\) is the beta function. The shape parameters for the distribution are \(a\) and \(b\) and enter into the model according to the following transformations, \[ a = \mu\cdot\phi \\ b = (1-\mu)\cdot\phi \]

If \(g_1(\cdot)\) is some link function, the specification of the shape parameters, \(\mu = g_1^{-1}(\mathbf{X}\boldsymbol{\beta})\), where \(\boldsymbol{X}\) is a \(N\times K\) dimensional matrix of predictors, and \(\boldsymbol{\beta}\), is a \(K\) dimensional vector of parameters associated with each predictor. \(\phi\) is a scalar parameter. With the shape parameter values included, the likelihood takes the form: \[ f(y_i | \mu, \phi) = \frac{y_i^{(\mu\phi-1)}(1-y_i)^{((1-\mu)\phi-1)}}{B(\mu\phi,(1-\mu)\phi)} \]

Bayesian analysis requires specifying prior distributions \(f(\boldsymbol{\beta})\) and \(f(\phi)\) for the vector of regression coefficients and \(\phi\). When modeling \(\phi\) with a linear predictor a full Bayesian analysis requires specifying the prior distributions \(f(\boldsymbol{\beta})\) and \(f(\boldsymbol{\gamma})\).

With a single set of explanatory variables, the posterior distribution of \(\boldsymbol{\beta}\) and \(\phi\) is proportional to the product of the likelihood contributions, the \(K\) priors on the \(\beta_k\) parameters, and \(\phi\), \[ f(\boldsymbol{\beta},\phi|\mathbf{y},\mathbf{X}) \propto \prod_{i=1}^N f(y_i | a, b) \times \prod_{k=1}^K f(\beta_k) \times f(\phi) \]

With two sets of explanatory variables, the posterior distribution of \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) is proportional to the product of the likelihood contribution, the \(K\) priors on the \(\beta_k\) parameters, and the \(J\) priors on the \(\gamma_j\) parameters,

\[ f(\boldsymbol{\beta},\boldsymbol{\gamma}|\mathbf{y},\mathbf{X}) \propto \prod_{i=1}^N f(y_i | a, b) \times \prod_{k=1}^K f(\beta_k) \times \prod_{j=1}^J f(\gamma_j) \]

#### Data import and munging
# read in willow. These data are the output of extensive munging accomplished using code in "EVMP_munging_20210131.Rmd"

willow <- read_csv("./output/exported_data/willow_mcro.csv")

## create timeClass var
willow.ht <- willow %>% 
  filter(!is.na(PLANT_HT_CM)) %>% 
  mutate(timeClass = case_when(yr == 2008 ~ "BL",
                              yr == 2009 ~ "BL", 
                              TRUE ~ as.character(yr))) 

## chr to factor conversions
willow.ht <- willow.ht %>% 
  mutate(timeClass = as.factor(timeClass)) %>% 
  mutate(timeClass = fct_relevel(timeClass, "2018", "2013", "BL")) %>% 
  mutate(FENCED = as.factor(FENCED))


## clean names
willow.ht <- willow.ht %>% 
  janitor::clean_names()

## filter out just the willow species
willow.ht <- willow.ht %>%
  filter(str_detect(species_code, pattern = "SA"))

## filter out 0
willow.ht <- willow.ht %>% 
  filter(plant_ht_cm > 0)

## reorder the time class
willow.ht <- willow.ht %>% 
  mutate(time_class = fct_rev(time_class))

## factor reverse: fenced
willow.ht <- willow.ht %>%
  mutate(fenced = fct_rev(fenced))

## create separate df for WC, WNC, combined WC & WNC, and KV


willow.ht.wc <- willow.ht %>% 
  filter(site_type == "WC") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>% 
  mutate(time_class = fct_drop(time_class))


# non core
willow.ht.nc <- willow.ht %>% 
  # distinct(site_type)
  filter(site_type == "WNC") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>% 
  mutate(time_class = fct_drop(time_class))

# combined WC and WNC
willow.ht.wc.wnc <- willow.ht %>% 
  # distinct(site_type)
  filter(site_type == "WNC" | site_type == "WC") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>% 
  mutate(time_class = fct_drop(time_class))

# kw (note: different munging workflow needed b/c of different sampling yrs)
# willow.ht.wk <- willow %>%
#   clean_names() %>%
#   filter(str_detect(species_code, pattern = "SA")) %>% 
#   filter(site_type == "WK") %>% 
#   mutate(fenced = as.factor(fenced)) %>% 
#   filter(plant_ht_cm > 0) %>% 
#   mutate(time_class = as.factor(time_class)) %>% 
#   mutate(fenced = fct_rev(fenced)) %>% 
#   mutate(time_class = fct_drop(time_class))

#   Making time class == year
willow.ht.wk <- willow %>%
  clean_names() %>%
  filter(str_detect(species_code, pattern = "SA")) %>% 
  filter(site_type == "WK") %>% 
  mutate(fenced = as.factor(fenced)) %>% 
  filter(plant_ht_cm > 0) %>% 
  mutate(time_class = as.factor(yr)) %>% 
  mutate(fenced = fct_rev(fenced)) %>% 
  mutate(time_class = fct_drop(time_class))

Results

Willow Height - Combined Winter Range Macroplot Data

Histograms of Willow Height

## histogram
willow.ht.wc.wnc %>% 
  ggplot() +
  geom_histogram(aes(plant_ht_cm, fill = fenced), color = "grey50") +
  # viridis::scale_fill_viridis(discrete=TRUE) +
  scale_fill_manual(values = colfunc3(2)) +
  labs(x="Plant height (cm)", y = "Count", fill = "", caption = "Core and noncore winter range plots. Macroplot heights") +
  theme_minimal() +
  facet_wrap(~time_class) +
  coord_flip()

ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_wcwnc.png", dpi = 300, width = 6.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_wcwnc.pdf", width = 6.75, height = 3.75)

Modeling and Posterior Description

## run with gamma
wc.wnc.stmod_TCxF <- stan_glmer(plant_ht_cm ~ time_class*fenced +  (1 | site_id), data = willow.ht.wc.wnc,
                      family=Gamma(link="log"),
                      iter = 8000,
                      )

Prior summary: plant_ht_cm ~ time_class*fenced + (1 | site_id)

# Get the prior summary
prior_summary(wc.wnc.stmod_TCxF)

## Priors for model 'wc.wnc.stmod_TCxF' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.70,5.00,5.75,...])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
ppc.plot(wc.wnc.stmod_TCxF) +
  labs(caption = "WC+WNC: willow_ht~TCxF")

ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wcwnc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wcwnc.pdf", width = 4.75, height = 3.75)

# Extract posteriors
posteriors.wc.wnc.ht <- insight::get_parameters(wc.wnc.stmod_TCxF)

## tidy
posteriors.wc.wnc.ht.tidy <- posteriors.wc.wnc.ht %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

# posteriors.wc.wnc.ht.tidy %>%
#   ggplot() +
#   geom_density(aes(fill = parameter, x = values),color="grey50", alpha = .25) +
#   # scale_fill_viridis(discrete = TRUE) +
#   scale_fill_manual(values = colfunc2(10)) +
#   # scale_color_manual(values = colfunc2(ncol(posteriors.wc.wnc.ht.tidy))) +
#   theme_minimal()

Effect Description

# use function defined above
fun.sexit.gt(wc.wnc.stmod_TCxF)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	4.60	0.9	4.47	4.73	1	1	1.00
time_class2013	0.14	0.9	0.09	0.18	1	1	0.00
time_class2018	0.22	0.9	0.18	0.26	1	1	0.00
fencedFenced	-0.57	0.9	-0.79	-0.34	1	1	0.97
time_class2013:fencedFenced	0.34	0.9	0.25	0.43	1	1	0.75
time_class2018:fencedFenced	0.72	0.9	0.63	0.80	1	1	1.00

## save as RTF
fun.sexit.gt(wc.wnc.stmod_TCxF) %>%
 tab_header(title = "Combined WC and WNC", subtitle = "Willow Height Posterior Description") %>%
 gt::gtsave(filename = "./output/tables/willow_height/postdesc_willowht_TCxF_wcwnc.rtf")

Probability of Direction

pd.wc.wnc <- p_direction(wc.wnc.stmod_TCxF)
# Visualize the pd
plot(pd.wc.wnc) +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow height, Combined Core and Non-core Winter Range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wcwnc.png", width = 6.75, height =3.75, dpi = 300)
ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wcwnc.png", width = 6.75, height =3.75)

Region of Practical Equivalence (ROPE)

rope.wcwnc <- rope(wc.wnc.stmod_TCxF, ci=0.9)

rope.wcwnc.gt <- rope.wcwnc %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 2,
    suffixing = TRUE
  ) 
rope.wcwnc.gt

Percent in ROPE
Parameter	CI	ROPE_low	ROPE_high	ROPE_Percentage	Effects	Component
(Intercept)	0.9	-0.1	0.1	0.00	fixed	conditional
time_class2013	0.9	-0.1	0.1	0.04	fixed	conditional
time_class2018	0.9	-0.1	0.1	0.00	fixed	conditional
fencedFenced	0.9	-0.1	0.1	0.00	fixed	conditional
time_class2013:fencedFenced	0.9	-0.1	0.1	0.00	fixed	conditional
time_class2018:fencedFenced	0.9	-0.1	0.1	0.00	fixed	conditional

# rope.wcwnc.gt %>%
#   gt::gtsave(filename = "./output/tables/willow_height/rope_willowht_wcwnc.rtf")

Contrasts

# Results are given on the log (not the response) scale.
# main effects
## main effect: time class
wc.wnc.emm.tc <- emmeans(wc.wnc.stmod_TCxF, ~ time_class)
wc.wnc.pairs.timeClass <- pairs(wc.wnc.emm.tc)

plot(wc.wnc.pairs.timeClass) +
  theme_minimal() +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Core and Noncore Winter Range, Willow Height")

## main effect: fencing
wc.wnc.emm.fenced <- emmeans(wc.wnc.stmod_TCxF, ~ fenced)
wc.wnc.pairs.fenced <- pairs(wc.wnc.emm.fenced)

# plot(wc.wnc.pairs.fenced) +
#   theme_minimal() +
#   labs(subtitle = "Point estimate displayed: median 
# HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Core and Noncore Winter Range, Willow Height")

## the above is not used b/c sig of interactions

## Effect: TC x fencing
emmeans(wc.wnc.stmod_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse=TRUE)

## fenced = Unfenced:
##  contrast    estimate lower.HPD upper.HPD
##  2013 - BL     0.1363    0.0825     0.187
##  2018 - BL     0.2163    0.1693     0.262
##  2018 - 2013   0.0801    0.0356     0.125
## 
## fenced = Fenced:
##  contrast    estimate lower.HPD upper.HPD
##  2013 - BL     0.4726    0.3799     0.573
##  2018 - BL     0.9362    0.8501     1.023
##  2018 - 2013   0.4632    0.3863     0.539
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95

# emmeans(wc.wnc.stmod_TCxF, ~ time_class | fenced) %>% 
#   pairs(reverse=TRUE) %>% 
#   plot(type = "response") +
#   theme_minimal() +
#   geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
#   labs(subtitle = "Point estimate displayed: median 
# HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Core and Noncore Winter Range, Willow Height")
# 
# # ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wcwnc.pdf", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wcwnc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts
emmeans(wc.wnc.stmod_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>%
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log scale")

Pairwise contrasts
Results are given on the log scale
contrast	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unfenced	0.14	0.08	0.19
2018 - BL	Unfenced	0.22	0.17	0.26
2018 - 2013	Unfenced	0.08	0.04	0.13
2013 - BL	Fenced	0.47	0.38	0.57
2018 - BL	Fenced	0.94	0.85	1.02
2018 - 2013	Fenced	0.46	0.39	0.54

# table of contrasts. Results are given on the response scale.
emmeans(wc.wnc.stmod_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")

Pairwise contrasts
Results are given on the response scale.
contrast	fenced	ratio	lower.HPD	upper.HPD
2013 / BL	Unfenced	1.15	1.09	1.21
2018 / BL	Unfenced	1.24	1.18	1.30
2018 / 2013	Unfenced	1.08	1.04	1.13
2013 / BL	Fenced	1.60	1.46	1.77
2018 / BL	Fenced	2.55	2.33	2.77
2018 / 2013	Fenced	1.59	1.47	1.71

emm.func.tc.x.fenced(wc.wnc.stmod_TCxF) +
  xlim(-0.5, 1.25) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Core and Noncore Winter Range, Willow Height \n emm_willowht_TCxF_wcwnc.png")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wcwnc.png", width = 6.5, height = 3.55, dpi=300)

color_scheme_set("teal")

plot.mcmcareas_ht.wcwnc <- mcmc_areas(posteriors.wc.wnc.ht,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ht.wcwnc

ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wcwnc.png", dpi = 300, width = 6.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wcnc.pdf", width = 6.75, height = 3.75)

Effects of Burning

#### Fenced and burned

stmod_ht.wc.wnc.TCxFxB <- stan_glmer(plant_ht_cm ~ time_class*fenced*burned + (1 | site_id), 
                      data = willow.ht.wc.wnc,
                      family=Gamma(link="log"),
                      iter = 8000,
                      seed = 1234
                      )

Prior summary height ~ time_class * burned * fenced + (1 | site_id)

prior_summary(stmod_ht.wc.wnc.TCxFxB)

## Priors for model 'stmod_ht.wc.wnc.TCxFxB' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.70,5.00,5.75,...])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

pp_check(stmod_ht.wc.wnc.TCxFxB) +
  labs(x = "Value", y = "Density") +
  theme_minimal()

# 
ggsave("./output/figures_202108/willow_cover_stats/wc_wnc_willowcover_ppc_TCxFxB.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202108/willow_cover_stats/wc_wnc_willowcover_ppc_TCxFxB.pdf", width = 4.75, height = 3.75)

# Extract posteriors
posteriors.wc.wnc.ht <- insight::get_parameters(stmod_ht.wc.wnc.TCxFxB)

##
color_scheme_set("teal")
# color_scheme_set("teal")

plot.mcmcareas_cov.wcwnc <- mcmc_areas(stmod_ht.wc.wnc.TCxFxB,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  labs(caption="mcmcarea_willowcov_TC_wcwnc.png")

plot.mcmcareas_cov.wcwnc

ggsave("./output/figures_202108/willow_cover_stats/mcmcarea_wilwht_TCxFxB_wcwnc.png", dpi = 300, width = 6.75, height = 11.75)

Effect Description

fun.sexit.gt(stmod_ht.wc.wnc.TCxFxB$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	4.26	0.9	3.96	4.59	1.00	1.00	1.00
time_class2013	-0.41	0.9	-0.64	-0.18	1.00	0.99	0.78
time_class2018	-0.18	0.9	-0.40	0.04	0.90	0.83	0.19
fencedFenced	-0.47	0.9	-0.88	-0.06	0.97	0.95	0.76
burnedUnburned	0.43	0.9	0.09	0.77	0.98	0.97	0.75
time_class2013:fencedFenced	0.70	0.9	0.43	0.99	1.00	1.00	0.99
time_class2018:fencedFenced	1.20	0.9	0.92	1.46	1.00	1.00	1.00
time_class2013:burnedUnburned	0.56	0.9	0.33	0.80	1.00	1.00	0.97
time_class2018:burnedUnburned	0.41	0.9	0.17	0.63	1.00	1.00	0.78
fencedFenced:burnedUnburned	0.02	0.9	-0.47	0.51	0.53	0.46	0.18
time_class2013:fencedFenced:burnedUnburned	-0.26	0.9	-0.56	0.03	0.92	0.88	0.41
time_class2018:fencedFenced:burnedUnburned	-0.56	0.9	-0.85	-0.28	1.00	1.00	0.93
shape	4.11	0.9	3.96	4.26	1.00	1.00	1.00
mean_PPD	109.33	0.9	106.94	111.61	1.00	1.00	1.00
log-posterior	-20389.94	0.9	-20406.75	-20373.21	1.00	1.00	1.00

## Effect: TC x fencing
wcwnc.emm.TCxFxB <- emmeans(stmod_ht.wc.wnc.TCxFxB, ~ time_class | fenced | burned)
# wcwnc.emm.TCxFxB <- pairs(wcwnc.emm.TCxFxB, reverse = TRUE)

wcwnc.emm.TCxFxB

## fenced = Unfenced, burned = Burned:
##  time_class emmean lower.HPD upper.HPD
##  BL           4.26      3.87      4.64
##  2013         3.86      3.50      4.24
##  2018         4.08      3.72      4.46
## 
## fenced = Fenced, burned = Burned:
##  time_class emmean lower.HPD upper.HPD
##  BL           3.79      3.48      4.12
##  2013         4.09      3.78      4.40
##  2018         4.81      4.51      5.11
## 
## fenced = Unfenced, burned = Unburned:
##  time_class emmean lower.HPD upper.HPD
##  BL           4.70      4.56      4.85
##  2013         4.86      4.71      4.99
##  2018         4.93      4.79      5.07
## 
## fenced = Fenced, burned = Unburned:
##  time_class emmean lower.HPD upper.HPD
##  BL           4.25      3.95      4.53
##  2013         4.85      4.55      5.13
##  2018         5.12      4.83      5.40
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95

emmeans(stmod_ht.wc.wnc.TCxFxB, ~ time_class | fenced | burned) %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") + 
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Winter Range, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxFxB_wcwnc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts
emmeans(wcwnc.emm.TCxFxB, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>%
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log scale")

Pairwise contrasts
Results are given on the log scale
contrast	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unfenced	-0.13	-0.26	0.02
2018 - BL	Unfenced	0.02	-0.11	0.16
2018 - 2013	Unfenced	0.15	0.04	0.26
2013 - BL	Fenced	0.45	0.34	0.55
2018 - BL	Fenced	0.94	0.84	1.04
2018 - 2013	Fenced	0.49	0.42	0.57

# table of contrasts. Results are given on the response scale.
emmeans(wcwnc.emm.TCxFxB, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")

Pairwise contrasts
Results are given on the response scale.
contrast	fenced	ratio	lower.HPD	upper.HPD
2013 / BL	Unfenced	0.88	0.76	1.01
2018 / BL	Unfenced	1.02	0.90	1.17
2018 / 2013	Unfenced	1.16	1.03	1.29
2013 / BL	Fenced	1.56	1.40	1.74
2018 / BL	Fenced	2.56	2.30	2.81
2018 / 2013	Fenced	1.64	1.52	1.77

Willow Height - Core Winter Range Macroplot Data

Histograms of Willow Height

willow.ht.wc %>% 
  mutate(time_class = fct_rev(time_class)) %>% 
  ggplot() +
  # geom_density(aes(plant_ht_cm, fill = fenced), color = "black") +
  geom_histogram(aes(plant_ht_cm, fill = fenced), color = "grey50") +
  # viridis::scale_fill_viridis(discrete=TRUE) +
  scale_fill_manual(values = colfunc3(2)) +
  labs(x="Plant height (cm)", y = "Count", fill = "", caption = "Core winter range plots. Macroplot heights") +
  theme_minimal() +
  facet_wrap(~time_class) +
  coord_flip()

ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TCxF_wc.png", dpi = 300, width = 6.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TCxF_wc.pdf", width = 6.75, height = 3.75)

Modeling and Posterior Description

#### null model
# Not currently run; run in past as model comparison  
# Previously run for comparisons but not otherwise reported
stmod_ht.0 <- stan_glmer(plant_ht_cm ~ 1 + (1 | site_id), data = willow.ht.wc,
                      family=Gamma(link="log"),
                      iter = 8000,
                      seed = 1234
                      )

prior_summary(stmod_ht.0)
pp_check(stmod_ht.0)

#### main factors only
## Run, but not used for reporting after model comparisons
stmod_ht1 <- stan_glmer(plant_ht_cm ~ time_class + fenced + (1 | site_id), data = willow.ht.wc,
                      family=Gamma(link="log"),
                      iter = 8000,
                      seed = 1234
                      )

prior_summary(stmod_ht1)
pp_check(stmod_ht1)

#### main factors + interactions
# After model comparisons, this is the selected model reported
stmod_ht2 <- stan_glmer(plant_ht_cm ~ time_class*fenced + (1 | site_id), data = willow.ht.wc,
                      family=Gamma(link="log"),
                      iter = 8000,
                      seed = 1234
                      )

## extract posterior
posteriors.wc.ht <- insight::get_parameters(stmod_ht2)

Prior summary: plant_ht_cm ~ time_class*fenced + (1 | site_id)

# summary(stmod_ht2)
prior_summary(stmod_ht2)

## Priors for model 'stmod_ht2' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.65,5.01,5.06,...])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model comparison.
## Run originally and code retained. stmod_ht2 selected
loo1 <- loo(stmod_ht1,
            k_threshold = 0.7) 
loo2 <- loo(stmod_ht2,
            k_threshold = 0.7) 
comp <- loo_compare(loo2, loo1)
print(comp, simplify = FALSE, digits = 2)

### model checks
ppc.plot(stmod_ht2) +
  labs(x = "Value", y = "Density") +
  theme_minimal()

ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wc.png", width = 4.75, height = 3.75, dpi = 300)
ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wc.pdf", width = 4.75, height = 3.75)

# full posterior description
fun.desc.post1.gt.rope(stmod_ht2) %>% 
  tab_header(title = "Core Winter Range", subtitle = "Willow Height Posterior Description")

Core Winter Range
Willow Height Posterior Description
Parameter	Component	Median	MAD	CI	CI_low	CI_high	pd	ps	ROPE_CI	ROPE_low	ROPE_high	ROPE_Percentage	Rhat	ESS
(Intercept)	conditional	4.4	0.1	0.9	4.3	4.6	1.0	1.0	0.9	-0.1	0.1	0.0	1	782.7
time_class2013	conditional	0.0	0.0	0.9	-0.1	0.1	0.5	0.0	0.9	-0.1	0.1	1.0	1	5622.0
time_class2018	conditional	0.1	0.0	0.9	0.0	0.2	1.0	0.5	0.9	-0.1	0.1	0.5	1	5660.1
fencedFenced	conditional	-0.4	0.1	0.9	-0.6	-0.2	1.0	1.0	0.9	-0.1	0.1	0.0	1	820.8
time_class2013:fencedFenced	conditional	0.5	0.1	0.9	0.4	0.6	1.0	1.0	0.9	-0.1	0.1	0.0	1	5579.5
time_class2018:fencedFenced	conditional	0.8	0.1	0.9	0.7	0.9	1.0	1.0	0.9	-0.1	0.1	0.0	1	5535.1
shape	distributional	4.6	0.1	0.9	4.4	4.9	1.0	1.0	0.9	-0.1	0.1	0.0	NA	NA

#### PP check plot
ppc.plot(stmod_ht2) +
  labs(caption = "WC: willow_ht~TCxF")

# ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wc.pdf", width = 4.75, height = 3.75)

Effect Description

fun.sexit.gt(stmod_ht2)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	4.43	0.9	4.27	4.58	1.00	1.00	1.00
time_class2013	0.00	0.9	-0.06	0.07	0.53	0.13	0.00
time_class2018	0.10	0.9	0.04	0.17	1.00	0.92	0.00
fencedFenced	-0.40	0.9	-0.63	-0.16	1.00	0.99	0.76
time_class2013:fencedFenced	0.47	0.9	0.37	0.57	1.00	1.00	1.00
time_class2018:fencedFenced	0.83	0.9	0.74	0.92	1.00	1.00	1.00

## save as RTF
fun.sexit.gt(stmod_ht2) %>% 
  tab_header(title = "Core Winter Range", subtitle = "Willow Height Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_height/xpostdesc_willowht_TCxF_wc.rtf")

Probability of Direction

pd.wc <- p_direction(stmod_ht2)

# Visualize the pd
plot(pd.wc) +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow height, core winter range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wc.pdf", width = 6.75, height =3.75)

Region of Practical Equivalence (ROPE)

# Region Of Practical Equivalence (ROPE)
rope.wc <- rope(stmod_ht2, ci=0.9)

rope.wc.gt <- rope.wc %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 
rope.wc.gt

Percent in ROPE
Parameter	CI	ROPE_low	ROPE_high	ROPE_Percentage	Effects	Component
(Intercept)	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2013	0.9	-0.1	0.1	1.0	fixed	conditional
time_class2018	0.9	-0.1	0.1	0.5	fixed	conditional
fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2013:fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2018:fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional

# rope.wc.gt %>%
#   gt::gtsave(filename = "./output/tables/willow_height/rope_willowht_wc.rtf")

Contrasts

Results are given on the response scale.

## Effect: TC x fencing
wc.emm.TCxF <- emmeans(stmod_ht2, ~ time_class | fenced)
wc.pairs.TCxF <- pairs(wc.emm.TCxF, reverse = TRUE)

emmeans(stmod_ht2, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") + 
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wc.pdf", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wc.png", width = 4.5, height = 3.55) # save plot

## alt per hanem
emm.wc.ht <- emmeans(stmod_ht2, ~ fenced | time_class) %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") + 
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Height")

# emm.wc.ht
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wc_new.pdf", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TCxF_wc_new.png", width = 4.5, height = 3.55) # save plot

# table of contrasts
emmeans(stmod_ht2, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>%
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log scale")

Pairwise contrasts
Results are given on the log scale
contrast	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unfenced	0.00	-0.07	0.08
2018 - BL	Unfenced	0.10	0.03	0.18
2018 - 2013	Unfenced	0.10	0.04	0.16
2013 - BL	Fenced	0.47	0.38	0.56
2018 - BL	Fenced	0.94	0.86	1.02
2018 - 2013	Fenced	0.46	0.39	0.53

# table of contrasts. Results are given on the response scale.
emmeans(stmod_ht2, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")

Pairwise contrasts
Results are given on the response scale.
contrast	fenced	ratio	lower.HPD	upper.HPD
2013 / BL	Unfenced	1.00	0.93	1.08
2018 / BL	Unfenced	1.11	1.03	1.19
2018 / 2013	Unfenced	1.11	1.04	1.17
2013 / BL	Fenced	1.60	1.46	1.75
2018 / BL	Fenced	2.55	2.35	2.76
2018 / 2013	Fenced	1.59	1.47	1.70

color_scheme_set("teal")

plot.mcmcareas_ht.wc <- mcmc_areas(posteriors.wc.ht,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ht.wc

ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wc.png", dpi = 300, width = 6.75, height = 3.75)
ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wc.pdf", width = 6.75, height = 3.75)

bayesplot::mcmc_areas(posteriors.wc.ht,
           pars = c("time_class2013", "time_class2018", "time_class2013:fencedFenced","time_class2018:fencedFenced"),
           prob = 0.9) +
  theme_minimal() +
  labs(title = "Posterior distributions", subtitle = "with medians and 90% intervals", caption = "Willow height, WC")

ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wc.png", dpi = 300, width = 6.75, height = 3.75)
ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wc.pdf", width = 6.75, height = 3.75)

Willow Height - Noncore Winter Range Macroplot Data

Histograms of Willow Height

## histogram
willow.ht.nc %>% 
  ggplot() +
  geom_histogram(aes(plant_ht_cm, fill = fenced), color = "grey50") +
  # viridis::scale_fill_viridis(discrete=TRUE) +
  scale_fill_manual(values = colfunc3(1)) +
  labs(x="Plant height (cm)", y = "Count", fill = "", caption = "Core winter range plots. Macroplot heights") +
  theme_minimal() +
  facet_wrap(~time_class) +
  coord_flip()

ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_wnc.png", dpi = 300, width = 6.75, height = 3.75)
ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_wnc.pdf", width = 6.75, height = 3.75)

Modeling and Posterior Description

## model
#### main factors only
## run with gamma
nc.stmod_ht1 <- stan_glmer(plant_ht_cm ~ time_class +  (1 | site_id), data = willow.ht.nc,
                       family=Gamma(link="log"),
                      iter = 8000,
                      )
posteriors.nc.ht <- insight::get_parameters(nc.stmod_ht1)

Prior summary: plant_ht_cm ~ time_class + (1 | site_id)

prior_summary(nc.stmod_ht1)

## Priors for model 'nc.stmod_ht1' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0], scale = [2.5,2.5])
##   Adjusted prior:
##     ~ normal(location = [0,0], scale = [5.77,5.02])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
ppc.plot(nc.stmod_ht1) +
  labs(caption = "WNC: willow_ht~TC")

ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wnc.png", dpi = 300, width = 4.75, height = 3.75)

Effect Description

# use function defined above
fun.sexit.gt(nc.stmod_ht1)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	4.85	0.9	4.67	5.03	1	1	1.00
time_class2013	0.22	0.9	0.16	0.28	1	1	0.01
time_class2018	0.28	0.9	0.22	0.33	1	1	0.25

## save as RTF
fun.sexit.gt(nc.stmod_ht1) %>% 
  tab_header(title = "Noncore Winter Range", subtitle = "Willow Height Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_height/xposterior_willowht_wnc.rtf")

# tidy the posteriors
posteriors.nc.ht.tidy <- posteriors.nc.ht %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

Probability of Direction

# Visualize the pd
p_direction(nc.stmod_ht1) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow height, Noncore winter range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wnc.png", width = 6.75, height =3.75, dpi = 300)

Region of Practical Equivalence (ROPE)

rope.wnc <- rope(nc.stmod_ht1, ci=0.9)

rope.wnc.gt <- rope.wnc %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.wnc.gt

Percent in ROPE
Parameter	CI	ROPE_low	ROPE_high	ROPE_Percentage	Effects	Component
(Intercept)	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2013	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2018	0.9	-0.1	0.1	0.0	fixed	conditional

rope.wnc.gt %>%
  gt::gtsave(filename = "./output/tables/willow_height/rope_willowht_wnc.rtf")

Contrasts

Results are given on the the response scale.

### Contrasts - WC time class
nc.emm <- emmeans(nc.stmod_ht1, ~ time_class)

nc.emm %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Non-core Winter Range, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wnc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts
emmeans(nc.stmod_ht1, ~ time_class) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>%
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log scale")

Pairwise contrasts
Results are given on the log scale
contrast	estimate	lower.HPD	upper.HPD
2013 - BL	0.22	0.14	0.29
2018 - BL	0.28	0.21	0.34
2018 - 2013	0.06	-0.01	0.13

# table of contrasts. Results are given on the response scale.
emmeans(nc.stmod_ht1, ~ time_class) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")

Pairwise contrasts
Results are given on the response scale.
contrast	ratio	lower.HPD	upper.HPD
2013 / BL	1.24	1.15	1.34
2018 / BL	1.32	1.24	1.41
2018 / 2013	1.06	0.99	1.14

color_scheme_set("teal")
# color_scheme_set("teal")

plot.mcmcareas_ht.wnc <- mcmc_areas(posteriors.nc.ht,
           prob = 0.9) +
  ggtitle("Noncore Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ht.wnc

ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TCxF_wnc.png", dpi = 300, width = 6.75, height = 3.75)

# bayesplot::mcmc_areas(posteriors.nc.ht,
#            pars = c("time_class2013","time_class2018"),
#            prob = 0.9) +
#   theme_minimal() +
#   labs(title = "Posterior distributions", subtitle = "with medians and 90% intervals", caption = "Willow height, WNC")
# 
# ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TC_wnc.png", dpi = 300, width = 6.75, height = 3.75)
# ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TC_wnc.pdf", width = 6.75, height = 3.75)

# Revised plot 2/15/21
# willow height combined plot per revisions
## main effect: fencing
wc.emm.fenced <- emmeans(stmod_ht2, ~ fenced)
wc.pairs.fenced <- pairs(wc.emm.fenced)

plot(wc.pairs.fenced) +
  theme_minimal() +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Height")

nc.emm %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Non-core Winter Range, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wnc.png", width = 4.5, height = 3.55) # save plot

emm.wc <- plot(wc.pairs.fenced) +
  theme_minimal() +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Height")

# emm.wnc.ht
emm.wnc <- nc.emm %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Non-core Winter Range, Willow Height") +
  xlim(-0.5, 1.5)
emm.wnc

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wnc.png", width = 7.5, height = 3.55) # save plot


# emm.wc.ht + emm.wnc +
#   plot_layout(ncol = 2) +
#   plot_annotation(caption = "emm_willowht_TC_wcwnc.png")

## two panel
xxx.wc <- emm.func.tc.x.fenced(stmod_ht2) +
  labs(x = "Estimate", y = "Contrast", Title = "Core Range")


xxx.wnc <- nc.emm %>% 
  pairs(reverse = TRUE) %>% 
  plot(type="response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Noncore Range", x = "Estimate", y = "Contrast") +
  xlim(-0.5, 1.5)

xxx.wc + xxx.wnc +
  plot_layout(ncol = 1, heights = c(2,1)) +
  plot_annotation(caption = "Willow height \n emm_willowht_wcwnc_2panel.png")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_wcwnc_2panel.png", width = 4.5, height = 6.55) # save plot

Willow Height - Kawuneeche Valley Macroplot Data

Histograms of Willow Height

willow.ht.wk <- willow.ht.wk %>% 
  mutate(time_class = fct_relevel(time_class, "BL"))

willow.ht.wk %>% 
  ggplot() +
  geom_histogram(aes(plant_ht_cm, fill = fenced), color = "grey50") +
  # viridis::scale_fill_viridis(discrete=TRUE) +
  scale_fill_manual(values = colfunc3(2)) +
  labs(x="Plant height (cm)", y = "Count", fill = "", caption = "KV. Macroplot heights") +
  theme_minimal() +
  facet_wrap(~time_class) +
  coord_flip() +
  theme(legend.position = "none")

ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_WK.png", dpi = 300, width = 6.75, height = 4.75)

willow.ht.wk %>% 
  ggplot() +
  geom_histogram(aes(plant_ht_cm, fill = fenced), color = "grey50") +
  # viridis::scale_fill_viridis(discrete=TRUE) +
  # scale_fill_manual(values = colfunc3(2)) +
  scale_fill_manual(values=c("grey80", "grey40")) +
  labs(x="Plant height (cm)", y = "Count", fill = "", caption = "KV. Macroplot heights") +
  theme_minimal() +
  facet_grid(time_class ~ fenced) +
  # facet_wrap(~time_class) +
  coord_flip() +
  theme(legend.position = "none")

ggsave("./output/figures_202108/willow_height_stats/hist_willowht_TC_WK_bw.png", dpi = 300, width = 6.75, height = 6.75)

Modeling and Posterior Description

#### TC only (no fenced)
# stmod_ht_wk_TC <- stan_glmer(plant_ht_cm ~ time_class + (1 | site_id), data = willow.ht.wk,
#                       family=Gamma(link="log"),
#                       iter = 8000,
#                       seed = 1234
#                       )
# TCxF need the full data for fencing
stmod_ht_wk_TCxF <- stan_glmer(plant_ht_cm ~ time_class * fenced + (1 | site_id), data = willow.ht.wk,
                      family=Gamma(link="log"),
                      iter = 8000,
                      seed = 1234
                      )

## extract posterior
posteriors.wk.ht <- insight::get_parameters(stmod_ht_wk_TCxF)
# posteriors.wk.ht.tc.f <- insight::get_parameters(stmod_ht_wk_TCxF)

Prior summary: plant_ht_cm ~ time_class * fenced + (1 | site_id)

## prior summary
prior_summary(stmod_ht_wk_TCxF)

## Priors for model 'stmod_ht_wk_TCxF' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [6.16,6.10,5.12,...])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
ppc.plot(stmod_ht_wk_TCxF) +
  labs(caption = "WK: willow_ht~TCxF")

ggsave("./output/figures_202108/willow_height_stats/ppc_willowht_wk.png", dpi = 300, width = 4.75, height = 3.75)

Effect Description

# use function defined above
fun.sexit.gt(stmod_ht_wk_TCxF)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	4.01	0.9	3.74	4.29	1.00	1.00	1.00
time_class2016	-0.14	0.9	-0.25	-0.02	0.98	0.89	0.01
time_class2017	-0.21	0.9	-0.32	-0.09	1.00	0.99	0.09
time_class2018	-0.45	0.9	-0.56	-0.35	1.00	1.00	0.99
fencedFenced	-0.14	0.9	-0.59	0.31	0.71	0.64	0.27
time_class2016:fencedFenced	0.58	0.9	0.40	0.76	1.00	1.00	0.99
time_class2017:fencedFenced	0.76	0.9	0.58	0.94	1.00	1.00	1.00
time_class2018:fencedFenced	1.14	0.9	0.96	1.30	1.00	1.00	1.00

fun.sexit.gt(stmod_ht_wk_TCxF) %>% 
  tab_header(title = "Kawuneeche Valley", subtitle = "Willow Height Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_height/xpostdesc_willowht_TC_wk.rtf")

Probability of Direction

pdkv <- p_direction(stmod_ht_wk_TCxF)
# Visualise the pd
plot(pdkv) +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow height, KV \n pd_willowht_wk.png") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_height_stats/pd_willowht_wk.png", width = 6.75, height =3.75, dpi = 300)

Region of Practical Equivalence (ROPE)

p.in.rope.kv <- rope(stmod_ht_wk_TCxF, ci=0.9)

pd.tab.kv <- p.in.rope.kv %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 
pd.tab.kv

Percent in ROPE
Parameter	CI	ROPE_low	ROPE_high	ROPE_Percentage	Effects	Component
(Intercept)	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2016	0.9	-0.1	0.1	0.3	fixed	conditional
time_class2017	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2018	0.9	-0.1	0.1	0.0	fixed	conditional
fencedFenced	0.9	-0.1	0.1	0.3	fixed	conditional
time_class2016:fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2017:fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional
time_class2018:fencedFenced	0.9	-0.1	0.1	0.0	fixed	conditional

pd.tab.kv %>%
  gt::gtsave(filename = "./output/tables/willow_height/rope_willowht_wk.rtf")

Contrasts

## contrasts
# emmeans(stmod_ht_wk_TCxF, ~ time_class) %>% 
#   pairs() %>% 
#   # plot() + # on the lm predict scale
#   plot(type = "response") +
#   theme_minimal() +
#   geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
#   labs(subtitle = "Point estimate displayed: median 
# HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "KV, Willow Height")
# 
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk.png", width = 4.5, height = 3.55) # save plot

emmeans(stmod_ht_wk_TCxF, ~ fenced | time_class)

## time_class = 2011:
##  fenced   emmean lower.HPD upper.HPD
##  Unfenced   4.01      3.66      4.35
##  Fenced     3.87      3.43      4.29
## 
## time_class = 2016:
##  fenced   emmean lower.HPD upper.HPD
##  Unfenced   3.88      3.54      4.23
##  Fenced     4.31      3.87      4.73
## 
## time_class = 2017:
##  fenced   emmean lower.HPD upper.HPD
##  Unfenced   3.80      3.47      4.15
##  Fenced     4.42      3.98      4.84
## 
## time_class = 2018:
##  fenced   emmean lower.HPD upper.HPD
##  Unfenced   3.56      3.22      3.89
##  Fenced     4.55      4.11      4.95
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95

emmeans(stmod_ht_wk_TCxF, ~ fenced | time_class) %>% 
  pairs(reverse=TRUE) %>% 
  # plot() + # on the lm predict scale
  plot(type = "response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  xlim(-1, 5) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "KV, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk.png", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk.pdf", width = 6.5, height = 4.55) # save plot

emmeans(stmod_ht_wk_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse=TRUE) %>% 
  # plot() + # on the lm predict scale
  plot(type = "response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  xlim(-2, 5) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "KV, Willow Height")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk_alt.png", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk_alt.pdf", width = 6.5, height = 4.55) # save plot

## custom
emmeans(stmod_ht_wk_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  ggplot() +
  geom_linerange(aes(y=contrast, xmin=lower.HPD, xmax=upper.HPD), color="#babefa", size=2.5) +
  geom_point(aes(x=estimate, y=contrast)) +
  theme_minimal() +
  geom_vline(xintercept = 0, lty="dashed", size=1) +
  facet_wrap(~fenced, ncol=1) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "KV, Willow Height")

  # labs(x= "Estimate", y="Contrast", caption = "KV")

ggsave("./output/figures_202108/willow_height_stats/emm_willowht_TC_wk_revised.png", width = 4.5, height = 3.75) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_ht_wk_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log (not the response) scale.")

Pairwise contrasts
Results are given on the log (not the response) scale.
contrast	fenced	estimate	lower.HPD	upper.HPD
2016 - 2011	Unfenced	-0.14	-0.28	-0.01
2017 - 2011	Unfenced	-0.21	-0.34	-0.08
2017 - 2016	Unfenced	-0.07	-0.20	0.05
2018 - 2011	Unfenced	-0.45	-0.58	-0.32
2018 - 2016	Unfenced	-0.31	-0.43	-0.19
2018 - 2017	Unfenced	-0.24	-0.36	-0.13
2016 - 2011	Fenced	0.44	0.27	0.61
2017 - 2011	Fenced	0.55	0.38	0.73
2017 - 2016	Fenced	0.11	-0.04	0.27
2018 - 2011	Fenced	0.69	0.52	0.84
2018 - 2016	Fenced	0.24	0.10	0.39
2018 - 2017	Fenced	0.13	-0.01	0.28

# table of contrasts. Results are given on the response scale.
emmeans(stmod_ht_wk_TCxF, ~ time_class) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")

Pairwise contrasts
Results are given on the response scale.
contrast	ratio	lower.HPD	upper.HPD
2016 / 2011	1.16	1.04	1.29
2017 / 2011	1.19	1.06	1.32
2017 / 2016	1.02	0.92	1.13
2018 / 2011	1.12	1.01	1.24
2018 / 2016	0.97	0.88	1.06
2018 / 2017	0.95	0.86	1.04

# EMM summaries with type = "response", the tests and confidence intervals are done before back-transforming. The ratios estimated here are actually ratios of geometric means. In general, a model with a log response is in fact a model for relative effects of any of its linear predictors, and this back-transformation to ratios goes hand-in-hand with that.

color_scheme_set("teal")
# color_scheme_set("teal")

plot.mcmcareas_ht.wk <- mcmc_areas(posteriors.wk.ht,
           prob = 0.9) +
  ggtitle("Kawuneeche Valley") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ht.wk

ggsave("./output/figures_202108/willow_height_stats/mcmcarea_willowht_TC_wk.png", dpi = 300, width = 6.75, height = 3.75)

Willow Cover - Combined Winter Range Macroplot Data

# read in data
# canopy.all.willow.sum <- read_csv("./output/exported_data/willow_cover_20200725.csv") %>% 
#   janitor::clean_names()

canopy.all.willow.sum <- read_csv("./output/exported_data/willow_macroplot_cover_BLto2018b.csv") %>% 
  janitor::clean_names()

# munge
canopy.all.willow.sum <- canopy.all.willow.sum %>% 
  mutate(burned = case_when(burned == "Y" ~ "Burned",
                            burned == "N" ~ "Unburned",
                            is.na(burned) ~ "Unburned",
                            TRUE ~ burned)) %>% 
  mutate(fenced = case_when(site_id == "WK01" ~ "Fenced",
                            site_id == "WK02" ~ "Fenced",
                            site_id == "WK03" ~ "Unfenced",
                            site_id == "WK04" ~ "Unfenced",
                            site_id == "WK05" ~ "Fenced",
                            site_id == "WK06" ~ "Unfenced",
                            site_id == "WK07" ~ "Unfenced",
                            site_id == "WK08" ~ "Unfenced",
                            site_id == "WK09" ~ "Unfenced",
                            site_id == "WK10" ~ "Fenced",
                            TRUE ~ fenced)) %>% 
  mutate(time_class = as_factor(time_class)) %>%
  mutate(fenced = as_factor(fenced)) %>%
  mutate(burned = as_factor(burned)) %>% 
  mutate(time_class = fct_rev(time_class))

# canopy.all.willow.sum %>% 
#   tabyl(site_type,time_class)

### create subsets of data for modeling
# pull out WK, note different years of data
cov.wk <- canopy.all.willow.sum %>% 
  filter(site_type == "WK") %>% 
  filter(!is.na(cover_allwillow)) %>% 
  mutate(cover_allwillow = cover_allwillow/100)

# pull out the core plus noncore winter range plots
cov.wc.wnc <- canopy.all.willow.sum %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  filter(site_type == "WNC" | site_type == "WC") %>% 
  mutate(cover_allwillow = cover_allwillow/100) %>% 
  filter(!is.na(cover_allwillow))

# pull out the core winter range plots
cov.wc <- canopy.all.willow.sum %>% 
  filter(site_type == "WC") %>%
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  mutate(cover_allwillow = cover_allwillow/100) %>%
  filter(!is.na(cover_allwillow))

# pull out the noncore winter range plots
cov.wnc <- canopy.all.willow.sum %>% 
  filter(site_type == "WNC") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  mutate(cover_allwillow = cover_allwillow/100) %>% 
  filter(!is.na(cover_allwillow))

Summary Tables

#### summary tables: time class and burned
cov.wc.wnc %>% 
  group_by(time_class, fenced) %>% 
  summarytools::descr(cover_allwillow,stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "Combined WC and WNC macroplot data: willow cover - fenced + time_class")

Combined WC and WNC macroplot data: willow cover - fenced + time_class
time_class	fenced	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Fenced	cover_allwillow	0.17	0.24	0.00	0.11	1	17	100
BL	Unfenced	cover_allwillow	0.33	0.27	0.00	0.27	1	51	100
2018	Fenced	cover_allwillow	0.55	0.35	0.05	0.56	1	28	100
2018	Unfenced	cover_allwillow	0.45	0.35	0.00	0.41	1	66	100
2013	Fenced	cover_allwillow	0.29	0.30	0.00	0.18	1	28	100
2013	Unfenced	cover_allwillow	0.31	0.30	0.00	0.25	1	60	100

## time class, fenced and burned
cov.wc.wnc %>% 
  group_by(time_class, fenced, burned) %>% 
  summarytools::descr(cover_allwillow, stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "Combined WC and WNC macroplot data: willow cover ~ burned + fenced + time_class")

Combined WC and WNC macroplot data: willow cover ~ burned + fenced + time_class
time_class	fenced	burned	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Fenced	Unburned	cover_allwillow	0.27	0.27	0.01	0.20	1.00	10	100
BL	Fenced	Burned	cover_allwillow	0.03	0.03	0.00	0.02	0.09	7	100
BL	Unfenced	Unburned	cover_allwillow	0.34	0.27	0.00	0.27	1.00	50	100
BL	Unfenced	Burned	cover_allwillow	0.03	NA	0.03	0.03	0.03	1	100
2018	Fenced	Unburned	cover_allwillow	0.69	0.33	0.14	0.76	1.00	15	100
2018	Fenced	Burned	cover_allwillow	0.39	0.29	0.05	0.26	1.00	13	100
2018	Unfenced	Unburned	cover_allwillow	0.49	0.33	0.00	0.47	1.00	58	100
2018	Unfenced	Burned	cover_allwillow	0.15	0.34	0.00	0.03	1.00	8	100
2013	Fenced	Unburned	cover_allwillow	0.43	0.30	0.02	0.48	1.00	15	100
2013	Fenced	Burned	cover_allwillow	0.12	0.21	0.00	0.04	0.79	13	100
2013	Unfenced	Unburned	cover_allwillow	0.34	0.29	0.00	0.27	1.00	51	100
2013	Unfenced	Burned	cover_allwillow	0.11	0.23	0.01	0.02	0.72	9	100

cov.wk %>% 
  group_by(time_class, fenced) %>% 
  summarytools::descr(cover_allwillow, stats = "common", na.rm=TRUE) %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "WK macroplot data: willow cover ~ fenced * time_class")

WK macroplot data: willow cover ~ fenced * time_class
time_class	fenced	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Fenced	cover_allwillow	0.00	0.00	0.00	0.00	0.00	3	100
BL	Unfenced	cover_allwillow	0.00	0.00	0.00	0.00	0.00	5	100
2018	Fenced	cover_allwillow	0.46	0.27	0.14	0.48	0.75	4	100
2018	Unfenced	cover_allwillow	0.26	0.25	0.02	0.17	0.53	5	100
2017	Fenced	cover_allwillow	0.34	0.25	0.17	0.23	0.62	3	100
2017	Unfenced	cover_allwillow	0.23	0.24	0.02	0.16	0.58	4	100
2016	Fenced	cover_allwillow	0.31	0.20	0.19	0.20	0.54	3	100
2016	Unfenced	cover_allwillow	0.23	0.23	0.04	0.18	0.54	4	100

Histograms of Willow Cover

cov.wc.wnc %>% 
  ggplot(aes(cover_allwillow)) +
  geom_histogram(color="grey50") +
  facet_grid(time_class~site_type) +
  theme_minimal() +
  labs(x="Cover", y = "Count")

ggsave("./output/figures_202108/willow_cover_stats/hist_willowht_TC_wcwnc.png", dpi = 300, width = 6.75, height = 3.75)

Modeling and Posterior Description

# cover_allwillow \~ time_class \* fenced + (1 \| site_id)
## prep for modeling
## subtract from max to conform
#### Model
cov.wc.wnc <- cov.wc.wnc %>% 
  mutate(cover_allwillow = case_when(cover_allwillow == 1 ~ 0.999999,
                                     TRUE ~ cover_allwillow))

## model
stmod_cov.wc.wnc <- stan_glmer(cover_allwillow ~ time_class * fenced + (1 | site_id), 
                      data = cov.wc.wnc,
                      family = mgcv::betar,
                      prior_aux = exponential(1),
                      iter = 8000,
                      seed = 1234)

Prior summary: cover_allwillow ~ time_class * fenced + (1 | site_id)

prior_summary(stmod_cov.wc.wnc)

## Priors for model 'stmod_cov.wc.wnc' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.15,5.22,5.49,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
# prior_summary(stmod_cov.wc.wnc)

pp_check(stmod_cov.wc.wnc) +
  labs(x = "Value", y = "Density") +
  theme_minimal() +
  labs(caption = "WC & WNC: cover_allwillow ~ time_class * fenced + (1 | site_id)")

ggsave("./output/figures_202108/willow_cover_stats/ppc_willowcov_wcwnc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202108/willow_cover_stats/ppc_willowcov_wcwnc.pdf", width = 4.75, height = 3.75)

## extract posterior
posteriors.cov.wcwnc <- insight::get_parameters(stmod_cov.wc.wnc)

Effect Description

fun.sexit.gt(stmod_cov.wc.wnc$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-1.69	0.9	-2.46	-0.97	1.00	1.00	1.00
time_class2018	2.61	0.9	2.03	3.15	1.00	1.00	1.00
time_class2013	0.64	0.9	0.10	1.21	0.97	0.96	0.84
fencedUnfenced	0.88	0.9	0.02	1.75	0.95	0.94	0.87
time_class2018:fencedUnfenced	-1.80	0.9	-2.44	-1.17	1.00	1.00	1.00
time_class2013:fencedUnfenced	-0.29	0.9	-0.93	0.35	0.77	0.73	0.49
(phi)	4.54	0.9	3.64	5.51	1.00	1.00	1.00
mean_PPD	0.41	0.9	0.39	0.44	1.00	1.00	1.00
log-posterior	245.12	0.9	225.27	263.47	1.00	1.00	1.00

fun.sexit.gt(stmod_cov.wc.wnc$stanfit) %>%
  tab_header(title = "Combined Core and Noncore Winter Range", subtitle = "Willow Cover Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_cover/xposterior_willowht_wcwnc.rtf")

color_scheme_set("teal")
# color_scheme_set("teal")

plot.mcmcareas_cov.wcwnc <- mcmc_areas(posteriors.cov.wcwnc,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  labs(caption="mcmcarea_willowcov_TC_wcwnc.png")

plot.mcmcareas_cov.wcwnc

ggsave("./output/figures_202108/willow_cover_stats/mcmcarea_willowcov_TC_wcwnc.png", dpi = 300, width = 6.75, height = 3.75)

Probability of Direction

# Visualize the pd
p_direction(stmod_cov.wc.wnc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow cover, Core and noncore winter range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wcwnc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wcwnc.pdf", width = 6.75, height =3.75)

Region of Practical Equivalence (ROPE)

rope.cov.wcwnc.gt <- rope(stmod_cov.wc.wnc$stanfit, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 2,
    suffixing = TRUE
  ) 

rope.cov.wcwnc.gt %>%
  gt::gtsave(filename = "./output/tables/willow_cover/rope_willowcov_wcwnc.rtf")

Contrasts

# contrast - time class and fencing
emmeans(stmod_cov.wc.wnc, ~ time_class | fenced) %>% 
  pairs(reverse=TRUE) %>% 
  # plot() + # Results are given on the logit (not the response) scale.
  plot(type = "response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Core and Noncore Winter Range, Willow Cover")

ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wcwnc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wcwnc.pdf", width = 6.75, height =3.75)


# table of contrasts
emmeans(stmod_cov.wc.wnc, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log odds ratio, not the response scale")

Pairwise contrasts
Results are given on the log odds ratio, not the response scale
contrast	fenced	estimate	lower.HPD	upper.HPD
2018 - BL	Fenced	2.61	1.91	3.26
2013 - BL	Fenced	0.64	0.00	1.31
2013 - 2018	Fenced	-1.97	-2.49	-1.45
2018 - BL	Unfenced	0.81	0.42	1.19
2013 - BL	Unfenced	0.35	-0.06	0.73
2013 - 2018	Unfenced	-0.46	-0.80	-0.13

emmeans(stmod_cov.wc.wnc, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale")

Pairwise contrasts
Results are given on the response scale
contrast	fenced	odds.ratio	lower.HPD	upper.HPD
2018 / BL	Fenced	13.55	6.14	24.58
2013 / BL	Fenced	1.89	0.89	3.41
2013 / 2018	Fenced	0.14	0.08	0.22
2018 / BL	Unfenced	2.24	1.44	3.17
2013 / BL	Unfenced	1.41	0.88	2.02
2013 / 2018	Unfenced	0.63	0.44	0.86

# bayesplot::mcmc_areas(posteriors.cov.wcwnc,
#            pars = c("time_class2013", "time_class2018", "time_class2013:fencedFenced","time_class2018:fencedFenced"),
#            prob = 0.9) +
#   theme_minimal() +
#   labs(title = "Posterior distributions", subtitle = "with medians and 90% intervals", caption = "Willow height, WC & WNC")
# 
# ggsave("./output/figures_202108/willow_cover_stats/mcmcarea_willowcov_TCxF_wcwnc.png", dpi = 300, width = 6.75, height = 3.75)
# ggsave("./output/figures_202108/willow_cover_stats/mcmcarea_willowcov_TCxF_wcwnc.pdf", width = 6.75, height = 3.75)

Effects of Burning - Combined Range

#### Fenced and burned

stmod_cov.wc.wnc.TCxFxB <- stan_glmer(cover_allwillow ~ time_class * burned * fenced + (1 | site_id), 
                      data = cov.wc.wnc,
                      family = mgcv::betar,
                      prior_aux = exponential(1),
                      iter = 8000,
                      seed = 1234)

Prior summary cover_allwillow ~ time_class * burned * fenced + (1 | site_id)

prior_summary(stmod_cov.wc.wnc.TCxFxB)

## Priors for model 'stmod_cov.wc.wnc.TCxFxB' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.15,5.22,6.19,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

pp_check(stmod_cov.wc.wnc.TCxFxB) +
  labs(x = "Value", y = "Density") +
  theme_minimal()

# 
ggsave("./output/figures_202107/wc_wnc_willowcover_ppc_TCxFxB.png", dpi = 300, width = 4.75, height = 3.75)

## extract posterior
posteriors.cov.wcwnc.TCxFxB <- insight::get_parameters(stmod_cov.wc.wnc.TCxFxB)

Effect Description

fun.sexit.gt(stmod_cov.wc.wnc.TCxFxB$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-1.02	0.9	-1.98	-0.08	0.96	0.95	0.90
time_class2018	3.00	0.9	2.32	3.75	1.00	1.00	1.00
time_class2013	0.78	0.9	0.09	1.47	0.97	0.96	0.87
burnedBurned	-1.36	0.9	-2.79	0.00	0.95	0.94	0.89
fencedUnfenced	0.42	0.9	-0.65	1.44	0.75	0.72	0.57
time_class2018:burnedBurned	-0.91	0.9	-2.04	0.15	0.91	0.90	0.82
time_class2013:burnedBurned	-0.44	0.9	-1.58	0.66	0.75	0.72	0.58
time_class2018:fencedUnfenced	-2.14	0.9	-2.93	-1.38	1.00	1.00	1.00
time_class2013:fencedUnfenced	-0.38	0.9	-1.17	0.36	0.80	0.76	0.57
burnedBurned:fencedUnfenced	0.43	0.9	-2.42	3.37	0.59	0.58	0.53
time_class2018:burnedBurned:fencedUnfenced	-0.10	0.9	-2.78	2.56	0.53	0.51	0.45
time_class2013:burnedBurned:fencedUnfenced	-0.70	0.9	-3.52	2.07	0.66	0.65	0.59
(phi)	4.68	0.9	3.70	5.64	1.00	1.00	1.00
mean_PPD	0.41	0.9	0.39	0.44	1.00	1.00	1.00
log-posterior	239.94	0.9	220.32	258.70	1.00	1.00	1.00

color_scheme_set("teal")
# color_scheme_set("teal")

plot.mcmcareas_cov.wcwnc.TCxFxB <- mcmc_areas(posteriors.cov.wcwnc.TCxFxB,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_cov.wcwnc.TCxFxB

ggsave("./output/figures_202108/willow_cover_stats/mcmcarea_willowcov_wcwnc_TCxFxB.png", dpi = 300, width = 6.75, height = 3.75)

Willow Cover - Core Winter Range Macroplot Data

Summary Tables

#### summary tables: time class and burned
cov.wc %>% 
  group_by(time_class) %>% 
  summarytools::descr(cover_allwillow,stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "WC macroplot data: willow cover ~ WC: all willow cover - time_class")

WC macroplot data: willow cover ~ WC: all willow cover - time_class
time_class	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	cover_allwillow	0.22	0.24	0	0.18	1	36	100
2018	cover_allwillow	0.44	0.35	0	0.33	1	62	100
2013	cover_allwillow	0.24	0.27	0	0.14	1	61	100

## time class, fenced and burned
cov.wc %>% 
  group_by(time_class, burned) %>% 
  summarytools::descr(cover_allwillow, stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "WC macroplot data: All willow cover - burned + time_class")

WC macroplot data: All willow cover - burned + time_class
time_class	burned	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Unburned	cover_allwillow	0.27	0.25	0	0.22	1.00	28	100
BL	Burned	cover_allwillow	0.03	0.03	0	0.02	0.09	8	100
2018	Unburned	cover_allwillow	0.51	0.35	0	0.46	1.00	41	100
2018	Burned	cover_allwillow	0.30	0.33	0	0.20	1.00	21	100
2013	Unburned	cover_allwillow	0.31	0.27	0	0.27	1.00	39	100
2013	Burned	cover_allwillow	0.12	0.22	0	0.03	0.79	22	100

Modeling and Posterior Description

#### Model
## add to min and subtract from max to conform
cov.wc <- cov.wc %>% 
  mutate(cover_allwillow = case_when(cover_allwillow == 1 ~ 0.999999,
                                     TRUE ~ cover_allwillow))


stmod_cov.wc <- stan_glmer(cover_allwillow ~ time_class * fenced + (1 | site_id), 
                      data = cov.wc,
                      family = mgcv::betar,
                      prior_aux = exponential(2),
                      iter = 8000,
                      seed = 1234)

Prior summary: cover_allwillow ~ time_class * fenced + (1 | site_id)

prior_summary(stmod_cov.wc)

## Priors for model 'stmod_cov.wc' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.11,5.12,5.00,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
pp_check(stmod_cov.wc) +
  labs(x = "Value", y = "Density") +
  theme_minimal() +
  labs(caption = "WC: cover_allwillow ~ time_class * fenced + (1 | site_id)")

ggsave("./output/figures_202108/willow_cover_stats/ppc_willowcov_wc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202108/willow_cover_stats/ppc_willowcov_wc.pdf", width = 4.75, height = 3.75)

Effect Description

fun.sexit.gt(stmod_cov.wc$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-1.78	0.9	-2.52	-1.02	1.00	1.00	1.00
time_class2018	2.73	0.9	2.14	3.28	1.00	1.00	1.00
time_class2013	0.68	0.9	0.13	1.22	0.98	0.97	0.88
fencedUnfenced	0.51	0.9	-0.43	1.45	0.81	0.79	0.64
time_class2018:fencedUnfenced	-2.23	0.9	-2.94	-1.51	1.00	1.00	1.00
time_class2013:fencedUnfenced	-0.58	0.9	-1.30	0.14	0.91	0.89	0.74
(phi)	5.31	0.9	3.88	6.82	1.00	1.00	1.00
mean_PPD	0.35	0.9	0.33	0.38	1.00	1.00	1.00
log-posterior	150.97	0.9	134.54	166.26	1.00	1.00	1.00

fun.sexit.gt(stmod_cov.wc$stanfit) %>%
  tab_header(title = "Core Winter Range", subtitle = "Willow Cover Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_cover/posterior_willowht_wc.rtf")

Probability of Direction

# Visualize the pd
p_direction(stmod_cov.wc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow cover, Core winter range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wc.pdf", width = 6.75, height =3.75)

Region of Practical Equivalence (ROPE)

rope.cov.wc.gt <- rope(stmod_cov.wc$stanfit, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 2,
    suffixing = TRUE
  ) 

rope.cov.wcwnc.gt %>%
  gt::gtsave(filename = "./output/tables/willow_cover/rope_willowcov_wc.rtf")

Contrasts

Results are given on the the response scale.

# library(emmeans)
wc.cov.emm.TCxF <- emmeans(stmod_cov.wc, ~ time_class | fenced)

# pairs(wc.cov.emm.TCxF, reverse = TRUE, type="response") %>% 
#   plot() +
#   theme_minimal() +
#   geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
#   labs(subtitle = "Point estimate displayed: median 
# HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Cover")

emm.func.tc.x.fenced(stmod_cov.wc) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Willow Cover")

ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wc.pdf", width = 6.75, height =3.75)

# table of contrasts
emmeans(stmod_cov.wc, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log odds ratio, not the response scale")

Pairwise contrasts
Results are given on the log odds ratio, not the response scale
contrast	fenced	estimate	lower.HPD	upper.HPD
2018 - BL	Fenced	2.73	2.05	3.41
2013 - BL	Fenced	0.68	0.04	1.35
2013 - 2018	Fenced	-2.04	-2.54	-1.51
2018 - BL	Unfenced	0.49	-0.08	1.02
2013 - BL	Unfenced	0.10	-0.49	0.66
2013 - 2018	Unfenced	-0.39	-0.84	0.09

emmeans(stmod_cov.wc, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale")

Pairwise contrasts
Results are given on the response scale
contrast	fenced	odds.ratio	lower.HPD	upper.HPD
2018 / BL	Fenced	15.26	6.65	27.70
2013 / BL	Fenced	1.98	0.88	3.51
2013 / 2018	Fenced	0.13	0.07	0.21
2018 / BL	Unfenced	1.63	0.85	2.66
2013 / BL	Unfenced	1.11	0.57	1.86
2013 / 2018	Unfenced	0.68	0.39	1.03

Willow Cover - Noncore Winter Range Macroplot Data

#### summary tables: time class and fenced
cov.wnc %>% 
  group_by(time_class, fenced) %>% 
  summarytools::descr(cover_allwillow,stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "WNC macroplot data: willow cover")

WNC macroplot data: willow cover
time_class	fenced	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Unfenced	cover_allwillow	0.37	0.28	0.03	0.30	1	32	100
2018	Unfenced	cover_allwillow	0.56	0.33	0.01	0.53	1	32	100
2013	Unfenced	cover_allwillow	0.43	0.32	0.02	0.37	1	27	100

Modeling and Posterior Description

#### Model
## add to min and subtract from max to conform
cov.wnc <- cov.wnc %>% 
  mutate(cover_allwillow = case_when(cover_allwillow == 1 ~ 0.999999,
                                     TRUE ~ cover_allwillow))

stmod_cov.wnc <- stan_glmer(cover_allwillow ~ time_class + (1 | site_id), 
                      data = cov.wnc,
                      family = mgcv::betar,
                      # prior_aux = exponential(2),
                      iter = 10000,
                      seed = 1234)

Prior summary: cover_allwillow ~ time_class + (1 | site_id)

prior_summary(stmod_cov.wnc)

## Priors for model 'stmod_cov.wnc' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0], scale = [2.5,2.5])
##   Adjusted prior:
##     ~ normal(location = [0,0], scale = [5.21,5.44])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

Effect Description

fun.sexit.gt(stmod_cov.wnc$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-0.22	0.9	-0.83	0.40	0.73	0.68	0.41
time_class2018	1.02	0.9	0.55	1.50	1.00	1.00	0.99
time_class2013	0.54	0.9	0.04	1.02	0.97	0.95	0.79
(phi)	3.04	0.9	2.12	3.99	1.00	1.00	1.00
mean_PPD	0.52	0.9	0.47	0.57	1.00	1.00	1.00
log-posterior	82.17	0.9	71.49	92.84	1.00	1.00	1.00

fun.sexit.gt(stmod_cov.wnc$stanfit) %>% 
  tab_header(title = "Noncore Winter Range", subtitle = "Willow Cover Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_cover/xposterior_willowht_wnc.rtf")

Probability of Direction

# Visualize the pd
p_direction(stmod_cov.wnc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow cover, Nnoncore winter range") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wnc.png", width = 6.75, height =3.75, dpi = 300)
# ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wnc.pdf", width = 6.75, height =3.75)

Region of Practical Equivalence (ROPE)

rope.cov.wnc.gt <- rope(stmod_cov.wnc$stanfit, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 2,
    suffixing = TRUE
  ) 

rope.cov.wnc.gt %>%
  gt::gtsave(filename = "./output/tables/willow_cover/rope_willowcov_wnc.rtf")

Contrasts

# library(emmeans)
wnc.cov.emm.TC <- emmeans(stmod_cov.wnc, ~ time_class)

pairs(wnc.cov.emm.TC, reverse = TRUE, type="response") %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Noncore Winter Range, Willow Cover. Results are given on the response scale.")

# ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wnc.png", width = 6.75, height =3.75, dpi = 300)

# table of contrasts
emmeans(stmod_cov.wnc, ~ time_class) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log odds ratio, not the response scale")

Pairwise contrasts
Results are given on the log odds ratio, not the response scale
contrast	estimate	lower.HPD	upper.HPD
2018 - BL	1.02	0.46	1.60
2013 - BL	0.54	-0.04	1.12
2013 - 2018	-0.49	-1.02	0.04

emmeans(stmod_cov.wnc, ~ time_class) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale")

Pairwise contrasts
Results are given on the response scale
contrast	odds.ratio	lower.HPD	upper.HPD
2018 / BL	2.78	1.40	4.61
2013 / BL	1.71	0.84	2.85
2013 / 2018	0.61	0.33	0.99

Willow Cover - Kawuneeche Valley Macroplot Data

#### summary tables: time class and fenced
cov.wk %>% 
  group_by(time_class, fenced) %>% 
  summarytools::descr(cover_allwillow,stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "WNC macroplot data: willow cover")

WNC macroplot data: willow cover
time_class	fenced	variable	mean	sd	min	med	max	n.valid	pct.valid
BL	Fenced	cover_allwillow	0.00	0.00	0.00	0.00	0.00	3	100
BL	Unfenced	cover_allwillow	0.00	0.00	0.00	0.00	0.00	5	100
2018	Fenced	cover_allwillow	0.46	0.27	0.14	0.48	0.75	4	100
2018	Unfenced	cover_allwillow	0.26	0.25	0.02	0.17	0.53	5	100
2017	Fenced	cover_allwillow	0.34	0.25	0.17	0.23	0.62	3	100
2017	Unfenced	cover_allwillow	0.23	0.24	0.02	0.16	0.58	4	100
2016	Fenced	cover_allwillow	0.31	0.20	0.19	0.20	0.54	3	100
2016	Unfenced	cover_allwillow	0.23	0.23	0.04	0.18	0.54	4	100

Modeling and Posterior Description

#### Model
## add to min and subtract from max to conform
cov.wk <- cov.wk %>% 
  mutate(cover_allwillow = cover_allwillow/100) %>% 
  mutate(cover_allwillow = case_when(cover_allwillow == 1 ~ 0.999999,
                                     cover_allwillow == 0 ~ 0.000001,
                                     TRUE ~ cover_allwillow)) 
  

stmod_cov.wk <- stan_glmer(cover_allwillow ~ time_class*fenced + (1 | site_id), 
                      data = cov.wk,
                      family = mgcv::betar,
                      # prior_aux = exponential(2),
                      iter = 8000,
                      seed = 1234)

Prior summary

prior_summary(stmod_cov.wk)

## Priors for model 'stmod_cov.wk' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.42,5.88,5.88,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

Effect Description

fun.sexit.gt(stmod_cov.wk$stanfit)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-4.17	0.9	-5.43	-3.02	1.00	1.00	1.00
time_class2018	1.14	0.9	-0.19	2.44	0.93	0.92	0.86
time_class2017	1.05	0.9	-0.28	2.54	0.89	0.88	0.82
time_class2016	1.04	0.9	-0.35	2.44	0.89	0.88	0.81
fencedUnfenced	0.00	0.9	-1.26	1.29	0.50	0.47	0.35
time_class2018:fencedUnfenced	-0.21	0.9	-1.90	1.45	0.58	0.56	0.46
time_class2017:fencedUnfenced	-0.14	0.9	-1.93	1.70	0.55	0.53	0.44
time_class2016:fencedUnfenced	-0.11	0.9	-1.91	1.67	0.54	0.52	0.43
(phi)	5.24	0.9	1.94	9.48	1.00	1.00	1.00
mean_PPD	0.04	0.9	0.01	0.07	1.00	0.26	0.00
log-posterior	132.75	0.9	126.83	137.90	1.00	1.00	1.00

fun.sexit.gt(stmod_cov.wk$stanfit) %>% 
  tab_header(title = "Kawuneeche Valley", subtitle = "Willow Cover Posterior Description") %>% 
  gt::gtsave(filename = "./output/tables/willow_cover/xposterior_willowht_wk.rtf")

Probability of Direction

# Visualize the pd
p_direction(stmod_cov.wk) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow cover, KV") +
  scale_fill_manual(values = colfunc3(2))

ggsave("./output/figures_202108/willow_cover_stats/pd_willowht_wk.png", width = 6.75, height =3.75, dpi = 300)

Region of Practical Equivalence (ROPE)

rope.cov.wk.gt <- rope(stmod_cov.wk$stanfit, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 2,
    suffixing = TRUE
  ) 

rope.cov.wk.gt %>%
  gt::gtsave(filename = "./output/tables/willow_cover/rope_willowcov_wk.rtf")

Contrasts

# library(emmeans)

emm.func.tc.x.fenced(stmod_cov.wk) +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Kawuneeche Valley, Willow Cover. Results are given on the response scale.")

ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wk.png", width = 6.75, height =3.75, dpi = 300)

emm.func.tc.x.fenced(stmod_cov.wk) +
  labs(caption = "Point estimate displayed: median 
HPD interval probability: 0.95 \n Willow Cover. Results are given on the response scale.", title = "Kawuneeche Valley", x = "Estimate", y = "Contrast")

ggsave("./output/figures_202108/willow_cover_stats/emm_willowht_TCxF_wk_alt.png", width = 5.75, height =5.75, dpi = 300)

Wild Basin Comparisons

# read in data and preprocess for joining with EVMP data
wbasin <- read_csv("./data/wild_basin/Wild Basin willow data survey 2018.csv")

wbasin <- wbasin %>%
  janitor::clean_names()
  
wbasin.sel <- wbasin %>% 
  rename(plant_ht_cm = max_ht_cm) %>% 
  mutate(yr = 2018) %>% 
  mutate(fenced = "N") %>% 
  mutate(site_type = "WB") %>% 
  mutate(site_id = paste0("WB",point)) %>% 
  dplyr::select(c(yr, site_type, site_id, plant_ht_cm, species, fenced, location)) %>% 
  distinct()

## process the evmp into a bindable format with wb data
willow.mcro.evmp <- willow.ht %>% 
  rename(species = species_code) %>% 
  dplyr::select(yr, site_type, species, site_id, plant_ht_cm, fenced, location) %>% 
  distinct() %>% 
  filter(yr == 2018)

#### bind wb with 2018 evmp
wb.evmp18 <- bind_rows(wbasin.sel, willow.mcro.evmp) %>% 
  distinct()

### Willow height comparison among Wild Basin and EVMP plots
## compare WILLOW height

wb.evmp18 %>%
  filter(str_detect(species, "^S")) %>% 
  ggplot() +
  geom_density(aes(fill = site_type, color = site_type, lty = site_type, x = plant_ht_cm), size=1, alpha = .25) +
  theme_minimal() +
  scale_fill_manual(values = colfunc1(4)) +
  scale_color_grey() +
  # scale_fill_manual(values = colfunc3_RdBu(4)) +
  # scale_color_manual(values = colfunc3_RdBu(4)) +
  labs(x = "Willow height (cm)", y  = "Density", lty= "Site type", color = "Site type", fill = "Site type", caption = "willow height, WB + EVMP, Salix only")

ggsave("./output/figures_202107/WB_EVMP_willow_height.png", dpi = 300, width = 4.75, height = 3.75)

wb.evmp18 <- wb.evmp18 %>%
  mutate(site_type.lbl = case_when(site_type == "WB" ~ "Wild Basin",
                               site_type == "WC" ~ "Core winter range",
                               site_type == "WNC" ~ "Noncore winter range",
                               site_type == "WK" ~ "Kawuneeche Valley",
                               TRUE ~ "other")) %>% 
  mutate(fenced = case_when(fenced == "N" ~ "Unfenced",
                            fenced == "Y" ~ "Fenced",
                            TRUE ~ fenced)) %>% 
  mutate(site.type.fenced = paste0(site_type.lbl,"\n",fenced))

# wb.evmp18 %>% distinct(site.type.fenced)
wb.evmp18.willows <- wb.evmp18 %>% 
  filter(str_detect(species, "^S"))

wb.evmp18 %>%
  filter(str_detect(species, "^S")) %>% 
  ggplot(aes(x=reorder(site.type.fenced, plant_ht_cm), plant_ht_cm)) +
  geom_boxplot(aes(fill=site_type)) +
  # geom_density(aes(fill = site_type, color = site_type, lty = site_type, x = plant_ht_cm), alpha = 0.05) +
  theme_minimal() +
  scale_fill_manual(values = colfunc2(4)) +
  scale_color_manual(values = colfunc2(4)) +
  theme(legend.position = "none") +
  # theme(axis.text.x = element_text(angle = 0, hjust = 1)) +
  labs(x = "", y  = "Height (cm)", lty= "Site type", color = "Site type", fill = "Site type", caption = "willow height, WB + EVMP, Salix only")

ggsave("./output/figures_202107/WB_EVMP_willow_height_boxplot.png", dpi = 300, width = 6.75, height = 4.75)

wb.evmp18 %>%
  group_by(site.type.fenced) %>%
  descr(plant_ht_cm, stat="common") %>% 
  tb() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  arrange(-mean) %>%
  datatable(caption = "All willow species combined")

wb.evmp18 %>%
  group_by(site.type.fenced, species) %>%
  descr(plant_ht_cm) %>% 
  tb() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  arrange(-mean) %>%
  datatable(caption = "By willow species")

wb.evmp18.willows.loc <- wb.evmp18.willows %>% 
  select(-yr) %>% 
  group_by(site_id,location) %>% 
  summarytools::descr(stats = "common") %>% 
  tb() %>%
  select(c(location,mean)) %>% 
  group_by(location) %>% 
  descr(mean,stats = "common") %>% 
  tb() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  select(-pct.valid)

wb.evmp18.willows.loc %>%
  arrange(-mean) %>% 
  gt()

location	variable	mean	sd	min	med	max	n.valid
Fishermans Nook	mean	290.0	7.1	285.0	290.0	295.0	2
Hollowell Park	mean	223.6	124.0	83.5	194.9	508.8	12
Wild Basin	mean	216.6	112.5	10.0	210.0	413.3	27
Lower Hidden Valley	mean	216.4	NA	216.4	216.4	216.4	1
McGraw Ranch	mean	183.1	76.6	94.3	170.8	309.2	8
Hidden Valley	mean	160.5	104.8	67.2	116.7	377.5	9
Upper Beaver Meadows	mean	154.7	90.2	55.0	146.6	280.0	8
Horseshoe Park	mean	143.2	63.2	49.0	168.6	226.8	12
Moraine Park	mean	120.1	62.3	30.0	122.0	270.0	33
Endovalley	mean	112.0	69.8	57.0	90.7	287.2	9
Kawuneeche	mean	67.0	32.0	25.0	73.5	111.7	9

# %>%
#   gt::gtsave(filename = "./output/tables/WB_vs_EVMP_willowHt_loc.rtf")

#### range type
wb.evmp18.willows.rt <- wb.evmp18.willows %>% 
  select(-yr) %>% 
  group_by(site_id,site_type) %>% 
  summarytools::descr(stats = "common") %>% 
  tb() %>%
  select(c(site_type,mean)) %>% 
  group_by(site_type) %>% 
  descr(mean,stats = "common") %>% 
  tb() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  select(-pct.valid)

wb.evmp18.willows.rt %>%
  arrange(-mean) %>%
  select(-variable) %>%
  gt() #%>%

site_type	mean	sd	min	med	max	n.valid
WB	216.6	112.5	10.0	210.0	413.3	27
WNC	199.7	104.4	67.2	179.8	508.8	32
WC	127.9	67.4	30.0	118.7	287.2	62
WK	67.0	32.0	25.0	73.5	111.7	9

  # gt::gtsave(filename = "./output/tables/WB_vs_EVMP_willowHt_site_type.rtf")

wb.evmp18 %>%
  filter(str_detect(species, "^S")) %>% 
  ggplot(aes(x=reorder(site.type.fenced, plant_ht_cm), plant_ht_cm)) +
  geom_boxplot(aes(fill=site_type)) +
  # geom_density(aes(fill = site_type, color = site_type, lty = site_type, x = plant_ht_cm), alpha = 0.05) +
  theme_minimal() +
  scale_fill_manual(values = colfunc2(4)) +
  scale_color_manual(values = colfunc2(4)) +
  theme(legend.position = "none") +
  # theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "", y  = "Height (cm)", lty= "Site type", color = "Site type", fill = "Site type", caption = "willow height, WB + EVMP, Salix only")

ggsave("./output/figures_202107/WB_EVMP_willow_height_boxplot.png", dpi = 300, width = 5.75, height = 4.75)

# 2 panel
pl.boxht.wc.f <- wb.evmp18 %>%
  filter(str_detect(species, "^S")) %>% 
  filter(site_type == "WC" & fenced == "Fenced") %>% 
  ggplot(aes(x=reorder(site.type.fenced, plant_ht_cm), plant_ht_cm)) +
  geom_boxplot(aes(fill=site_type)) +
  # geom_density(aes(fill = site_type, color = site_type, lty = site_type, x = plant_ht_cm), alpha = 0.05) +
  theme_minimal() +
  scale_fill_manual(values = colfunc2(5)) +
  scale_color_manual(values = colfunc2(5)) +
  theme(legend.position = "none") +
  # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y  = "Height (cm)", lty= "Site type", color = "Site type", fill = "Site type", caption = "willow height,  all willow spp") +
  ylim(0,600) +
  facet_wrap(~fenced, scales = "free_x")

pl.boxht.uf.all <- wb.evmp18 %>%
  filter(str_detect(species, "^S")) %>% 
  filter(fenced == "Unfenced") %>% 
  ggplot(aes(x=reorder(site.type.fenced, plant_ht_cm), plant_ht_cm)) +
  geom_boxplot(aes(fill=site_type)) +
  # geom_density(aes(fill = site_type, color = site_type, lty = site_type, x = plant_ht_cm), alpha = 0.05) +
  theme_minimal() +
  scale_fill_manual(values = colfunc2(4)) +
  scale_color_manual(values = colfunc2(4)) +
  theme(legend.position = "none") +
  # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y  = "", lty= "Site type", color = "Site type", fill = "Site type") +
  ylim(0,600) +
  facet_wrap(~fenced, scales = "free_x")


pl.boxht.wc.f + pl.boxht.uf.all + plot_layout(widths = c(.55,2))

ggsave("./output/figures_202107/WB_EVMP_willow_height_boxplot2.png", dpi = 300, width = 6.75, height = 4.75)

Modeling and Posterior Description

wb_rt2stan <- wb.evmp18.willows %>% 
  select(-yr) %>% 
  group_by(site_id,site_type) %>% 
  summarytools::descr(stats = "common") %>% 
  tb() %>%
  select(c(site_type,mean))

wb_stanaov <- stan_aov(mean ~ site_type, data = wb_rt2stan, 
                  prior = R2(location = 0.25), adapt_delta = 0.999,
                  iter = 8000,
                  seed = 12345)
wb_stanaov

## stan_aov
##  family:       gaussian [identity]
##  formula:      mean ~ site_type
##  observations: 130
##  predictors:   4
## ------
##              Median MAD_SD
## (Intercept)   214.1   16.4
## site_typeWC   -84.9   19.7
## site_typeWK  -143.0   33.2
## site_typeWNC  -16.2   21.9
## 
## Auxiliary parameter(s):
##               Median MAD_SD
## R2             0.2    0.1  
## log-fit_ratio  0.0    0.1  
## sigma         86.7    5.4  
## 
## ANOVA-like table:
##                   Median  MAD_SD 
## Mean Sq site_type 89042.0 29010.3
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

# report::report(wb_stanaov)

wb_stanglmer <- stan_lmer(mean ~ 1 + (1|site_type),
                   data = wb_rt2stan, prior_intercept = cauchy(),
                   iter=8000,
                   prior_covariance = decov(shape = 2, scale = 2),
                   adapt_delta = 0.999, seed = 12345)

summary(wb_stanglmer)#

## 
## Model Info:
##  function:     stan_lmer
##  family:       gaussian [identity]
##  formula:      mean ~ 1 + (1 | site_type)
##  algorithm:    sampling
##  sample:       16000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 130
##  groups:       site_type (4)
## 
## Estimates:
##                                            mean     sd       10%      50%   
## (Intercept)                                  11.8     39.1     -5.1      0.6
## b[(Intercept) site_type:WB]                 202.9     42.6    169.4    211.7
## b[(Intercept) site_type:WC]                 115.7     40.4     89.4    124.9
## b[(Intercept) site_type:WK]                  54.3     46.8      6.7     60.4
## b[(Intercept) site_type:WNC]                186.4     41.9    156.3    195.2
## sigma                                        87.3      5.5     80.4     87.0
## Sigma[site_type:(Intercept),(Intercept)]  56309.2  61804.0  13914.1  37753.1
##                                            90%   
## (Intercept)                                  35.3
## b[(Intercept) site_type:WB]                 236.3
## b[(Intercept) site_type:WC]                 142.7
## b[(Intercept) site_type:WK]                 101.6
## b[(Intercept) site_type:WNC]                218.2
## sigma                                        94.5
## Sigma[site_type:(Intercept),(Intercept)] 117135.8
## 
## Fit Diagnostics:
##            mean   sd    10%   50%   90%
## mean_PPD 158.7   10.9 144.8 158.5 172.6
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                                          mcse   Rhat   n_eff
## (Intercept)                                 1.4    1.0   760
## b[(Intercept) site_type:WB]                 1.4    1.0   873
## b[(Intercept) site_type:WC]                 1.4    1.0   825
## b[(Intercept) site_type:WK]                 1.4    1.0  1177
## b[(Intercept) site_type:WNC]                1.4    1.0   859
## sigma                                       0.1    1.0  7279
## Sigma[site_type:(Intercept),(Intercept)] 1158.7    1.0  2845
## mean_PPD                                    0.1    1.0 16701
## log-posterior                               0.1    1.0  1048
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

describe_posterior(
  wb_stanglmer,
  effects = "all",
  component = "all",
  test = c("p_direction", "p_significance"),
  centrality = "all"
) %>% 
  gt()

Parameter	Effects	Component	Median	Mean	MAP	CI	CI_low	CI_high	pd	ps	Rhat	ESS
(Intercept)	fixed	conditional	6.474931e-01	11.76255	-6.225121e-03	0.95	-28.35513	118.91765	0.5757500	0.1733125	1.004323	759.9500
b[(Intercept) site_type:WB]	random	conditional	2.116868e+02	202.92684	2.132987e+02	0.95	94.91439	260.35942	0.9959375	0.9943125	1.003650	872.5293
b[(Intercept) site_type:WC]	random	conditional	1.249110e+02	115.65689	1.258237e+02	0.95	10.18839	165.68727	0.9643125	0.9571875	1.004226	824.7230
b[(Intercept) site_type:WK]	random	conditional	6.043394e+01	54.28535	6.462747e+01	0.95	-48.29569	137.09116	0.9135000	0.8911250	1.002778	1176.6734
b[(Intercept) site_type:WNC]	random	conditional	1.952310e+02	186.42151	1.994495e+02	0.95	76.97108	240.74741	0.9933125	0.9912500	1.004171	859.4127
Sigma[site_type:(Intercept),(Intercept)]	random	conditional	3.775306e+04	56309.22616	2.521362e+04	0.95	858.97442	161849.86698	1.0000000	1.0000000	1.000574	2845.2345
sigma	fixed	sigma	8.698569e+01	87.26022	8.631719e+01	0.95	76.70337	97.97795	1.0000000	1.0000000	NA	NA

Session info

R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)
Run date: 2021 August 28

Package Information
Package	Version	Reference
bayesplot	1.8.1	Gabry J, Mahr T (2021). "bayesplot: Plotting for Bayesian Models." Rpackage version 1.8.1, <URL: https://mc-stan.org/bayesplot/>.
bayestestR	0.10.5	Makowski, D., Ben-Shachar, M., & Lüdecke, D. (2019). bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework. Journal of Open Source Software, 4(40), 1541. doi:10.21105/joss.01541
dplyr	1.0.7	Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2021). dplyr: A Grammar of Data Manipulation. R package version 1.0.7. https://CRAN.R-project.org/package=dplyr
DT	0.18	Yihui Xie, Joe Cheng and Xianying Tan (2021). DT: A Wrapper of the JavaScript Library 'DataTables'. R package version 0.18. https://CRAN.R-project.org/package=DT
emmeans	1.6.2.1	Russell V. Lenth (2021). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.6.2-1. https://CRAN.R-project.org/package=emmeans
forcats	0.5.1	Hadley Wickham (2021). forcats: Tools for Working with Categorical Variables (Factors). R package version 0.5.1. https://CRAN.R-project.org/package=forcats
ggplot2	3.3.5	H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
ggridges	0.5.3	Claus O. Wilke (2021). ggridges: Ridgeline Plots in 'ggplot2'. R package version 0.5.3. https://CRAN.R-project.org/package=ggridges
gt	0.3.1	Richard Iannone, Joe Cheng and Barret Schloerke (2021). gt: Easily Create Presentation-Ready Display Tables. R package version 0.3.1. https://CRAN.R-project.org/package=gt
janitor	2.1.0	Sam Firke (2021). janitor: Simple Tools for Examining and Cleaning Dirty Data. R package version 2.1.0. https://CRAN.R-project.org/package=janitor
knitr	1.33	Yihui Xie (2021). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.33.
lubridate	1.7.10	Garrett Grolemund, Hadley Wickham (2011). Dates and Times Made Easy with lubridate. Journal of Statistical Software, 40(3), 1-25. URL https://www.jstatsoft.org/v40/i03/.
patchwork	1.1.1	Thomas Lin Pedersen (2020). patchwork: The Composer of Plots. R package version 1.1.1. https://CRAN.R-project.org/package=patchwork
purrr	0.3.4	Lionel Henry and Hadley Wickham (2020). purrr: Functional Programming Tools. R package version 0.3.4. https://CRAN.R-project.org/package=purrr
R	4.0.3	R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
Rcpp	1.0.7	Dirk Eddelbuettel and Romain Francois (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. URL https://www.jstatsoft.org/v40/i08/.
readr	2.0.1	Hadley Wickham and Jim Hester (2021). readr: Read Rectangular Text Data. R package version 2.0.1. https://CRAN.R-project.org/package=readr
rstanarm	2.21.1	Goodrich B, Gabry J, Ali I & Brilleman S. (2020). rstanarm: Bayesian applied regression modeling via Stan. R package version 2.21.1 https://mc-stan.org/rstanarm.
stringr	1.4.0	Hadley Wickham (2019). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.4.0. https://CRAN.R-project.org/package=stringr
summarytools	1.0.0	Dominic Comtois (2021). summarytools: Tools to Quickly and Neatly Summarize Data. R package version 1.0.0. https://CRAN.R-project.org/package=summarytools
tibble	3.1.3	Kirill Müller and Hadley Wickham (2021). tibble: Simple Data Frames. R package version 3.1.3. https://CRAN.R-project.org/package=tibble
tidyr	1.1.3	Hadley Wickham (2021). tidyr: Tidy Messy Data. R package version 1.1.3. https://CRAN.R-project.org/package=tidyr
tidyverse	1.3.1	Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
viridis	0.6.1	Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco Sciaini, and Cédric Scherer (2021). Rvision - Colorblind-Friendly Color Maps for R. R package version 0.6.1.
viridisLite	0.4.0	Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco Sciaini, and Cédric Scherer (2021). Rvision - Colorblind-Friendly Color Maps for R. R package version 0.4.0.

References

Brilleman, SL, MJ Crowther, M Moreno-Betancur, J Buros Novik, and R Wolfe. 2018. “StatCon 2018. 10-12 Jan 2018.” In. Pacific Grove, CA, USA. https://github.com/stan-dev/stancon_talks/.

Gabry, Jonah, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. 2019. “Visualization in Bayesian Workflow.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 182 (2): 389–402. https://doi.org/10.1111/rssa.12378.

Goodrich, B, J Gabry, I Ali, and S Brilleman. 2020. “Rstanarm: Bayesian Applied Regression Modeling via Stan. R Package Version 2.19.3 Https://Mc-Stan.org/Rstanarm.”

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, Alex Hayes, Lionel Henry, and Jim Hester. 2019. “Welcome to the Tidyverse.” Journal of Open Source Software 4 (43): 1686.

Zeigenfuss, Linda C., and Therese L. Johnson. 2015. “Monitoring of Vegetation Response to Elk Population and Habitat Management in Rocky Mountain National Park, 200814.” Reston, VA. https://doi.org/10.3133/ofr20151216.

Zeigenfuss, Linda C, Therese L Johnson, and Zachary Wiebe. 2011. “Monitoring Plan for Vegetation Responses to Elk Management in Rocky Mountain National Park.”

Appendix D. EVMP Willow Statistical Analyses

Introduction

Methods

Results

Willow Height - Combined Winter Range Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Effects of Burning

Effect Description

Willow Height - Core Winter Range Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Willow Height - Noncore Winter Range Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Willow Height - Kawuneeche Valley Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Willow Cover - Combined Winter Range Macroplot Data

Summary Tables

Histograms of Willow Cover

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Effects of Burning - Combined Range

Effect Description

Willow Cover - Core Winter Range Macroplot Data

Summary Tables

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Willow Cover - Noncore Winter Range Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Willow Cover - Kawuneeche Valley Macroplot Data

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Wild Basin Comparisons

Modeling and Posterior Description

Session info

References