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

# 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 Elk Vegetation Management Plan (EVMP) data collected through the 2018 sampling season. Analyses are limited to data from plots established in aspen communities in areas of the elk winter range (core and noncore areas) and the Kawuneeche Valley (i.e., “aspen” plots, see Zeigenfuss et al. 2011) and is 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 this 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 start from and build upon derived data produced in a separate code document focused on ingesting, compiling, and cleaning raw data provided by RMNP staff. 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., willow 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 repeated measures analyses were fit separately for combined core and noncore winter range, core winter range, and noncore winter range plots with plots treated as random factors using the ‘stan_glmer’ function in the rstanrm package (Goodrich et al. 2020)(Brilleman et al. 2018). In addition, models were fit to particular subgroups of data (e.g., short saplings vs tall saplings).

Bayesian estimation was performed via MCMC adding independent weakly informed priors specific to the data type being modeled. Count data like aspen stem counts were modeled as poisson processes, while proportion data (e.g., cover) were modeled using a beta distribution (Ferrari and Cribari-Neto 2004), which constrains values from 0 to 1.

Bayesian models view a model parameter \(\theta\) as a [random variable]. In contrast, 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 in hand 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 analytic solution, Markov chain Monte Carlo (MCMC) methods are used to numerically solve the posterior density by directly generating random draws of parameters.
General steps in a Bayesian analysis include:

Specify a joint distribution for the outcome(s) and unknowns, which is 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 analyses were fit separately for combined core and non-core winter range, core winter range (WC) and non-core winter range plots (WNC) with plots treated as random factors using the ‘stan_glmer’ function in the rstanrm package (Goodrich et al. 2020)(Brilleman et al. 2018). Bayesian estimation was performed via MCMC adding independent weakly informed priors specific to the data type being modeled.

Likelihood

If the outcome for an observation \(y\) is assumed to follow a Poisson distribution (e.g., aspen stem count data), the likelihood for one observation can be written as:

\[\tfrac{1}{y!} \lambda^y e^{-\lambda},\]

where \(\lambda = E(y | \mathbf{x}) = g^{-1}(\eta)\) and \(\eta = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) is a linear predictor. For the Poisson distribution it is also true that \(\lambda = Var(y | \mathbf{x})\), i.e. the mean and variance are both \(\lambda\). The rate parameter \(\lambda\) must be positive, so with a Poisson GLM, the link function \(g\) maps between the positive real numbers \(\mathbb{R}^+\) (thesupport of \(\lambda\)) and the set of all real numbers \(\mathbb{R}\). When applied to a linear predictor \(\eta\) with values in \(\mathbb{R}\), the inverse link function \(g^{-1}(\eta)\) returns a positive real number.

The standard link function for a Poisson GLM is the log link \(g(x) = \ln{(x)}\). With the log link, the inverse link function is the exponential function and the likelihood for a single observation is:

\[\frac{g^{-1}(\eta)^y}{y!} e^{-g^{-1}(\eta)} = \frac{e^{\eta y}}{y!} e^{-e^\eta}.\]

Posterior

With independent prior distributions, the joint posterior distribution for \(\alpha\) and \(\boldsymbol{\beta}\) in the Poisson model is proportional to the product of the priors and the \(N\) likelihood contributions:

\[f\left(\alpha,\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 { \frac{g^{-1}(\eta_i)^{y_i}}{y_i!} e^{-g^{-1}(\eta_i)}}.\]

This is the posterior distribution drawn from when using MCMC.
The Probability of Direction (PD) was used to represent certainty associated with the most probable direction (positive or negative) of the effect (e.g., time class, fencing). The PD is correlated with the frequentist p-value, with a two-sided p-value of respectively .1, .05, and .01 approximated by a PD of 95%, 97.5%, and 99.5% (Makowski, Ben-Shachar, and Lüdecke 2019). The “region of practical equivalence” (ROPE) was used to evaluate the probability of a parameter being outside a range that can be considered as “practically no effect,” i.e., a region enclosing values that are equivalent to the null value for practical purposes (Kruschke 2018). The proportion of the 95% HDI credible interval falling within the ROPE was used as an index for an analog to frequentist “null-hypothesis” testing.

Aspen count

## read in data
# !!! check the read in below - a major change

# Aspen dbh
# ht_perc1 <- read_csv("./output/exported_data/asp_ht_perc1_20200309.csv")
dbh_perc1 <- read_csv("./output/exported_data/asp_dbh_perc1_20210601.csv")

ht_perc1 <- read_csv("./output/exported_data/asp_ht_perc1_20200309.csv")


# Data munging

## filter out in-between years
ht_perc1 <- ht_perc1 %>%
  filter(timeClass %in% c("BL","2013", "2018")) 

### factor relevel
ht_perc1 <- ht_perc1 %>% 
  mutate(timeClass = as.factor(timeClass)) %>% 
  mutate(timeClass = fct_relevel(timeClass, "BL", "2013", "2018")) %>% 
  mutate(RANGE_TYPE = as.factor(RANGE_TYPE)) %>% 
  mutate(RANGE_TYPE = fct_relevel(RANGE_TYPE, "core winter range","non-core winter range","Kawuneeche Valley"))   

dbh_perc1 <- dbh_perc1 %>% 
  filter(timeClass %in% c("BL","2013", "2018")) %>% 
  mutate(timeClass = as.factor(timeClass)) %>% 
  mutate(timeClass = fct_relevel(timeClass, "BL", "2013", "2018")) %>% 
  mutate(RANGE_TYPE = as.factor(RANGE_TYPE)) %>% 
  mutate(RANGE_TYPE = fct_relevel(RANGE_TYPE, "core winter range","non-core winter range","Kawuneeche Valley"))

# factor conversions...

ht_perc1 <- ht_perc1 %>% 
  clean_names() %>% 
  rename(htclass2 = h_tclass2) %>% 
  rename(htclass = h_tclass)

## set factor levels for burned and fenced
ht_perc1 <- ht_perc1 %>%
  mutate(fenced.long = case_when(fenced == "N" ~ "Unfenced",
                                 fenced == "Y" ~ "Fenced",
                                 TRUE ~ fenced)) %>% 
  mutate(burned = case_when(burned == "Not burned" ~ "Unburned",
                            TRUE ~ burned)) %>%
  mutate(fenced.long = as_factor(fenced.long)) %>%
  mutate(burned = as_factor(burned))

# convert 'fenced' to factor
ht_perc1 <- ht_perc1 %>%
  mutate(fenced = as_factor(fenced)) %>% 
  mutate(burned = as_factor(burned))

# levels(ht_perc1$fenced)
# levels(ht_perc1$burned)

## subset height data to all LIVE stems
asp.acanc <- ht_perc1 %>% 
  filter(site_type == "AC" | site_type == "ANC") %>% 
  filter(live_dead == "LIVE")
  
asp.ac <- ht_perc1 %>% 
  filter(site_type == "AC") %>% 
  filter(live_dead == "LIVE")

# non core
asp.anc <- ht_perc1 %>% 
  filter(site_type == "ANC") %>% 
  filter(live_dead == "LIVE")

# KV
asp.ak <- ht_perc1 %>% 
  filter(site_type == "AK") %>% 
  filter(live_dead == "LIVE")

### subset to small saplings <1.5m
## subset to all LIVE stems
ht_perc1 <- ht_perc1 %>% 
  mutate(sap_size = case_when(htclass2 == "0-50cm" | htclass2 == "51-100cm" | htclass2 == "101-150cm" ~ "small saplings",
                              htclass2 == "151-200cm" | htclass2 == "201-250cm" ~ "tall saplings")) 

asp.sm.acanc <- ht_perc1 %>%
  filter(sap_size == "small saplings") %>% 
  filter(site_type == "AC" | site_type == "ANC") %>% 
  filter(live_dead == "LIVE")
  
asp.sm.ac <- ht_perc1 %>% 
  filter(sap_size == "small saplings") %>%
  filter(site_type == "AC") %>% 
  filter(live_dead == "LIVE")

# non core
asp.sm.anc <- ht_perc1 %>% 
  filter(sap_size == "small saplings") %>%
  filter(site_type == "ANC") %>% 
  filter(live_dead == "LIVE")

# KV
asp.sm.ak <- ht_perc1 %>% 
  filter(sap_size == "small saplings") %>%
  filter(site_type == "AK") %>% 
  filter(live_dead == "LIVE")

### subset to large saplings 1.5m-2.5
## subset to all LIVE stems
asp.lg.acanc <- ht_perc1 %>% 
  filter(sap_size == "tall saplings") %>% 
  filter(site_type == "AC" | site_type == "ANC") %>% 
  filter(live_dead == "LIVE")
  
asp.lg.ac <- ht_perc1 %>% 
  filter(sap_size == "tall saplings") %>%
  filter(site_type == "AC") %>% 
  filter(live_dead == "LIVE")

# non core
asp.lg.anc <- ht_perc1 %>% 
  filter(sap_size == "tall saplings") %>%
  filter(site_type == "ANC") %>% 
  filter(live_dead == "LIVE")

# KV
asp.lg.ak <- ht_perc1 %>% 
  filter(sap_size == "tall saplings") %>%
  filter(site_type == "AK") %>% 
  filter(live_dead == "LIVE")

### dbh
dbh_perc1 <- dbh_perc1 %>% 
  clean_names()

dbh_perc1 <- dbh_perc1 %>%
  mutate(fenced.long = case_when(fenced == "N" ~ "Unfenced",
                                 fenced == "Y" ~ "Fenced",
                                 TRUE ~ fenced)) %>% 
  mutate(burned = case_when(burned == "Not burned" ~ "Unburned",
                            TRUE ~ burned)) %>%
  mutate(fenced.long = as_factor(fenced.long)) %>%
  mutate(burned = as_factor(burned))

Results

Comparison of Sapling Classes

ht_perc1 %>%
  filter(live_dead == "LIVE") %>%
  group_by(sap_size, site_id, time_class, fenced, burned) %>% 
  summarise(sum.st = sum(stem_tally)) %>% 
  ungroup() %>% 
  ggplot() +
  geom_col(aes(x=site_id, y=sum.st, fill=sap_size)) +
  facet_grid(time_class ~ fenced)

tally.sapsize <- ht_perc1 %>%
  filter(live_dead == "LIVE") %>%
  group_by(sap_size, site_id, site_type, time_class, fenced, burned) %>% 
  summarise(sum.st = sum(stem_tally)) %>% 
  ungroup() 


sapsize_summaryTCxSiteType <- tally.sapsize %>%
  group_by(time_class, site_type, sap_size) %>%  
  summarytools::descr(var = sum.st, round.digits = 1) %>% 
  summarytools::tb() %>%
  mutate(across(where(is.numeric), round, 1)) %>% 
  select(-c(variable, se.skewness, kurtosis, pct.valid,cv, skewness))
  
sapsize_summaryTCxSiteType %>% 
  rename(n = n.valid, 'Time class' = time_class) %>% 
  arrange(sap_size,site_type,'time_class') %>% 
  gt() %>% 
  tab_header(title = "Aspen Saplings")#

Aspen Saplings
Time class	site_type	sap_size	mean	sd	min	q1	med	q3	max	mad	iqr	n
BL	AC	small saplings	25.4	30.3	0	7.0	17.5	30.0	139	16.3	21.5	40
2013	AC	small saplings	45.0	82.8	0	11.0	22.0	45.0	523	22.2	34.0	45
2018	AC	small saplings	25.0	29.3	0	5.0	11.5	34.0	111	12.6	27.2	46
BL	AK	small saplings	10.1	6.6	3	4.5	9.0	14.5	22	6.7	9.0	8
2013	AK	small saplings	12.2	9.8	3	4.0	11.5	16.0	32	8.9	12.0	8
2018	AK	small saplings	14.8	11.3	0	6.0	15.0	20.5	35	8.9	11.8	8
BL	ANC	small saplings	11.9	12.4	0	3.5	7.5	16.0	45	6.7	11.8	20
2013	ANC	small saplings	14.4	12.3	1	4.0	14.0	20.0	53	14.8	16.0	21
2018	ANC	small saplings	21.0	14.8	2	6.0	18.0	33.0	51	19.3	27.0	21
BL	AC	tall saplings	0.1	0.3	0	0.0	0.0	0.0	1	0.0	0.0	40
2013	AC	tall saplings	3.6	10.1	0	0.0	0.0	1.0	49	0.0	1.0	45
2018	AC	tall saplings	5.3	12.6	0	0.0	0.0	3.0	73	0.0	3.0	46
BL	AK	tall saplings	0.0	0.0	0	0.0	0.0	0.0	0	0.0	0.0	8
2013	AK	tall saplings	0.0	0.0	0	0.0	0.0	0.0	0	0.0	0.0	8
2018	AK	tall saplings	0.0	0.0	0	0.0	0.0	0.0	0	0.0	0.0	8
BL	ANC	tall saplings	0.8	2.9	0	0.0	0.0	0.0	13	0.0	0.0	20
2013	ANC	tall saplings	1.0	2.3	0	0.0	0.0	1.0	10	0.0	1.0	21
2018	ANC	tall saplings	1.9	2.8	0	0.0	0.0	3.0	10	0.0	3.0	21

## barchart
sapsize_summaryTCxSiteType %>% 
  mutate(site_type = as_factor(site_type)) %>% 
  mutate(site_type = fct_relevel(site_type, "AK", after = 2)) %>% 
  ggplot(aes(x=time_class, y=mean)) +
  geom_col(aes(fill = sap_size), position = 'dodge') +
  scale_fill_manual(values=colfunc3(2)) +
  facet_wrap(~site_type) +
  theme_minimal() +
  labs(x = "Time class", y = "Mean stem tally", fill = "")

ggsave("./output/figures_exported/aspen_stem_count/asp_saplingmean_3panel.png", dpi = 300, width = 8.75, height = 4.75)

## separate panels
pl.sapmean.ac <- tally.sapsize %>%
  filter(site_type == "AC") %>% 
  group_by(time_class, sap_size, fenced) %>%  
  summarytools::descr(var = sum.st, round.digits = 1) %>% 
  summarytools::tb() %>%
  mutate(across(where(is.numeric), round, 1)) %>% 
  ggplot(aes(x=time_class, y=mean)) +
  geom_col(aes(fill = fenced), position = 'dodge') +
  scale_fill_manual(values=colfunc3(2)) +
  facet_wrap(~sap_size) +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(title = "Core Winter Range", x = "Time class", y = "Mean stem tally", fill = "")

## anc
pl.sapmean.anc <- tally.sapsize %>%
  filter(site_type == "ANC") %>% 
  group_by(time_class, sap_size, fenced) %>%  
  summarytools::descr(var = sum.st, round.digits = 1) %>% 
  summarytools::tb() %>%
  mutate(across(where(is.numeric), round, 1)) %>% 
  ggplot(aes(x=time_class, y=mean)) +
  geom_col(aes(fill = fenced), position = 'dodge') +
  scale_fill_manual(values=colfunc3(2)) +
  facet_wrap(~sap_size) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Noncore Winter Range", x = "Time class", y = "Mean stem tally", fill = "")

## ak
pl.sapmean.ak <- tally.sapsize %>%
  filter(site_type == "AK") %>% 
  group_by(time_class, sap_size, fenced) %>%  
  summarytools::descr(var = sum.st, round.digits = 1) %>% 
  summarytools::tb() %>%
  mutate(across(where(is.numeric), round, 1)) %>% 
  ggplot(aes(x=time_class, y=mean)) +
  geom_col(aes(fill = fenced), position = 'dodge') +
  scale_fill_manual(values=colfunc3(2)) +
  facet_wrap(~sap_size) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Kawuneeche Valley", x = "Time class", y = "Mean stem tally", fill = "")

pl.sapmean.ac + pl.sapmean.anc + pl.sapmean.ak + patchwork::plot_layout(ncol=3)

# save plots
ggsave("./output/figures_exported/aspen_stem_count/asp_saplingmean_3panelalt.png", dpi = 300, width = 5.75, height = 6.75)

asp.sm.acanc %>%
  filter(LIVE_DEAD == 'LIVE') %>% 
  filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>% 
  mutate(timeClass = fct_drop(timeClass)) %>% 
  mutate(timeClass = fct_rev(timeClass)) %>%
  group_by(timeClass, SITE_TYPE) %>% 
  summarytools::descr(var = stemDen.ha,
                      stats = "common")

Short Saplings (<=1.5m Tall)) - Combined Winter Range (AC+ANC)

asp.sm.acanc <- asp.sm.acanc %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001) 

# reverse factor levels
# asp.sm.acanc <- asp.sm.acanc %>% 
#   mutate(fenced = fct_rev(fenced))
  # mutate(time_class = fct_rev(time_class))

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.acanc.sm <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.sm.acanc,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.acanc.sm)

## Priors for model 'stmod_stally2.acanc.sm' 
## ------
## 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.27,5.25,6.16,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.acanc.sm)

# extract posteriors
posteriors.acanc.stcnt.sm <- insight::get_parameters(stmod_stally2.acanc.sm) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) %>%
  rename(Fenced = fencedFenced) %>%
  rename('2013:Fenced' = `time_class2013:fencedFenced`) %>%
  rename('2018:Fenced' = `time_class2018:fencedFenced`)

# names(posteriors.acanc.stcnt.sm)

# tidy
posteriors.acanc.stcnt.sm.tidy <- posteriors.acanc.stcnt.sm %>%
  pivot_longer(cols = c(`(Intercept)`,"2013","2018","Fenced","2013:Fenced", "2018:Fenced"),names_to = "parameter", values_to = "values")
  # pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_acanc.stcnt.sm <- mcmc_areas(posteriors.acanc.stcnt.sm,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_acanc.stcnt.sm

Effect Description

fun.sexit.gt(stmod_stally2.acanc.sm)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.42	0.9	1.16	1.67	1	1	1.00
time_class2013	0.51	0.9	0.43	0.58	1	1	1.00
time_class2018	0.22	0.9	0.15	0.30	1	1	0.05
fencedFenced	1.02	0.9	0.48	1.58	1	1	0.98
time_class2013:fencedFenced	-1.04	0.9	-1.18	-0.90	1	1	1.00
time_class2018:fencedFenced	-1.94	0.9	-2.13	-1.76	1	1	1.00

# names(stmod_stally2.acanc.sm)
# stmod_stally2.acanc.sm$family 
# stmod_stally2.acanc.sm$linear.predictors

fun.sexit.gt(stmod_stally2.acanc.sm) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Small Saplings, Combined Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_acanc_smsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
plot.pdif.acanc.sm <- p_direction(stmod_stally2.acanc.sm) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Small Saplings, Combined Winter Range") +
  scale_fill_manual(values = colfunc3(2))

plot.pdif.acanc.sm

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_acanc_smsap.png", dpi = 300, width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.acanc.sm <- rope(stmod_stally2.acanc.sm, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.acanc.sm %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_acanc_smap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.acanc.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Winter Range, Aspen Stem Count, Small Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc_smsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc_smsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.acanc.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  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
BL - 2013	Unfenced	-0.51	-0.60	-0.42
BL - 2018	Unfenced	-0.22	-0.32	-0.13
2013 - 2018	Unfenced	0.28	0.22	0.35
BL - 2013	Fenced	0.53	0.40	0.68
BL - 2018	Fenced	1.72	1.52	1.92
2013 - 2018	Fenced	1.19	0.99	1.39

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.acanc.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	Unfenced	0.60	0.55	0.66
BL / 2018	Unfenced	0.80	0.73	0.87
2013 / 2018	Unfenced	1.33	1.24	1.42
BL / 2013	Fenced	1.70	1.47	1.95
BL / 2018	Fenced	5.59	4.57	6.82
2013 / 2018	Fenced	3.29	2.68	3.98

# 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.

plot.emmip.smsap.acanc <- emmip(stmod_stally2.acanc.sm, fenced ~ time_class, type = "reponse",
  breaks = seq(0.50, 1, by = 0.10)) +
  theme_minimal() +
  scale_color_manual(values = colfunc3(2)) +
  labs(color="") +
  theme(legend.position = c(.15, .15))

plot.emmip.smsap.acanc

Short Saplings (<=1.5m Tall) - Core Winter Range (AC)

asp.sm.ac <- asp.sm.ac %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.ac.sm <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.sm.ac,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.ac.sm)

## Priors for model 'stmod_stally2.ac.sm' 
## ------
## 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.26,5.23,5.42,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.ac.sm)

# extract posteriors
posteriors.ac.stcnt.sm <- insight::get_parameters(stmod_stally2.ac.sm) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) %>%
  rename(Fenced = fencedFenced) %>%
  rename('2013:Fenced' = `time_class2013:fencedFenced`) %>%
  rename('2018:Fenced' = `time_class2018:fencedFenced`)

# names(posteriors.ac.stcnt.sm)

# tidy
posteriors.ac.stcnt.sm.tidy <- posteriors.ac.stcnt.sm %>%
  pivot_longer(cols = c(`(Intercept)`,"2013","2018","Fenced","2013:Fenced", "2018:Fenced"),names_to = "parameter", values_to = "values")

# tidy
# posteriors.ac.stcnt.sm.tidy <- posteriors.ac.stcnt.sm %>%
#   pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_ac.stcnt.sm <- mcmc_areas(posteriors.ac.stcnt.sm,
           prob = 0.9) +
  ggtitle("Core Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ac.stcnt.sm

Effect Description

fun.sexit.gt(stmod_stally2.ac.sm)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.61	0.9	1.27	1.97	1.00	1.00	1.00
time_class2013	0.57	0.9	0.48	0.65	1.00	1.00	1.00
time_class2018	0.09	0.9	0.00	0.19	0.94	0.76	0.00
fencedFenced	0.83	0.9	0.19	1.47	0.98	0.98	0.91
time_class2013:fencedFenced	-1.10	0.9	-1.25	-0.95	1.00	1.00	1.00
time_class2018:fencedFenced	-1.81	0.9	-2.01	-1.62	1.00	1.00	1.00

fun.sexit.gt(stmod_stally2.ac.sm) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Small Saplings, Core Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ac_smsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
p_direction(stmod_stally2.ac.sm) %>% #distinct(Parameter)
  # mutate(Parameter = case_when(Parameter == "time_class2013"~'2013',
  # Parameter == "time_class2018"~'2018',
  # Parameter == "time_class2013:fencedFenced"~'2013:Fenced',
  # Parameter == "time_class2018:fencedFenced"~'2018:Fenced',
  # Parameter == "fencedFenced"~'Fenced')) %>%
  plot() +
  theme_minimal() +
  scale_fill_manual(values = colfunc3(2)) +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Small Saplings, Core Winter Range") +
  scale_fill_manual(values = colfunc3(2))

# p_direction(stmod_stally2.ac.sm) %>%
#   tabyl(Parameter)

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac_smsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac_smsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.ac.sm <- rope(stmod_stally2.ac.sm, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.ac.sm %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ac_smap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.ac.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Aspen Stem Count, Small Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac_smsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac_smsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.ac.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  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
BL - 2013	Unfenced	-0.57	-0.67	-0.46
BL - 2018	Unfenced	-0.09	-0.20	0.02
2013 - 2018	Unfenced	0.48	0.40	0.55
BL - 2013	Fenced	0.53	0.39	0.67
BL - 2018	Fenced	1.72	1.53	1.94
2013 - 2018	Fenced	1.19	0.99	1.39

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.ac.sm, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	Unfenced	0.57	0.51	0.63
BL / 2018	Unfenced	0.91	0.82	1.02
2013 / 2018	Unfenced	1.61	1.49	1.74
BL / 2013	Fenced	1.70	1.47	1.94
BL / 2018	Fenced	5.59	4.55	6.85
2013 / 2018	Fenced	3.29	2.67	3.98

# 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.

Short Saplings (<=1.5m Tall) - Noncore Winter Range (ANC)

asp.sm.anc <- asp.sm.anc %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.anc.sm <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.sm.anc,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.anc.sm)

## Priors for model 'stmod_stally2.anc.sm' 
## ------
## 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.27,5.27])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.anc.sm)

# extract posteriors
posteriors.anc.stcnt.sm <- insight::get_parameters(stmod_stally2.anc.sm)

# tidy
posteriors.anc.stcnt.sm.tidy <- posteriors.anc.stcnt.sm %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")


posteriors.anc.stcnt.sm <- insight::get_parameters(stmod_stally2.anc.sm) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) 

# names(posteriors.anc.stcnt.sm)

# tidy
# posteriors.ac.stcnt.sm.tidy <- posteriors.ac.stcnt.sm %>%
#   pivot_longer(cols = c(`(Intercept)`,"2013","2018","Unfenced","2013:Fenced", "2018:Fenced"),names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_anc.stcnt.sm <- mcmc_areas(posteriors.anc.stcnt.sm,
           prob = 0.9) +
  ggtitle("Noncore Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_anc.stcnt.sm

Effect Description

fun.sexit.gt(stmod_stally2.anc.sm)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.12	0.9	0.84	1.42	1.00	1.00	1.00
time_class2013	0.20	0.9	0.07	0.35	0.99	0.96	0.14
time_class2018	0.58	0.9	0.46	0.72	1.00	1.00	1.00

fun.sexit.gt(stmod_stally2.anc.sm) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Small Saplings, Noncore Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_anc_smsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - anc
p_direction(stmod_stally2.anc.sm) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Small Saplings, Noncore Winter Range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc_smsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc_smsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.anc.sm <- rope(stmod_stally2.anc.sm, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.anc.sm %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_anc_smap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.anc.sm, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Noncore Winter Range, Aspen Stem Count, Small Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc_smsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc_smsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.anc.sm, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  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	estimate	lower.HPD	upper.HPD
BL - 2013	-0.20	-0.37	-0.03
BL - 2018	-0.58	-0.74	-0.43
2013 - 2018	-0.38	-0.53	-0.24

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.anc.sm, ~ time_class) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	0.82	0.68	0.96
BL / 2018	0.56	0.47	0.65
2013 / 2018	0.68	0.59	0.79

Short Saplings (<=1.5m Tall)) - Kawuneeche Valley (AK)

Modeling and Posterior Description

asp.sm.ak <- asp.sm.ak %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

#### main factors + interactions:
stmod_stally2.ak.sm <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.sm.ak,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.ak.sm)

## Priors for model 'stmod_stally2.ak.sm' 
## ------
## 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.27,5.27])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.ak.sm)

# extract posteriors
posteriors.ak.stcnt.sm <- insight::get_parameters(stmod_stally2.ak.sm) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) 

# tidy
posteriors.ak.stcnt.sm.tidy <- posteriors.ak.stcnt.sm %>%
  pivot_longer(cols = c('2013','2018'), names_to = "parameter", values_to = "values")
  # pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

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

plot.mcmcareas_ak.stcnt.sm <- mcmc_areas(posteriors.ak.stcnt.sm,
           prob = 0.9) +
  ggtitle("Kawuneeche Valley") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ak.stcnt.sm

Panel plot of above MCMC area plots

## rendered for print doc
## turned off for Rmd appx
# panel of ACNC AC ANC AK - small saplings

smsap.plot1 <- mcmc_areas(posteriors.acanc.stcnt.sm,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  geom_vline(aes(xintercept=0), lty="dashed") +
  theme(text = element_text(family = 'sans'))

smsap.plot2 <- mcmc_areas(posteriors.ac.stcnt.sm,
           prob = 0.9) +
  ggtitle("Core Winter Range") +
  theme_minimal() +
  geom_vline(aes(xintercept=0), lty="dashed") +
  theme(text = element_text(family = 'sans'))

smsap.plot3 <- mcmc_areas(posteriors.anc.stcnt.sm,
           prob = 0.9) +
  ggtitle("Noncore Winter Range") +
  theme_minimal() +
  geom_vline(aes(xintercept=0), lty="dashed") +
  theme(text = element_text(family = 'sans'))

smsap.plot4 <- mcmc_areas(posteriors.ak.stcnt.sm,
           prob = 0.9) +
  ggtitle("Kawuneeche Valley") +
  theme_minimal() +
  geom_vline(aes(xintercept=0), lty="dashed") +
  theme(text = element_text(family = 'sans'))

# combine
smsap.plot1 + smsap.plot2 + smsap.plot3 + smsap.plot4 + plot_annotation(
  title = 'Posterior Distributions: Small Aspen Stems')

# save plots
ggsave("./output/figures_exported/asp_small_sapling_panel_plot.png", dpi = 300, width = 8.75, height = 4.75)
ggsave("./output/figures_exported/asp_small_sapling_panel_plot.pdf", width = 8.75, height = 4.75)

Effect Description

fun.sexit.gt(stmod_stally2.ak.sm)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.08	0.9	0.66	1.47	1.0	1.00	1.00
time_class2013	0.19	0.9	-0.05	0.45	0.9	0.82	0.24
time_class2018	0.38	0.9	0.14	0.61	1.0	0.99	0.70

fun.sexit.gt(stmod_stally2.ak.sm) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Small Saplings, Kawuneeche Valley") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ak_smsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - ak
p_direction(stmod_stally2.ak.sm) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Small Saplings, Kawuneeche Valley") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak_smsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak_smsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.ak.sm <- rope(stmod_stally2.ak.sm, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.ak.sm %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ak_smap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.ak.sm, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Kawuneeche Valley, Aspen Stem Count, Small Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ak_smsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ak_smsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.ak.sm, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  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	estimate	lower.HPD	upper.HPD
BL - 2013	-0.19	-0.50	0.09
BL - 2018	-0.38	-0.66	-0.09
2013 - 2018	-0.18	-0.46	0.07

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.ak.sm, ~ time_class) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	0.83	0.59	1.08
BL / 2018	0.69	0.50	0.89
2013 / 2018	0.83	0.62	1.06

Effects of Burning on Short Saplings (<1.5m Tall) - Combined Winter Range (AC+ANC)

asp.sm.acanc %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  group_by(time_class, range_type, burned, fenced) %>% 
  summarytools::descr(var = stem_den_ac, stats = "fivenum") %>%
  summarytools::tb() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  gt() %>% 
  tab_header(title="Summary statistics")

Summary statistics
time_class	range_type	burned	fenced	variable	min	q1	med	q3	max
BL	core winter range	Unburned	Unfenced	stem_den_ac	0.0	0.0	161.9	1295.0	8741.2
BL	core winter range	Unburned	Fenced	stem_den_ac	0.0	0.0	809.4	1942.5	21529.3
BL	core winter range	Burned	Unfenced	stem_den_ac	0.0	0.0	161.9	161.9	161.9
BL	core winter range	Burned	Fenced	stem_den_ac	0.0	0.0	566.6	4046.9	4370.6
BL	non-core winter range	Unburned	Unfenced	stem_den_ac	0.0	0.0	323.7	809.4	6798.7
BL	non-core winter range	Burned	Unfenced	stem_den_ac	161.9	161.9	323.7	809.4	809.4
2013	core winter range	Unburned	Unfenced	stem_den_ac	0.0	0.0	485.6	2104.4	20234.3
2013	core winter range	Unburned	Fenced	stem_den_ac	0.0	161.9	485.6	1295.0	12302.5
2013	core winter range	Burned	Unfenced	stem_den_ac	0.0	1375.9	4532.5	18615.6	43058.6
2013	core winter range	Burned	Fenced	stem_den_ac	0.0	647.5	3075.6	3885.0	5341.9
2013	non-core winter range	Unburned	Unfenced	stem_den_ac	0.0	161.9	566.6	1133.1	5341.9
2013	non-core winter range	Burned	Unfenced	stem_den_ac	0.0	0.0	0.0	161.9	161.9
2018	core winter range	Unburned	Unfenced	stem_den_ac	0.0	0.0	485.6	1942.5	12140.6
2018	core winter range	Unburned	Fenced	stem_den_ac	0.0	0.0	161.9	647.5	2751.9
2018	core winter range	Burned	Unfenced	stem_den_ac	0.0	404.7	3642.2	6717.8	9712.5
2018	core winter range	Burned	Fenced	stem_den_ac	0.0	323.7	647.5	1133.1	1456.9
2018	non-core winter range	Unburned	Unfenced	stem_den_ac	0.0	323.7	1052.2	1699.7	4694.4
2018	non-core winter range	Burned	Unfenced	stem_den_ac	161.9	161.9	161.9	323.7	323.7

### STAN model
stmod_smsap_TCxFxB_acanc <- stan_glmer(stem_tally ~ fenced * time_class * burned + (1 | site_id), data = asp.sm.acanc,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )

fun.sexit.gt(stmod_smsap_TCxFxB_acanc)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.32	0.9	1.06	1.57	1.00	1.00	1.00
fencedFenced	1.12	0.9	0.55	1.70	1.00	1.00	0.99
time_class2013	0.33	0.9	0.25	0.41	1.00	1.00	0.70
time_class2018	0.42	0.9	0.34	0.50	1.00	1.00	1.00
burnedBurned	1.08	0.9	0.07	2.00	0.96	0.96	0.90
fencedFenced:time_class2013	-1.16	0.9	-1.32	-0.99	1.00	1.00	1.00
fencedFenced:time_class2018	-2.37	0.9	-2.58	-2.15	1.00	1.00	1.00
fencedFenced:burnedBurned	-1.29	0.9	-2.81	0.19	0.92	0.91	0.86
time_class2013:burnedBurned	0.43	0.9	-0.20	1.02	0.88	0.86	0.64
time_class2018:burnedBurned	-0.75	0.9	-1.34	-0.11	0.97	0.96	0.88
fencedFenced:time_class2013:burnedBurned	0.96	0.9	0.27	1.62	0.99	0.98	0.94
fencedFenced:time_class2018:burnedBurned	1.89	0.9	1.14	2.60	1.00	1.00	1.00

fun.sexit.gt(stmod_smsap_TCxFxB_acanc) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Combined Winter Range, Small Saplings") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/postdesc_stemcnt_smsap_TCxFxB_acanc.rtf")

prior_summary(stmod_smsap_TCxFxB_acanc)
pp_check(stmod_smsap_TCxFxB_acanc)

Probability of Direction

# Visualize the pd
p_direction(stmod_smsap_TCxFxB_acanc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, combined winter range, TCxFxB") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_smsap_TCxFxB_acanc.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_smsap_TCxFxB_acanc.pdf", width = 4.75, height = 3.75)

#### Contrasts
# linear predictors
# emmeans(stmod_smsap_TCxFxB_acanc, ~ time_class | burned | fenced) %>% 
#   pairs(reverse = TRUE) %>% 
#   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, Aspen Stem Count, TCxFxB")

## response scale
emmeans(stmod_smsap_TCxFxB_acanc, ~ time_class | burned | 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 Winter Range, Aspen Stem Count, TCxFxB")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_smsap_ac.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_smsap_ac.png", width = 4.5, height = 3.55) # save plot

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

Pairwise contrasts
Results are given on the response scale
contrast	burned	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unburned	Unfenced	1.3846566	1.2566761	1.5192577
2018 - BL	Unburned	Unfenced	1.5281534	1.3935401	1.6765504
2018 - 2013	Unburned	Unfenced	1.1029472	1.0124093	1.1953320
2013 - BL	Burned	Unfenced	2.1184675	1.0550338	4.4790006
2018 - BL	Burned	Unfenced	0.7223279	0.3611212	1.5339940
2018 - 2013	Burned	Unfenced	0.3406245	0.3009267	0.3885701
2013 - BL	Unburned	Fenced	0.4346064	0.3665325	0.5081867
2018 - BL	Unburned	Fenced	0.1431008	0.1115601	0.1799221
2018 - 2013	Unburned	Fenced	0.3290702	0.2592981	0.4191948
2013 - BL	Burned	Fenced	1.7398958	1.2665184	2.3957595
2018 - BL	Burned	Fenced	0.4497267	0.2919416	0.6739578
2018 - 2013	Burned	Fenced	0.2589152	0.1817707	0.3739754

Tall Saplings (1.5-2.5m Tall) - Combined Winter Range (AC+ANC)

asp.lg.acanc <- asp.lg.acanc %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001) 

# reverse factor levels
# asp.lg.acanc <- asp.lg.acanc %>% 
#   mutate(fenced = fct_rev(fenced))
  # mutate(time_class = fct_rev(time_class))

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.acanc.lg <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.lg.acanc,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.acanc.lg)

## Priors for model 'stmod_stally2.acanc.lg' 
## ------
## 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.26,5.24,6.16,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.acanc.lg)

# extract posteriors
posteriors.acanc.stcnt.lg <- insight::get_parameters(stmod_stally2.acanc.lg)

# tidy
posteriors.acanc.stcnt.lg.tidy <- posteriors.acanc.stcnt.lg %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_acanc.stcnt.lg <- mcmc_areas(posteriors.acanc.stcnt.lg,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_acanc.stcnt.lg

Effect Description

fun.sexit.gt(stmod_stally2.acanc.lg)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-3.77	0.9	-4.70	-2.89	1.00	1.00	1.00
time_class2013	-0.24	0.9	-0.72	0.32	0.77	0.72	0.43
time_class2018	1.58	0.9	1.16	2.01	1.00	1.00	1.00
fencedFenced	0.23	0.9	-1.58	2.04	0.58	0.56	0.47
time_class2013:fencedFenced	4.72	0.9	3.45	6.03	1.00	1.00	1.00
time_class2018:fencedFenced	2.68	0.9	1.47	3.98	1.00	1.00	1.00

fun.sexit.gt(stmod_stally2.acanc.lg) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Large Saplings, Combined Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_acanc_lgsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
p_direction(stmod_stally2.acanc.lg) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Large Saplings, Combined Winter Range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_acanc_lgsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_acanc_lgsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.acanc.lg <- rope(stmod_stally2.acanc.lg, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.acanc.lg %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_acanc_lgap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.acanc.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Winter Range, Aspen Stem Count, Large Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc_lgsap.pdf", width = 4.5, height = 3.55) # save plot

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc_lgsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.acanc.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  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
BL - 2013	Unfenced	0.24	-0.39	0.86
BL - 2018	Unfenced	-1.58	-2.09	-1.06
2013 - 2018	Unfenced	-1.81	-2.25	-1.41
BL - 2013	Fenced	-4.46	-6.04	-3.17
BL - 2018	Fenced	-4.25	-5.82	-2.95
2013 - 2018	Fenced	0.20	-0.03	0.44

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.acanc.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	Unfenced	1.27	0.59	2.18
BL / 2018	Unfenced	0.21	0.11	0.32
2013 / 2018	Unfenced	0.16	0.10	0.24
BL / 2013	Fenced	0.01	0.00	0.03
BL / 2018	Fenced	0.01	0.00	0.04
2013 / 2018	Fenced	1.23	0.96	1.53

# 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.

Tall Saplings (1.5-2.5m Tall) - Core Winter Range (AC)

asp.lg.ac <- asp.lg.ac %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.ac.lg <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.lg.ac,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.ac.lg)

## Priors for model 'stmod_stally2.ac.lg' 
## ------
## 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.25,5.23,5.42,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.ac.lg)

# extract posteriors
posteriors.ac.stcnt.lg <- insight::get_parameters(stmod_stally2.ac.lg) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) %>%
  rename(Fenced = fencedFenced) %>%
  rename('2013:Fenced' = `time_class2013:fencedFenced`) %>%
  rename('2018:Fenced' = `time_class2018:fencedFenced`)

names(posteriors.ac.stcnt.lg)

## [1] "(Intercept)" "2013"        "2018"        "Fenced"      "2013:Fenced"
## [6] "2018:Fenced"

# tidy
posteriors.ac.stcnt.lg.tidy <- posteriors.ac.stcnt.lg %>%
  pivot_longer(cols = c(`(Intercept)`,"2013","2018","Fenced","2013:Fenced", "2018:Fenced"),names_to = "parameter", values_to = "values")

# tidy
# posteriors.ac.stcnt.lg.tidy <- posteriors.ac.stcnt.lg %>%
#   pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_ac.stcnt.lg <- mcmc_areas(posteriors.ac.stcnt.lg,
           prob = 0.9) +
  ggtitle("Core Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ac.stcnt.lg

Effect Description

fun.sexit.gt(stmod_stally2.ac.lg)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-6.06	0.9	-7.92	-4.32	1.00	1.00	1.00
time_class2013	0.26	0.9	-1.41	2.04	0.60	0.58	0.48
time_class2018	3.63	0.9	2.23	5.22	1.00	1.00	1.00
fencedFenced	2.47	0.9	0.05	4.79	0.96	0.96	0.94
time_class2013:fencedFenced	4.27	0.9	2.09	6.35	1.00	1.00	1.00
time_class2018:fencedFenced	0.68	0.9	-1.36	2.52	0.72	0.70	0.63

fun.sexit.gt(stmod_stally2.ac.lg) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Large Saplings, Core Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ac_lgsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
p_direction(stmod_stally2.ac.lg) %>% #distinct(Parameter)
  # mutate(Parameter = case_when(Parameter == "time_class2013"~'2013',
  # Parameter == "time_class2018"~'2018',
  # Parameter == "time_class2013:fencedFenced"~'2013:Fenced',
  # Parameter == "time_class2018:fencedFenced"~'2018:Fenced',
  # Parameter == "fencedFenced"~'Fenced')) %>%
  plot() +
  theme_minimal() +
  scale_fill_manual(values = colfunc3(2)) +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Tall Saplings, Core Winter Range") +
  scale_fill_manual(values = colfunc3(2))

p_direction(stmod_stally2.ac.lg) %>%
  tabyl(Parameter)

##                    Parameter n   percent
##                  (Intercept) 1 0.1666667
##                 fencedFenced 1 0.1666667
##               time_class2013 1 0.1666667
##  time_class2013:fencedFenced 1 0.1666667
##               time_class2018 1 0.1666667
##  time_class2018:fencedFenced 1 0.1666667

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac_lgsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac_lgsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.ac.lg <- rope(stmod_stally2.ac.lg, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.ac.lg %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ac_lgap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.ac.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Aspen Stem Count, Tall Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac_lgsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac_lgsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.ac.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE) %>% 
  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
BL - 2013	Unfenced	-0.26	-2.47	1.69
BL - 2018	Unfenced	-3.63	-5.69	-2.06
2013 - 2018	Unfenced	-3.37	-4.43	-2.45
BL - 2013	Fenced	-4.50	-6.09	-3.18
BL - 2018	Fenced	-4.29	-5.94	-3.01
2013 - 2018	Fenced	0.20	-0.03	0.43

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.ac.lg, ~ time_class | fenced) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	Unfenced	0.77	0.01	3.66
BL / 2018	Unfenced	0.03	0.00	0.09
2013 / 2018	Unfenced	0.03	0.01	0.07
BL / 2013	Fenced	0.01	0.00	0.03
BL / 2018	Fenced	0.01	0.00	0.04
2013 / 2018	Fenced	1.22	0.96	1.52

# 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.

Tall Saplings (1.5-2.5m Tall) - Noncore Winter Range (ANC)

asp.lg.anc <- asp.lg.anc %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

Modeling and Posterior Description

#### main factors + interactions:
stmod_stally2.anc.lg <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.lg.anc,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.anc.lg)

## Priors for model 'stmod_stally2.anc.lg' 
## ------
## 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.26,5.26])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.anc.lg)

# extract posteriors
posteriors.anc.stcnt.lg <- insight::get_parameters(stmod_stally2.anc.lg)

# tidy
posteriors.anc.stcnt.lg.tidy <- posteriors.anc.stcnt.lg %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")


posteriors.anc.stcnt.lg <- insight::get_parameters(stmod_stally2.anc.lg) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) 

# names(posteriors.anc.stcnt.lg)

# tidy
# posteriors.ac.stcnt.lg.tidy <- posteriors.ac.stcnt.lg %>%
#   pivot_longer(cols = c(`(Intercept)`,"2013","2018","Unfenced","2013:Fenced", "2018:Fenced"),names_to = "parameter", values_to = "values")

color_scheme_set("teal")
plot.mcmcareas_anc.stcnt.lg <- mcmc_areas(posteriors.anc.stcnt.lg,
           prob = 0.9) +
  ggtitle("Noncore Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_acanc.stcnt.lg

Effect Description

fun.sexit.gt(stmod_stally2.anc.lg)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-2.52	0.9	-3.74	-1.46	1.00	1.00	1.00
time_class2013	0.21	0.9	-0.35	0.74	0.74	0.68	0.39
time_class2018	0.85	0.9	0.38	1.34	1.00	1.00	0.98

fun.sexit.gt(stmod_stally2.anc.lg) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Tall Saplings, Noncore Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_anc_lgsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - anc
p_direction(stmod_stally2.anc.lg) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Tall Saplings, Noncore Winter Range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc_lgsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc_lgsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.anc.lg <- rope(stmod_stally2.anc.lg, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.anc.lg %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_anc_lgsap.rtf")

Contrasts

#### Contrasts
# EMMs for later fanctor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.anc.lg, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Noncore Winter Range, Aspen Stem Count, Tall Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc_lgsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc_lgsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.anc.lg, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  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	estimate	lower.HPD	upper.HPD
BL - 2013	-0.21	-0.87	0.43
BL - 2018	-0.85	-1.46	-0.32
2013 - 2018	-0.65	-1.19	-0.12

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.anc.lg, ~ time_class) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	0.81	0.36	1.43
BL / 2018	0.43	0.21	0.70
2013 / 2018	0.52	0.28	0.84

Tall Saplings (1.5-2.5m Tall) - Kawuneeche Valley (AK)

Modeling and Posterior Description

asp.lg.ak <- asp.lg.ak %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

#### main factors + interactions:
stmod_stally2.ak.lg <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.lg.ak,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.ak.lg)

## Priors for model 'stmod_stally2.ak.lg' 
## ------
## 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.25,5.25])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.ak.lg)

# extract posteriors
posteriors.ak.stcnt.lg <- insight::get_parameters(stmod_stally2.ak.lg) %>% 
  rename('2013' = time_class2013) %>%
  rename('2018' = time_class2018) 

# tidy
posteriors.ak.stcnt.lg.tidy <- posteriors.ak.stcnt.lg %>%
  pivot_longer(cols = c('2013','2018'), names_to = "parameter", values_to = "values")
  # pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

color_scheme_set("teal")

plot.mcmcareas_ak.stcnt.lg <- mcmc_areas(posteriors.ak.stcnt.lg,
           prob = 0.9) +
  ggtitle("Kawuneeche Valley") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

plot.mcmcareas_ak.stcnt.lg

Panel plot of above MCMC area plots

## rendered for print doc
## turned off for Rmd appx
# panel of ACNC AC ANC AK - tall saplings
# glimpse(posteriors.acanc.stcnt.lg)

lgsap.plot1 <- mcmc_areas(posteriors.acanc.stcnt.lg,
           prob = 0.9) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  geom_vline(aes(xintercept=0), lty="dashed")

lgsap.plot2 <- mcmc_areas(posteriors.ac.stcnt.lg,
           prob = 0.9) +
  ggtitle("Core Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  geom_vline(aes(xintercept=0), lty="dashed")

lgsap.plot3 <- mcmc_areas(posteriors.anc.stcnt.lg,
           prob = 0.9) +
  ggtitle("Noncore Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  geom_vline(aes(xintercept=0), lty="dashed")

lgsap.plot4 <- mcmc_areas(posteriors.ak.stcnt.lg,
           prob = 0.9) +
  ggtitle("Kawuneeche Valley") +
  theme_minimal() +
  theme(text = element_text(family = 'sans')) +
  geom_vline(aes(xintercept=0), lty="dashed")

# combine
lgsap.plot1 + lgsap.plot2 + lgsap.plot3 + lgsap.plot4 + plot_annotation(
  title = 'Posterior Distributions: Tall Aspen Stems')

# save plots
ggsave("./output/figures_exported/asp_tall_sapling_panel_plot.png", dpi = 300, width = 8.75, height = 4.75)
ggsave("./output/figures_exported/asp_tall_sapling_panel_plot.pdf", width = 8.75, height = 4.75)

Effect Description

fun.sexit.gt(stmod_stally2.ak.lg)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-5.24	0.9	-8.77	-2.40	1.00	1.0	1.00
time_class2013	-0.04	0.9	-4.90	4.56	0.51	0.5	0.46
time_class2018	-0.08	0.9	-4.85	4.82	0.51	0.5	0.47

fun.sexit.gt(stmod_stally2.ak.lg) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Tall Saplings, Kawuneeche Valley") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ak_lgsap.rtf")

Probability of Direction

# Visualize the pd
#### PD - ak
p_direction(stmod_stally2.ak.lg) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Tall Saplings, Kawuneeche Valley") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak_lgsap.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak_lgsap.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.ak.lg <- rope(stmod_stally2.ak.lg, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.ak.lg %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ak_lgsap.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.ak.lg, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  # plot(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 = "Kawuneeche Valley, Aspen Stem Count, Tall Saplings. Results are given on the log scale")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ak_lgsap.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ak_lgsap.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.ak.lg, ~ time_class) %>% 
  pairs(reverse = FALSE) %>% 
  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	estimate	lower.HPD	upper.HPD
BL - 2013	0.04	-5.43	6.04
BL - 2018	0.08	-5.47	6.26
2013 - 2018	0.02	-6.85	6.61

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.ak.lg, ~ time_class) %>% 
  pairs(reverse = FALSE, 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
BL / 2013	1.04	0	140.43
BL / 2018	1.09	0	167.57
2013 / 2018	1.02	0	255.15

Effects of Burning on Tall Saplings (1.5-2.5m Tall) - Combined Winter Range (AC+ANC)

asp.lg.acanc %>% 
  filter(range_type == "core winter range") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  group_by(time_class, range_type, burned, fenced) %>% 
  summarytools::descr(stats = "fivenum") %>%
  summarytools::tb() %>% 
  gt() %>% 
  tab_header(title="Summary statistics")

Summary statistics
time_class	range_type	burned	fenced	variable	min	q1	med	q3	max
BL	core winter range	Unburned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Unburned	Unfenced	site_number	2.000000e+00	9.000000e+00	2.000000e+01	3.000000e+01	4.000000e+01
BL	core winter range	Unburned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Unburned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07
BL	core winter range	Unburned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Unburned	Unfenced	sum_stem	0.000000e+00	6.000000e+00	1.400000e+01	2.600000e+01	5.900000e+01
BL	core winter range	Unburned	Unfenced	utm_e_nad83	4.453480e+05	4.464030e+05	4.465360e+05	4.473970e+05	4.505480e+05
BL	core winter range	Unburned	Unfenced	utm_n_nad83	4.467366e+06	4.468377e+06	4.472409e+06	4.473033e+06	4.473650e+06
BL	core winter range	Unburned	Unfenced	year	2.006000e+03	2.006000e+03	2.007000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Unburned	Fenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	5.000000e-02
BL	core winter range	Unburned	Fenced	site_number	8.000000e+00	3.100000e+01	3.950000e+01	6.000000e+01	6.200000e+01
BL	core winter range	Unburned	Fenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.618744e+02
BL	core winter range	Unburned	Fenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	4.000000e+02
BL	core winter range	Unburned	Fenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.000000e+00
BL	core winter range	Unburned	Fenced	sum_stem	0.000000e+00	1.300000e+01	2.450000e+01	4.600000e+01	1.390000e+02
BL	core winter range	Unburned	Fenced	utm_e_nad83	4.466650e+05	4.481500e+05	4.483695e+05	4.495780e+05	4.500990e+05
BL	core winter range	Unburned	Fenced	utm_n_nad83	4.467680e+06	4.468965e+06	4.471327e+06	4.472678e+06	4.472916e+06
BL	core winter range	Unburned	Fenced	year	2.007000e+03	2.007000e+03	2.008000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Burned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	1.666667e-01	3.333333e-01	3.333333e-01
BL	core winter range	Burned	Unfenced	site_number	2.700000e+01	2.700000e+01	2.700000e+01	2.700000e+01	2.700000e+01
BL	core winter range	Burned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	8.093720e+01	1.618744e+02	1.618744e+02
BL	core winter range	Burned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	2.000000e+02	4.000000e+02	4.000000e+02
BL	core winter range	Burned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	5.000000e-01	1.000000e+00	1.000000e+00
BL	core winter range	Burned	Unfenced	sum_stem	3.000000e+00	3.000000e+00	3.000000e+00	3.000000e+00	3.000000e+00
BL	core winter range	Burned	Unfenced	utm_e_nad83	4.466150e+05	4.466150e+05	4.466150e+05	4.466150e+05	4.466150e+05
BL	core winter range	Burned	Unfenced	utm_n_nad83	4.467326e+06	4.467326e+06	4.467326e+06	4.467326e+06	4.467326e+06
BL	core winter range	Burned	Unfenced	year	2.009000e+03	2.009000e+03	2.009000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Burned	Fenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Burned	Fenced	site_number	1.000000e+01	1.000000e+01	1.450000e+01	1.900000e+01	1.900000e+01
BL	core winter range	Burned	Fenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Burned	Fenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07
BL	core winter range	Burned	Fenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
BL	core winter range	Burned	Fenced	sum_stem	2.700000e+01	2.700000e+01	2.950000e+01	3.200000e+01	3.200000e+01
BL	core winter range	Burned	Fenced	utm_e_nad83	4.469290e+05	4.469290e+05	4.470120e+05	4.470950e+05	4.470950e+05
BL	core winter range	Burned	Fenced	utm_n_nad83	4.467463e+06	4.467463e+06	4.467476e+06	4.467489e+06	4.467489e+06
BL	core winter range	Burned	Fenced	year	2.006000e+03	2.006000e+03	2.006000e+03	2.006000e+03	2.006000e+03
2013	core winter range	Unburned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	5.000000e-01
2013	core winter range	Unburned	Unfenced	site_number	2.000000e+00	9.000000e+00	2.000000e+01	3.000000e+01	4.000000e+01
2013	core winter range	Unburned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.618744e+02
2013	core winter range	Unburned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	4.000000e+02
2013	core winter range	Unburned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.000000e+00
2013	core winter range	Unburned	Unfenced	sum_stem	0.000000e+00	7.000000e+00	2.200000e+01	3.400000e+01	1.440000e+02
2013	core winter range	Unburned	Unfenced	utm_e_nad83	4.453480e+05	4.464030e+05	4.465360e+05	4.473970e+05	4.505480e+05
2013	core winter range	Unburned	Unfenced	utm_n_nad83	4.467366e+06	4.468377e+06	4.472409e+06	4.473033e+06	4.473650e+06
2013	core winter range	Unburned	Unfenced	year	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Unburned	Fenced	percent_tot	0.000000e+00	0.000000e+00	1.846591e-01	3.142857e-01	8.000000e-01
2013	core winter range	Unburned	Fenced	site_number	8.000000e+00	3.100000e+01	4.400000e+01	6.100000e+01	6.700000e+01
2013	core winter range	Unburned	Fenced	stem_den_ac	0.000000e+00	0.000000e+00	7.284348e+02	1.456870e+03	4.694358e+03
2013	core winter range	Unburned	Fenced	stem_den_ha	1.000000e-07	1.000000e-07	1.800000e+03	3.600000e+03	1.160000e+04
2013	core winter range	Unburned	Fenced	stem_tally	0.000000e+00	0.000000e+00	4.500000e+00	9.000000e+00	2.900000e+01
2013	core winter range	Unburned	Fenced	sum_stem	4.000000e+00	1.000000e+01	2.800000e+01	7.000000e+01	1.040000e+02
2013	core winter range	Unburned	Fenced	utm_e_nad83	4.466650e+05	4.477950e+05	4.483260e+05	4.495780e+05	4.500990e+05
2013	core winter range	Unburned	Fenced	utm_n_nad83	4.467489e+06	4.468593e+06	4.471308e+06	4.472678e+06	4.472916e+06
2013	core winter range	Unburned	Fenced	year	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Burned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.908397e-03
2013	core winter range	Burned	Unfenced	site_number	2.700000e+01	4.500000e+01	6.350000e+01	6.450000e+01	6.500000e+01
2013	core winter range	Burned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.618744e+02
2013	core winter range	Burned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	4.000000e+02
2013	core winter range	Burned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.000000e+00
2013	core winter range	Burned	Unfenced	sum_stem	1.700000e+01	9.250000e+01	1.685000e+02	3.465000e+02	5.240000e+02
2013	core winter range	Burned	Unfenced	utm_e_nad83	4.465450e+05	4.465800e+05	4.468505e+05	4.471050e+05	4.471240e+05
2013	core winter range	Burned	Unfenced	utm_n_nad83	4.467326e+06	4.467333e+06	4.467378e+06	4.467427e+06	4.467439e+06
2013	core winter range	Burned	Unfenced	year	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Burned	Fenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
2013	core winter range	Burned	Fenced	site_number	1.000000e+01	1.000000e+01	1.900000e+01	6.600000e+01	6.600000e+01
2013	core winter range	Burned	Fenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
2013	core winter range	Burned	Fenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07
2013	core winter range	Burned	Fenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00
2013	core winter range	Burned	Fenced	sum_stem	4.400000e+01	4.400000e+01	5.000000e+01	5.200000e+01	5.200000e+01
2013	core winter range	Burned	Fenced	utm_e_nad83	4.469290e+05	4.469290e+05	4.470050e+05	4.470950e+05	4.470950e+05
2013	core winter range	Burned	Fenced	utm_n_nad83	4.467463e+06	4.467463e+06	4.467489e+06	4.467494e+06	4.467494e+06
2013	core winter range	Burned	Fenced	year	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2018	core winter range	Unburned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	2.295082e-01
2018	core winter range	Unburned	Unfenced	site_number	2.000000e+00	1.000000e+01	2.100000e+01	3.100000e+01	6.800000e+01
2018	core winter range	Unburned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	2.266242e+03
2018	core winter range	Unburned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	1.000000e-07	5.600000e+03
2018	core winter range	Unburned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	1.400000e+01
2018	core winter range	Unburned	Unfenced	sum_stem	0.000000e+00	6.500000e+00	1.800000e+01	4.650000e+01	9.500000e+01
2018	core winter range	Unburned	Unfenced	utm_e_nad83	4.453480e+05	4.464100e+05	4.465395e+05	4.474950e+05	4.505480e+05
2018	core winter range	Unburned	Unfenced	utm_n_nad83	4.467366e+06	4.468958e+06	4.472412e+06	4.473022e+06	4.473650e+06
2018	core winter range	Unburned	Unfenced	year	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Unburned	Fenced	percent_tot	0.000000e+00	0.000000e+00	2.222222e-01	3.333333e-01	6.111111e-01
2018	core winter range	Unburned	Fenced	site_number	8.000000e+00	3.100000e+01	4.400000e+01	6.100000e+01	6.700000e+01
2018	core winter range	Unburned	Fenced	stem_den_ac	0.000000e+00	0.000000e+00	3.237488e+02	1.133121e+03	3.723111e+03
2018	core winter range	Unburned	Fenced	stem_den_ha	1.000000e-07	1.000000e-07	8.000000e+02	2.800000e+03	9.200000e+03
2018	core winter range	Unburned	Fenced	stem_tally	0.000000e+00	0.000000e+00	2.000000e+00	7.000000e+00	2.300000e+01
2018	core winter range	Unburned	Fenced	sum_stem	3.000000e+00	4.000000e+00	9.000000e+00	2.100000e+01	5.700000e+01
2018	core winter range	Unburned	Fenced	utm_e_nad83	4.466650e+05	4.477950e+05	4.483260e+05	4.495780e+05	4.500990e+05
2018	core winter range	Unburned	Fenced	utm_n_nad83	4.467489e+06	4.468593e+06	4.471308e+06	4.472678e+06	4.472916e+06
2018	core winter range	Unburned	Fenced	year	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Burned	Unfenced	percent_tot	0.000000e+00	0.000000e+00	0.000000e+00	1.114394e-01	2.303371e-01
2018	core winter range	Burned	Unfenced	site_number	2.700000e+01	4.500000e+01	6.350000e+01	6.450000e+01	6.500000e+01
2018	core winter range	Burned	Unfenced	stem_den_ac	0.000000e+00	0.000000e+00	0.000000e+00	2.994676e+03	6.636850e+03
2018	core winter range	Burned	Unfenced	stem_den_ha	1.000000e-07	1.000000e-07	1.000000e-07	7.400000e+03	1.640000e+04
2018	core winter range	Burned	Unfenced	stem_tally	0.000000e+00	0.000000e+00	0.000000e+00	1.850000e+01	4.100000e+01
2018	core winter range	Burned	Unfenced	sum_stem	2.000000e+00	3.950000e+01	9.650000e+01	1.470000e+02	1.780000e+02
2018	core winter range	Burned	Unfenced	utm_e_nad83	4.465450e+05	4.465800e+05	4.468505e+05	4.471050e+05	4.471240e+05
2018	core winter range	Burned	Unfenced	utm_n_nad83	4.467326e+06	4.467333e+06	4.467378e+06	4.467427e+06	4.467439e+06
2018	core winter range	Burned	Unfenced	year	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Burned	Fenced	percent_tot	7.142857e-02	1.428571e-01	2.169231e-01	3.333333e-01	4.400000e-01
2018	core winter range	Burned	Fenced	site_number	1.000000e+01	1.000000e+01	1.900000e+01	6.600000e+01	6.600000e+01
2018	core winter range	Burned	Fenced	stem_den_ac	1.618744e+02	3.237488e+02	1.052184e+03	1.780618e+03	2.104367e+03
2018	core winter range	Burned	Fenced	stem_den_ha	4.000000e+02	8.000000e+02	2.600000e+03	4.400000e+03	5.200000e+03
2018	core winter range	Burned	Fenced	stem_tally	1.000000e+00	2.000000e+00	6.500000e+00	1.100000e+01	1.300000e+01
2018	core winter range	Burned	Fenced	sum_stem	1.400000e+01	1.400000e+01	2.500000e+01	3.900000e+01	3.900000e+01
2018	core winter range	Burned	Fenced	utm_e_nad83	4.469290e+05	4.469290e+05	4.470050e+05	4.470950e+05	4.470950e+05
2018	core winter range	Burned	Fenced	utm_n_nad83	4.467463e+06	4.467463e+06	4.467489e+06	4.467494e+06	4.467494e+06
2018	core winter range	Burned	Fenced	year	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03

### STAN model
stmod_lgsap_TCxFxB_acanc <- stan_glmer(stem_tally ~ fenced * time_class * burned + (1 | site_id), data = asp.lg.acanc,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )

fun.sexit.gt(stmod_lgsap_TCxFxB_acanc)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	-4.07	0.9	-5.08	-3.15	1.00	1.00	1.00
fencedFenced	0.56	0.9	-1.42	2.52	0.68	0.67	0.58
time_class2013	0.30	0.9	-0.23	0.81	0.82	0.78	0.50
time_class2018	1.44	0.9	1.00	1.89	1.00	1.00	1.00
burnedBurned	2.32	0.9	-0.64	5.14	0.90	0.89	0.87
fencedFenced:time_class2013	4.23	0.9	2.94	5.65	1.00	1.00	1.00
fencedFenced:time_class2018	2.51	0.9	1.25	3.92	1.00	1.00	1.00
fencedFenced:burnedBurned	-9.16	0.9	-20.68	0.35	0.97	0.97	0.96
time_class2013:burnedBurned	-2.34	0.9	-4.87	0.39	0.92	0.91	0.89
time_class2018:burnedBurned	0.41	0.9	-1.76	2.70	0.63	0.62	0.54
fencedFenced:time_class2013:burnedBurned	-9.32	0.9	-32.35	10.31	0.80	0.80	0.79
fencedFenced:time_class2018:burnedBurned	7.50	0.9	-1.32	18.92	0.94	0.94	0.93

fun.sexit.gt(stmod_lgsap_TCxFxB_acanc) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Combined Winter Range, Tall Saplings") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/postdesc_stemcnt_lgsap_TCxFxB_acanc.rtf")

prior_summary(stmod_lgsap_TCxFxB_acanc)
pp_check(stmod_lgsap_TCxFxB_acanc)

Probability of Direction

# Visualize the pd
p_direction(stmod_lgsap_TCxFxB_acanc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, combined winter range, TCxFxB") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_lgsap_TCxFxB_acanc.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_lgsap_TCxFxB_acanc.pdf", width = 4.75, height = 3.75)

#### Contrasts
# linear predictors
# emmeans(stmod_lgsap_TCxFxB_acanc, ~ time_class | burned | fenced) %>% 
#   pairs(reverse = TRUE) %>% 
#   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, Aspen Stem Count, TCxFxB")

## response scale
emmeans(stmod_lgsap_TCxFxB_acanc, ~ time_class | burned | 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 Winter Range, Large Saplings, TCxFxB")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_lgsap_ac.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_lgsap_ac.png", width = 4.5, height = 3.55) # save plot

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

Pairwise contrasts
Results are given on the response scale
contrast	burned	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unburned	Unfenced	1.352523e+00	7.261301e-01	2.504175e+00
2018 - BL	Unburned	Unfenced	4.223233e+00	2.461819e+00	7.119706e+00
2018 - 2013	Unburned	Unfenced	3.131594e+00	1.942904e+00	4.925634e+00
2013 - BL	Burned	Unfenced	1.308466e-01	7.081628e-03	4.027584e+00
2018 - BL	Burned	Unfenced	6.338827e+00	5.962042e-01	1.265609e+02
2018 - 2013	Burned	Unfenced	4.803957e+01	1.274348e+01	2.695391e+02
2013 - BL	Unburned	Fenced	9.106910e+01	2.434602e+01	4.919535e+02
2018 - BL	Unburned	Fenced	5.180530e+01	1.346572e+01	2.764413e+02
2018 - 2013	Unburned	Fenced	5.641663e-01	4.386050e-01	7.338011e-01
2013 - BL	Burned	Fenced	1.115404e-03	3.311301e-16	1.143993e+07
2018 - BL	Burned	Fenced	1.647177e+05	1.828172e+01	1.562622e+11
2018 - 2013	Burned	Fenced	1.838080e+08	2.599885e+01	5.353741e+20

All Stem Diameters - Combined Winter Range (AC+ANC)

# purrr plots
asp.acanc  %>% 
  group_nest(range_type) %>%
  mutate(plots = map2(.y = range_type, .x = data, ~{ggplot(data = .x) +
      geom_boxplot(aes(y = stem_den_ac, x = time_class), fill="Gray80") +
      geom_jitter(aes(y = stem_den_ac, x = time_class), fill="blue", alpha=0.25) +
      labs(title = paste0("Range type: ", .y), x = "", y = "Stem density (stems/acre)") +
      theme_minimal()})) %>% 
  pull(plots)

## [[1]]

## 
## [[2]]

Modeling and Posterior Description

#### Main Factors Only
## Not run. Used in model comparisons, but model w/ interaction selected
asp.acanc <- asp.acanc %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

## poisson
stmod_stally1.acanc <- stan_glmer(stem_tally ~ time_class + fenced + (1 | site_id), data = asp.acanc,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )
prior_summary(stmod_stally1.acanc)
pp_check(stmod_stally1.acanc)

#### main factors + interactions:
stmod_stally2.acanc <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.acanc,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2.acanc)

## Priors for model 'stmod_stally2.acanc' 
## ------
## 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.27,5.25,6.16,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2.acanc)

# extract posteriors
posteriors.acanc.stcnt <- insight::get_parameters(stmod_stally2.acanc)

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

library("bayesplot")
color_scheme_set("teal")
mcmc_areas(posteriors.acanc.stcnt,
           prob = 0.8) +
  ggtitle("Combined Winter Range") +
  theme_minimal() +
  theme(text = element_text(family = 'sans'))

Effect Description

fun.sexit.gt(stmod_stally2.acanc)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	0.94	0.9	0.68	1.19	1	1	1.00
time_class2013	0.49	0.9	0.42	0.56	1	1	1.00
time_class2018	0.29	0.9	0.22	0.37	1	1	0.45
fencedFenced	0.92	0.9	0.39	1.43	1	1	0.97
time_class2013:fencedFenced	-0.58	0.9	-0.71	-0.46	1	1	1.00
time_class2018:fencedFenced	-1.19	0.9	-1.34	-1.04	1	1	1.00

fun.sexit.gt(stmod_stally2.acanc) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Core Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ac.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
p_direction(stmod_stally2.acanc) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Combined Winter Range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_acanc.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_acanc.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

rope.stemcnt.ac <- rope(stmod_stally2.acanc, ci=0.9) %>% 
  gt() %>%
  tab_header(title = "Percent in ROPE") %>% 
  fmt_number(
    columns = 5,
    decimals = 1,
    suffixing = TRUE
  ) 

rope.stemcnt.ac %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_acanc.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2.acanc, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Combined Winter Range, Aspen Stem Count, Results are given on the log scale")

# ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc.pdf", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_acanc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2.acanc, ~ 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
2013 - BL	Unfenced	0.49	0.40	0.58
2018 - BL	Unfenced	0.29	0.21	0.39
2018 - 2013	Unfenced	-0.20	-0.26	-0.13
2013 - BL	Fenced	-0.09	-0.22	0.03
2018 - BL	Fenced	-0.90	-1.05	-0.75
2018 - 2013	Fenced	-0.81	-0.96	-0.66

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2.acanc, ~ 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.63	1.50	1.78
2018 / BL	Unfenced	1.34	1.23	1.47
2018 / 2013	Unfenced	0.82	0.77	0.88
2013 / BL	Fenced	0.91	0.80	1.03
2018 / BL	Fenced	0.41	0.35	0.47
2018 / 2013	Fenced	0.45	0.38	0.52

# 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.

All Stem Diameters - Core Winter Range (AC)

Modeling and Posterior Description

#### Main Factors Only
## Not run. Used in model comparisons, but model w/ interaction selected
asp.ac <- asp.ac %>%
  mutate(stem_den_ha = stem_den_ha + 0.0000001)

## poisson
stmod_stally1 <- stan_glmer(stem_tally ~ time_class + fenced + (1 | site_id), data = asp.ac,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )
prior_summary(stmod_stally1)
pp_check(stmod_stally1)

#### main factors + interactions:
stmod_stally2 <- stan_glmer(stem_tally ~ time_class * fenced + (1 | site_id), data = asp.ac,
                      family= poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_stally2)

## Priors for model 'stmod_stally2' 
## ------
## 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.26,5.23,5.42,...])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_stally2)

# extract posteriors
posteriors.ac.stcnt <- insight::get_parameters(stmod_stally2)

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

### ppc plot
ppc.plot(stmod_stally2) +
  labs(caption = "AC: stem_cnt~TCxF")

# save plots
ggsave("./output/figures_exported/aspen_stem_count/ppc_stemcnt_ac.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/ppc_stemcnt_ac.pdf", width = 4.75, height = 3.75)

## model comparison. Run previously for model comparison;
# not currently set to run to minimize computational overhead

loo1 <- loo(stmod_stally1,
            k_threshold = 0.7) 
loo2 <- loo(stmod_stally2,
            k_threshold = 0.7) 
comp <- loo_compare(loo2, loo1)

## create table of comparisons
modcomp_ac <- print(comp, simplify = TRUE, digits = 2)

modcomp_ac %>% 
  gt()
# citation(package = "loo")

# show more details with simplify=FALSE
print(modcomp_ac, simplify = FALSE, digits = 3)

# model with interactions selected

# library("bayesplot")
color_scheme_set("teal")
mcmc_areas(stmod_stally2,
           prob = 0.8) +
  ggtitle("Posterior distributions",
                      "with medians and 80% credible intervals")

Effect Description

fun.sexit.gt(stmod_stally2)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	1.12	0.9	0.76	1.45	1.00	1.00	1.00
time_class2013	0.55	0.9	0.46	0.64	1.00	1.00	1.00
time_class2018	0.18	0.9	0.08	0.27	1.00	0.99	0.01
fencedFenced	0.73	0.9	0.12	1.36	0.97	0.97	0.88
time_class2013:fencedFenced	-0.64	0.9	-0.78	-0.51	1.00	1.00	1.00
time_class2018:fencedFenced	-1.08	0.9	-1.23	-0.92	1.00	1.00	1.00

fun.sexit.gt(stmod_stally2) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Core Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_ac.rtf")

Probability of Direction

# Visualize the pd
#### PD - AC
p_direction(stmod_stally2) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, core winter range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ac.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

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

rope.stemcnt.ac %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ac.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

# library(emmeans)
emmeans(stmod_stally2, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  # labs(title= "Core winter range", caption = "Point estimate displayed: median 
# HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "Contrast")
  labs(subtitle = "Point estimate displayed: median HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "Core Winter Range, Aspen Stem Count, Results are given on the log scale")

# ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac.pdf", width = 4.5, height = 3.55) # save plot
# ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_ac.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_stally2, ~ 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
2013 - BL	Unfenced	0.55	0.45	0.66
2018 - BL	Unfenced	0.18	0.07	0.29
2018 - 2013	Unfenced	-0.37	-0.45	-0.30
2013 - BL	Fenced	-0.09	-0.22	0.03
2018 - BL	Fenced	-0.90	-1.06	-0.75
2018 - 2013	Fenced	-0.81	-0.95	-0.66

# table of contrasts. Results are given on the response scale.
emmeans(stmod_stally2, ~ 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.73	1.56	1.92
2018 / BL	Unfenced	1.19	1.07	1.33
2018 / 2013	Unfenced	0.69	0.64	0.74
2013 / BL	Fenced	0.91	0.80	1.03
2018 / BL	Fenced	0.41	0.34	0.47
2018 / 2013	Fenced	0.45	0.38	0.51

Effects of Burning on All Stems - Core Winter Range (AC)

ht_perc1 %>% 
  filter(range_type == "core winter range") %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  group_by(time_class, range_type, burned, fenced) %>% 
  summarytools::descr(stats = "fivenum") %>%
  summarytools::tb() %>% 
  gt() %>% 
  tab_header(title="Summary statistics")

Summary statistics
time_class	range_type	burned	fenced	variable	min	q1	med	q3	max
BL	core winter range	Unburned	Unfenced	percent_tot	0	0.000000e+00	0.000000e+00	3.333333e-01	1.000000e+00
BL	core winter range	Unburned	Unfenced	site_number	2	9.000000e+00	2.000000e+01	3.000000e+01	4.000000e+01
BL	core winter range	Unburned	Unfenced	stem_den_ac	0	0.000000e+00	0.000000e+00	3.237488e+02	8.741218e+03
BL	core winter range	Unburned	Unfenced	stem_den_ha	0	0.000000e+00	0.000000e+00	8.000000e+02	2.160000e+04
BL	core winter range	Unburned	Unfenced	stem_tally	0	0.000000e+00	0.000000e+00	2.000000e+00	5.400000e+01
BL	core winter range	Unburned	Unfenced	sum_stem	0	6.000000e+00	1.400000e+01	2.600000e+01	5.900000e+01
BL	core winter range	Unburned	Unfenced	utm_e_nad83	445348	4.464030e+05	4.465360e+05	4.473970e+05	4.505480e+05
BL	core winter range	Unburned	Unfenced	utm_n_nad83	4467366	4.468377e+06	4.472409e+06	4.473033e+06	4.473650e+06
BL	core winter range	Unburned	Unfenced	year	2006	2.006000e+03	2.007000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Unburned	Fenced	percent_tot	0	0.000000e+00	3.597122e-02	3.333333e-01	1.000000e+00
BL	core winter range	Unburned	Fenced	site_number	8	3.100000e+01	3.950000e+01	6.000000e+01	6.200000e+01
BL	core winter range	Unburned	Fenced	stem_den_ac	0	0.000000e+00	0.000000e+00	9.712464e+02	2.152930e+04
BL	core winter range	Unburned	Fenced	stem_den_ha	0	0.000000e+00	0.000000e+00	2.400000e+03	5.320000e+04
BL	core winter range	Unburned	Fenced	stem_tally	0	0.000000e+00	0.000000e+00	6.000000e+00	1.330000e+02
BL	core winter range	Unburned	Fenced	sum_stem	0	1.300000e+01	2.450000e+01	4.600000e+01	1.390000e+02
BL	core winter range	Unburned	Fenced	utm_e_nad83	446665	4.481500e+05	4.483695e+05	4.495780e+05	4.500990e+05
BL	core winter range	Unburned	Fenced	utm_n_nad83	4467680	4.468965e+06	4.471327e+06	4.472678e+06	4.472916e+06
BL	core winter range	Unburned	Fenced	year	2007	2.007000e+03	2.008000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Burned	Unfenced	percent_tot	0	0.000000e+00	3.333333e-01	3.333333e-01	3.333333e-01
BL	core winter range	Burned	Unfenced	site_number	27	2.700000e+01	2.700000e+01	2.700000e+01	2.700000e+01
BL	core winter range	Burned	Unfenced	stem_den_ac	0	0.000000e+00	1.618744e+02	1.618744e+02	1.618744e+02
BL	core winter range	Burned	Unfenced	stem_den_ha	0	0.000000e+00	4.000000e+02	4.000000e+02	4.000000e+02
BL	core winter range	Burned	Unfenced	stem_tally	0	0.000000e+00	1.000000e+00	1.000000e+00	1.000000e+00
BL	core winter range	Burned	Unfenced	sum_stem	3	3.000000e+00	3.000000e+00	3.000000e+00	3.000000e+00
BL	core winter range	Burned	Unfenced	utm_e_nad83	446615	4.466150e+05	4.466150e+05	4.466150e+05	4.466150e+05
BL	core winter range	Burned	Unfenced	utm_n_nad83	4467326	4.467326e+06	4.467326e+06	4.467326e+06	4.467326e+06
BL	core winter range	Burned	Unfenced	year	2009	2.009000e+03	2.009000e+03	2.009000e+03	2.009000e+03
BL	core winter range	Burned	Fenced	percent_tot	0	0.000000e+00	0.000000e+00	2.187500e-01	1.000000e+00
BL	core winter range	Burned	Fenced	site_number	10	1.000000e+01	1.450000e+01	1.900000e+01	1.900000e+01
BL	core winter range	Burned	Fenced	stem_den_ac	0	0.000000e+00	0.000000e+00	1.133121e+03	4.370609e+03
BL	core winter range	Burned	Fenced	stem_den_ha	0	0.000000e+00	0.000000e+00	2.800000e+03	1.080000e+04
BL	core winter range	Burned	Fenced	stem_tally	0	0.000000e+00	0.000000e+00	7.000000e+00	2.700000e+01
BL	core winter range	Burned	Fenced	sum_stem	27	2.700000e+01	2.950000e+01	3.200000e+01	3.200000e+01
BL	core winter range	Burned	Fenced	utm_e_nad83	446929	4.469290e+05	4.470120e+05	4.470950e+05	4.470950e+05
BL	core winter range	Burned	Fenced	utm_n_nad83	4467463	4.467463e+06	4.467476e+06	4.467489e+06	4.467489e+06
BL	core winter range	Burned	Fenced	year	2006	2.006000e+03	2.006000e+03	2.006000e+03	2.006000e+03
2013	core winter range	Unburned	Unfenced	percent_tot	0	0.000000e+00	0.000000e+00	3.636364e-01	1.000000e+00
2013	core winter range	Unburned	Unfenced	site_number	2	9.000000e+00	2.000000e+01	3.000000e+01	4.000000e+01
2013	core winter range	Unburned	Unfenced	stem_den_ac	0	0.000000e+00	0.000000e+00	6.474976e+02	2.023430e+04
2013	core winter range	Unburned	Unfenced	stem_den_ha	0	0.000000e+00	0.000000e+00	1.600000e+03	5.000000e+04
2013	core winter range	Unburned	Unfenced	stem_tally	0	0.000000e+00	0.000000e+00	4.000000e+00	1.250000e+02
2013	core winter range	Unburned	Unfenced	sum_stem	0	7.000000e+00	2.200000e+01	3.400000e+01	1.440000e+02
2013	core winter range	Unburned	Unfenced	utm_e_nad83	445348	4.464030e+05	4.465360e+05	4.473970e+05	4.505480e+05
2013	core winter range	Unburned	Unfenced	utm_n_nad83	4467366	4.468377e+06	4.472409e+06	4.473033e+06	4.473650e+06
2013	core winter range	Unburned	Unfenced	year	2013	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Unburned	Fenced	percent_tot	0	2.884615e-02	1.666667e-01	2.788462e-01	9.411765e-01
2013	core winter range	Unburned	Fenced	site_number	8	3.100000e+01	4.400000e+01	6.100000e+01	6.700000e+01
2013	core winter range	Unburned	Fenced	stem_den_ac	0	1.618744e+02	4.856232e+02	1.456870e+03	1.230245e+04
2013	core winter range	Unburned	Fenced	stem_den_ha	0	4.000000e+02	1.200000e+03	3.600000e+03	3.040000e+04
2013	core winter range	Unburned	Fenced	stem_tally	0	1.000000e+00	3.000000e+00	9.000000e+00	7.600000e+01
2013	core winter range	Unburned	Fenced	sum_stem	4	1.000000e+01	2.800000e+01	7.000000e+01	1.040000e+02
2013	core winter range	Unburned	Fenced	utm_e_nad83	446665	4.477950e+05	4.483260e+05	4.495780e+05	4.500990e+05
2013	core winter range	Unburned	Fenced	utm_n_nad83	4467489	4.468593e+06	4.471308e+06	4.472678e+06	4.472916e+06
2013	core winter range	Unburned	Fenced	year	2013	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Burned	Unfenced	percent_tot	0	0.000000e+00	3.721374e-02	4.279044e-01	8.928571e-01
2013	core winter range	Burned	Unfenced	site_number	27	4.500000e+01	6.350000e+01	6.450000e+01	6.500000e+01
2013	core winter range	Burned	Unfenced	stem_den_ac	0	0.000000e+00	2.428116e+02	9.064966e+03	4.305859e+04
2013	core winter range	Burned	Unfenced	stem_den_ha	0	0.000000e+00	6.000000e+02	2.240000e+04	1.064000e+05
2013	core winter range	Burned	Unfenced	stem_tally	0	0.000000e+00	1.500000e+00	5.600000e+01	2.660000e+02
2013	core winter range	Burned	Unfenced	sum_stem	17	9.250000e+01	1.685000e+02	3.465000e+02	5.240000e+02
2013	core winter range	Burned	Unfenced	utm_e_nad83	446545	4.465800e+05	4.468505e+05	4.471050e+05	4.471240e+05
2013	core winter range	Burned	Unfenced	utm_n_nad83	4467326	4.467333e+06	4.467378e+06	4.467427e+06	4.467439e+06
2013	core winter range	Burned	Unfenced	year	2013	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2013	core winter range	Burned	Fenced	percent_tot	0	0.000000e+00	0.000000e+00	4.800000e-01	6.346154e-01
2013	core winter range	Burned	Fenced	site_number	10	1.000000e+01	1.900000e+01	6.600000e+01	6.600000e+01
2013	core winter range	Burned	Fenced	stem_den_ac	0	0.000000e+00	0.000000e+00	3.723111e+03	5.341855e+03
2013	core winter range	Burned	Fenced	stem_den_ha	0	0.000000e+00	0.000000e+00	9.200000e+03	1.320000e+04
2013	core winter range	Burned	Fenced	stem_tally	0	0.000000e+00	0.000000e+00	2.300000e+01	3.300000e+01
2013	core winter range	Burned	Fenced	sum_stem	44	4.400000e+01	5.000000e+01	5.200000e+01	5.200000e+01
2013	core winter range	Burned	Fenced	utm_e_nad83	446929	4.469290e+05	4.470050e+05	4.470950e+05	4.470950e+05
2013	core winter range	Burned	Fenced	utm_n_nad83	4467463	4.467463e+06	4.467489e+06	4.467494e+06	4.467494e+06
2013	core winter range	Burned	Fenced	year	2013	2.013000e+03	2.013000e+03	2.013000e+03	2.013000e+03
2018	core winter range	Unburned	Unfenced	percent_tot	0	0.000000e+00	5.681818e-03	3.809524e-01	1.000000e+00
2018	core winter range	Unburned	Unfenced	site_number	2	1.000000e+01	2.100000e+01	3.100000e+01	6.800000e+01
2018	core winter range	Unburned	Unfenced	stem_den_ac	0	0.000000e+00	0.000000e+00	9.712464e+02	1.214058e+04
2018	core winter range	Unburned	Unfenced	stem_den_ha	0	0.000000e+00	0.000000e+00	2.400000e+03	3.000000e+04
2018	core winter range	Unburned	Unfenced	stem_tally	0	0.000000e+00	0.000000e+00	6.000000e+00	7.500000e+01
2018	core winter range	Unburned	Unfenced	sum_stem	0	6.500000e+00	1.800000e+01	4.650000e+01	9.500000e+01
2018	core winter range	Unburned	Unfenced	utm_e_nad83	445348	4.464100e+05	4.465395e+05	4.474950e+05	4.505480e+05
2018	core winter range	Unburned	Unfenced	utm_n_nad83	4467366	4.468958e+06	4.472412e+06	4.473022e+06	4.473650e+06
2018	core winter range	Unburned	Unfenced	year	2018	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Unburned	Fenced	percent_tot	0	0.000000e+00	2.000000e-01	3.333333e-01	6.666667e-01
2018	core winter range	Unburned	Fenced	site_number	8	3.100000e+01	4.400000e+01	6.100000e+01	6.700000e+01
2018	core winter range	Unburned	Fenced	stem_den_ac	0	0.000000e+00	3.237488e+02	6.474976e+02	3.723111e+03
2018	core winter range	Unburned	Fenced	stem_den_ha	0	0.000000e+00	8.000000e+02	1.600000e+03	9.200000e+03
2018	core winter range	Unburned	Fenced	stem_tally	0	0.000000e+00	2.000000e+00	4.000000e+00	2.300000e+01
2018	core winter range	Unburned	Fenced	sum_stem	3	4.000000e+00	9.000000e+00	2.100000e+01	5.700000e+01
2018	core winter range	Unburned	Fenced	utm_e_nad83	446665	4.477950e+05	4.483260e+05	4.495780e+05	4.500990e+05
2018	core winter range	Unburned	Fenced	utm_n_nad83	4467489	4.468593e+06	4.471308e+06	4.472678e+06	4.472916e+06
2018	core winter range	Unburned	Fenced	year	2018	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Burned	Unfenced	percent_tot	0	0.000000e+00	1.329911e-01	4.083072e-01	5.086207e-01
2018	core winter range	Burned	Unfenced	site_number	27	4.500000e+01	6.350000e+01	6.450000e+01	6.500000e+01
2018	core winter range	Burned	Unfenced	stem_den_ac	0	0.000000e+00	7.284348e+02	6.394039e+03	9.712464e+03
2018	core winter range	Burned	Unfenced	stem_den_ha	0	0.000000e+00	1.800000e+03	1.580000e+04	2.400000e+04
2018	core winter range	Burned	Unfenced	stem_tally	0	0.000000e+00	4.500000e+00	3.950000e+01	6.000000e+01
2018	core winter range	Burned	Unfenced	sum_stem	2	3.950000e+01	9.650000e+01	1.470000e+02	1.780000e+02
2018	core winter range	Burned	Unfenced	utm_e_nad83	446545	4.465800e+05	4.468505e+05	4.471050e+05	4.471240e+05
2018	core winter range	Burned	Unfenced	utm_n_nad83	4467326	4.467333e+06	4.467378e+06	4.467427e+06	4.467439e+06
2018	core winter range	Burned	Unfenced	year	2018	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03
2018	core winter range	Burned	Fenced	percent_tot	0	7.692308e-02	2.000000e-01	2.857143e-01	5.000000e-01
2018	core winter range	Burned	Fenced	site_number	10	1.000000e+01	1.900000e+01	6.600000e+01	6.600000e+01
2018	core winter range	Burned	Fenced	stem_den_ac	0	3.237488e+02	8.093720e+02	1.294995e+03	2.104367e+03
2018	core winter range	Burned	Fenced	stem_den_ha	0	8.000000e+02	2.000000e+03	3.200000e+03	5.200000e+03
2018	core winter range	Burned	Fenced	stem_tally	0	2.000000e+00	5.000000e+00	8.000000e+00	1.300000e+01
2018	core winter range	Burned	Fenced	sum_stem	14	1.400000e+01	2.500000e+01	3.900000e+01	3.900000e+01
2018	core winter range	Burned	Fenced	utm_e_nad83	446929	4.469290e+05	4.470050e+05	4.470950e+05	4.470950e+05
2018	core winter range	Burned	Fenced	utm_n_nad83	4467463	4.467463e+06	4.467489e+06	4.467494e+06	4.467494e+06
2018	core winter range	Burned	Fenced	year	2018	2.018000e+03	2.018000e+03	2.018000e+03	2.018000e+03

### STAN model
stmod_stally_TCxFxB_ac <- stan_glmer(stem_tally ~ fenced * time_class * burned + (1 | site_id), data = asp.ac,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )

fun.sexit.gt(stmod_stally_TCxFxB_ac)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	0.92	0.9	0.58	1.27	1.00	1.00	1.00
fencedFenced	0.92	0.9	0.30	1.56	0.99	0.99	0.95
time_class2013	0.37	0.9	0.28	0.47	1.00	1.00	0.89
time_class2018	0.37	0.9	0.28	0.46	1.00	1.00	0.90
burnedBurned	0.57	0.9	-0.80	1.95	0.75	0.73	0.63
fencedFenced:time_class2013	-0.61	0.9	-0.76	-0.46	1.00	1.00	1.00
fencedFenced:time_class2018	-1.48	0.9	-1.66	-1.31	1.00	1.00	1.00
fencedFenced:burnedBurned	-0.76	0.9	-2.52	1.05	0.76	0.75	0.67
time_class2013:burnedBurned	1.32	0.9	0.31	2.39	0.99	0.99	0.97
time_class2018:burnedBurned	0.46	0.9	-0.53	1.54	0.78	0.75	0.60
fencedFenced:time_class2013:burnedBurned	-0.48	0.9	-1.63	0.54	0.78	0.75	0.61
fencedFenced:time_class2018:burnedBurned	0.62	0.9	-0.45	1.74	0.81	0.80	0.68

fun.sexit.gt(stmod_stally_TCxFxB_ac) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Core Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxFxB_ac.rtf")

prior_summary(stmod_stally_TCxFxB_ac)
pp_check(stmod_stally_TCxFxB_ac)

# Visualize the pd
p_direction(stmod_stally_TCxFxB_ac) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, core winter range, TCxFxB")

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_TCxFxB_ac.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_TCxFxB_ac.pdf", width = 4.75, height = 3.75)

#### Contrasts
# linear predictors
# emmeans(stmod_stally_TCxFxB_ac, ~ time_class | burned | fenced) %>% 
#   pairs(reverse = TRUE) %>% 
#   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, Aspen Stem Count, TCxFxB")

## response scale
emmeans(stmod_stally_TCxFxB_ac, ~ time_class | burned | 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, Aspen Stem Count, TCxFxB")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_ac.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxFxB_ac.png", width = 4.5, height = 3.55) # save plot

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

Pairwise contrasts
Results are given on the response scale
contrast	burned	fenced	estimate	lower.HPD	upper.HPD
2013 - BL	Unburned	Unfenced	1.4481053	1.2944900	1.6177177
2018 - BL	Unburned	Unfenced	1.4496198	1.2932847	1.6129241
2018 - 2013	Unburned	Unfenced	1.0016386	0.9067273	1.1104585
2013 - BL	Burned	Unfenced	5.4275555	1.6872482	20.3252189
2018 - BL	Burned	Unfenced	2.2955855	0.7465244	9.0412182
2018 - 2013	Burned	Unfenced	0.4250671	0.3740706	0.4779338
2013 - BL	Unburned	Fenced	0.7907530	0.6920066	0.9123482
2018 - BL	Unburned	Fenced	0.3287545	0.2745928	0.3912501
2018 - 2013	Unburned	Fenced	0.4158567	0.3484068	0.4959015
2013 - BL	Burned	Fenced	1.8364509	1.3582380	2.5291972
2018 - BL	Burned	Fenced	0.9786425	0.6978295	1.3925061
2018 - 2013	Burned	Fenced	0.5331809	0.4076841	0.7058706

All Stem Diameters - Noncore Winter Range (ANC)

Modeling and Posterior Description

stmod_anc_stally1 <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.anc,
                      family=poisson,
                      iter = 8000,
                      seed = 1234
                      )

# get posteriors
posteriors.anc.stmcnt <- insight::get_parameters(stmod_anc_stally1)
# tidy
posteriors.anc.stcnt.tidy <- posteriors.anc.stmcnt %>%
  pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")

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

prior_summary(stmod_anc_stally1)

## Priors for model 'stmod_anc_stally1' 
## ------
## 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.27,5.27])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_anc_stally1)

#### PP check plot
ppc.plot(stmod_anc_stally1) +
  labs(caption = "ANC: stem_tally~TC")

ggsave("./output/figures_exported/aspen_stem_count/ppc_willowht_anc.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/ppc_willowht_anc.pdf", width = 4.75, height = 3.75)

# library("bayesplot")
color_scheme_set("teal")
mcmc_areas(stmod_anc_stally1,
           prob = 0.8) +
  ggtitle("Posterior distributions",
                      "with medians and 80% credible intervals")

Effect Description

fun.sexit.gt(stmod_anc_stally1)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	0.67	0.9	0.38	0.95	1.00	1.00	0.98
time_class2013	0.20	0.9	0.06	0.34	0.99	0.97	0.12
time_class2018	0.60	0.9	0.48	0.73	1.00	1.00	1.00

fun.sexit.gt(stmod_anc_stally1) %>%
 tab_header(title = "Aspen Count Posterior Description", subtitle = "Noncore Winter Range") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_anc.rtf")

Probability of Direction

# Visualize the pd
#### PD - ANC
p_direction(stmod_anc_stally1) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Noncore winter range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_anc.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

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

rope.stemcnt.ac %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_anc.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

emmeans(stmod_anc_stally1, ~ 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 = "Noncore Winter Range, Aspen Stem Count")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TCxF_anc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_anc_stally1, ~ 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 (not the response) scale.")

Pairwise contrasts
Results are given on the log (not the response) scale.
contrast	estimate	lower.HPD	upper.HPD
2013 - BL	0.2	0.04	0.37
2018 - BL	0.6	0.45	0.75
2018 - 2013	0.4	0.26	0.54

# table of contrasts. Results are given on the response scale.
emmeans(stmod_anc_stally1, ~ 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.23	1.04	1.44
2018 / BL	1.83	1.55	2.11
2018 / 2013	1.49	1.29	1.71

# For 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.

# comment 2.2 "provide analyses on single graphic...
## Plot combinations
contr.plot.ac <- emmeans(stmod_stally2, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  # plot(type = "response") +
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(title= "Core winter range", caption = "", x = "Estimate", y = "Contrast")

contr.plot.anc <- emmeans(stmod_anc_stally1, ~ time_class) %>% 
  pairs(reverse = TRUE) %>% 
  plot(type = "response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(title= "Noncore winter range", caption = "Point estimate displayed: median 
HPD interval probability: 0.95, Results on log scale", x = "Estimate", y = "")

# patchwork
contr.plot.ac + contr.plot.anc

ggsave("./output/figures_exported/revised_figs_aspen/emm_stemcnt_AC_ANC.png", width = 8, height = 3.75, dpi=300)
ggsave("./output/figures_exported/revised_figs_aspen/emm_stemcnt_AC_ANC.pdf", width = 8, height = 3.75)

All Stem Diameters - Kawuneeche Valley (AK)

Modeling and Posterior Description

stmod_ak_stally1 <- stan_glmer(stem_tally ~ time_class + (1 | site_id), data = asp.ak,
                       family=poisson,
                      iter = 8000,
                      seed = 1234
                      )

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

prior_summary(stmod_ak_stally1)

## Priors for model 'stmod_ak_stally1' 
## ------
## 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.28,5.28])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

## model summary
fun.model.summary(stmod_ak_stally1)

# extract posteriors
posteriors.ak.stcnt <- insight::get_parameters(stmod_ak_stally1)

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

prior_summary(stmod_ak_stally1)

## Priors for model 'stmod_ak_stally1' 
## ------
## 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.28,5.28])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

#### PP check plot
ppc.plot(stmod_ak_stally1) +
  labs(caption = "AK: stem_tally~TC")

ggsave("./output/figures_exported/aspen_stem_count/ppc_willowht_ak.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/ppc_willowht_ak.pdf", width = 4.75, height = 3.75)

# library("bayesplot")
color_scheme_set("teal")
mcmc_areas(stmod_ak_stally1,
           prob = 0.8) +
  ggtitle("Posterior distributions",
                      "with medians and 80% credible intervals")

Effect Description

fun.sexit.gt(stmod_ak_stally1)

Effect Existence and Significance Testing
Parameter	Median	CI	CI_low	CI_high	Direction	Significance	Large
(Intercept)	0.57	0.9	0.15	0.96	0.98	0.98	0.87
time_class2013	0.19	0.9	-0.05	0.45	0.90	0.83	0.24
time_class2018	0.38	0.9	0.14	0.61	1.00	0.99	0.71

fun.sexit.gt(stmod_ak_stally1) %>%
  tab_header(title = "Aspen Count Posterior Description", subtitle = "Kawuneeche Valley") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/xpostdesc_stemcnt_TCxF_anc.rtf")

### Model posterior parameters to table
fun.desc.post1.gt.rope(stmod_ak_stally1) %>% 
  tab_header(title = "Aspen Count Posterior Description", subtitle = "Kawuneeche Valley")

## save as RTF
fun.desc.post1.gt.rope(stmod_anc_stally1) %>% 
  tab_header(title = "Aspen Count Posterior Description", subtitle = "Kawuneeche Valley") %>%
  gt::gtsave(filename = "./output/tables/aspen_stem_count/postdesc_stemcnt_TCxF_anc.rtf")

Probability of Direction

# Visualize the pd
#### PD - ANC
p_direction(stmod_ak_stally1) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Aspen stem count, Noncore winter range") +
  scale_fill_manual(values = colfunc3(2))

# save plots
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak.png", dpi = 300, width = 4.75, height = 3.75)
ggsave("./output/figures_exported/aspen_stem_count/pd_stemcnt_ak.pdf", width = 4.75, height = 3.75)

Region of Practical Equivalence (ROPE)

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

rope.stemcnt.ak %>% 
  gt::gtsave(filename = "./output/tables/aspen_stem_count/rope_stemcnt_ak.rtf")

Contrasts

#### Contrasts
# EMMs for later factor levels are subtracted from those for earlier levels; if you want the comparisons to go in the other direction, use pairs.

emmeans(stmod_ak_stally1, ~ 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 = "Noncore Winter Range, Aspen Stem Count")

ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TC_anc.pdf", width = 4.5, height = 3.55) # save plot
ggsave("./output/figures_exported/aspen_stem_count/emm_stemcnt_TC_anc.png", width = 4.5, height = 3.55) # save plot

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(stmod_ak_stally1, ~ 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 (not the response) scale.")

Pairwise contrasts
Results are given on the log (not the response) scale.
contrast	estimate	lower.HPD	upper.HPD
2013 - BL	0.19	-0.10	0.49
2018 - BL	0.38	0.10	0.67
2018 - 2013	0.19	-0.09	0.45

# table of contrasts. Results are given on the response scale.
emmeans(stmod_ak_stally1, ~ 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.21	0.88	1.59
2018 / BL	1.46	1.06	1.90
2018 / 2013	1.20	0.90	1.55

# For 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.

## version for plot combination
# ANC 
emmeans(stmod_ak_stally1, ~ 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 = "Noncore Winter Range, Aspen Stem Count")

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.
shiny	1.6.0	Winston Chang, Joe Cheng, JJ Allaire, Carson Sievert, Barret Schloerke, Yihui Xie, Jeff Allen, Jonathan McPherson, Alan Dipert and Barbara Borges (2021). shiny: Web Application Framework for R. R package version 1.6.0. https://CRAN.R-project.org/package=shiny
shinystan	2.5.0	Jonah Gabry (2018). shinystan: Interactive Visual and Numerical Diagnostics and Posterior Analysis for Bayesian Models. R package version 2.5.0. https://CRAN.R-project.org/package=shinystan
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.”

Kruschke, John K. 2018. “Rejecting or Accepting Parameter Values in Bayesian Estimation.” Advances in Methods and Practices in Psychological Science 1 (2): 270–80. https://doi.org/10.1177/2515245918771304.

Makowski, Dominique, Mattan Ben-Shachar, and Daniel Lüdecke. 2019. “bayestestR: Describing Effects and Their Uncertainty, Existence and Significance Within the Bayesian Framework.” Journal of Open Source Software 4 (40): 1541.

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 B. EVMP Aspen Statistical Analyses

Introduction

Methods

Likelihood

Posterior

Aspen count

Results

Comparison of Sapling Classes

Short Saplings (<=1.5m Tall)) - Combined Winter Range (AC+ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Short Saplings (<=1.5m Tall) - Core Winter Range (AC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Short Saplings (<=1.5m Tall) - Noncore Winter Range (ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Short Saplings (<=1.5m Tall)) - Kawuneeche Valley (AK)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Effects of Burning on Short Saplings (<1.5m Tall) - Combined Winter Range (AC+ANC)

Probability of Direction

Tall Saplings (1.5-2.5m Tall) - Combined Winter Range (AC+ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Tall Saplings (1.5-2.5m Tall) - Core Winter Range (AC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Tall Saplings (1.5-2.5m Tall) - Noncore Winter Range (ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Tall Saplings (1.5-2.5m Tall) - Kawuneeche Valley (AK)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Effects of Burning on Tall Saplings (1.5-2.5m Tall) - Combined Winter Range (AC+ANC)

Probability of Direction

All Stem Diameters - Combined Winter Range (AC+ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

All Stem Diameters - Core Winter Range (AC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

Effects of Burning on All Stems - Core Winter Range (AC)

All Stem Diameters - Noncore Winter Range (ANC)

Modeling and Posterior Description

Effect Description

Probability of Direction

Region of Practical Equivalence (ROPE)

Contrasts

All Stem Diameters - Kawuneeche Valley (AK)