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

Introduction

This document is an appendix to 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” presenting additional analyses of line intercept data not presented in the main report. The code and analyses aim to characterize basic data structure, clean and reorganize as needed for plotting and modeling. Analyses are limited to line intercept data from plots established in areas of the elk winter range and Kawuneeche Valley (i.e., see Zeigenfuss et al. 2011). 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 build on derived data produced in a separate code document focused on ingesting, compiling, and cleaning raw data provided by RMNP staff (Appendix B). Extensive “munging” of raw files was required. Separate code documents were developed for other distinct elements of the larger EVMP (e.g., aspen and upland plots).
In general, packages included in the “Tidyverse” ecosystem of R packages were used for data transformation and visualization (Wickham et al. 2019), although specialized packages particular to specific tasks (e.g., the “bayesplot” plot for visualization of Bayesian posterior distributions) were also used (Gabry et al. 2019).

Bayesian repeated measures analyses were fit separately for combined core and non-core winter range, core winter range (WC) and non-core winter range plots (WNC). For willow height, models were fit using weekly informative priors and a gamma distribution with the “stan_lmer” function in the “rstanarm” package (Goodrich et al. 2020)(Brilleman et al. 2018).

For continuous proportions such as plant cover (i.e, those not derived from count data), a common approach to modeling is to transform data then use ordinary linear models, however this creates problems for interpretation and inference (Douma and Weedon 2019). Techniques relying on transformations make estimates on the transformed scale, requiring back‐transformation for reporting and interpretation. However, the relationship between the original and transformed proportions is often non‐linear, creating issues for interpretation (Douma and Weedon 2019). Alternative approaches using beta regression (Cribari-Neto and Zeileis 2009)(Ferrari and Cribari-Neto 2004) have been developed and are relied on here for modeling continuous proportions (e.g., willow cover). After fitting appropriate models using “rstanarm” functions, comparisons of estimated marginal means were made between time classes and management classes (e.g., fencing) using functions in the “emm” R package (Lenth 2020).

## read in data from the derived folder. See appx_
csv.all.lc.li.df <- read_csv("./data/EVMP_derived/line_intercept_willow_cleaned.csv")
### Counts of plots in fenced and unfenced contexts
## address missing attributes for pType
csv.all.lc.li.df <- csv.all.lc.li.df %>%
  mutate(pType = case_when(is.na(pType) ~ "willow",
                            TRUE ~ pType))

## address missing attributes for RANGE_TYPE
csv.all.lc.li.df <- csv.all.lc.li.df %>%
  mutate(RANGE_TYPE = case_when(is.na(RANGE_TYPE)  & SITE_TYPE == "WC" ~ "core winter range",
                            TRUE ~ RANGE_TYPE))
## Make NA unburned; presumes NA are "unburned"
## standardize encoding of BURNED
csv.all.lc.li.df <- csv.all.lc.li.df %>%
  mutate(BURNED = case_when(is.na(BURNED) ~ "Unburned",
                            BURNED == "Not burned" ~ "Unburned",
                            TRUE ~ BURNED))

## reverse factor levels
csv.all.lc.li.df <- csv.all.lc.li.df %>%
  mutate(timeClass = forcats::fct_rev(timeClass))

## Reorder levels
csv.all.lc.li.df <- csv.all.lc.li.df %>% 
  mutate(SITE_TYPE = as_factor(SITE_TYPE)) %>% 
  mutate(SITE_TYPE = forcats::fct_relevel(SITE_TYPE, "WK", after = 2))

csv.all.lc.li.df <- csv.all.lc.li.df %>% 
  mutate(timeClass = fct_relevel(timeClass, "BL", after=Inf)) %>% 
  mutate(timeClass = fct_rev(timeClass)) %>% 
  mutate(FENCED = fct_rev(FENCED))
## Address missing RANGE_TYPE attribution
csv.all.lc.li.df <- csv.all.lc.li.df %>%
  mutate(RANGE_TYPE = case_when(SITE_TYPE == "WC" ~ "core winter range",
                                TRUE ~ RANGE_TYPE)) %>% 
  distinct()

Results

Maximum Height

Descriptive Statistics and plots

All Shrubs Species

## Boxplot
csv.all.lc.li.df %>% 
  filter(pType == "willow") %>%
  mutate(yr = as.character(yr)) %>%
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill = FENCED), outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  ylim(0,300) +
  labs(x="", y= "Height (cm)", title = "Maximum shrub height", caption = "WC, All shrub species line intercept plots") +
  scale_fill_manual(values = c("grey90","grey50")) +
  theme_minimal() + 
  # theme(axis.text.x=element_text(angle = 45, hjust = 1)) #+
  facet_wrap(~SITE_TYPE)

# ggsave("./output/figures_exported/WCWNCWK_LI_shrubHt_boxplot.png", width = 6.5, height = 4.875)
#### Table of mean max height
csv.all.lc.li.df %>% 
  filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>% 
  group_by(timeClass, SITE_TYPE) %>%
  descr(MAX_HEIGHT_CM, stats = "common") %>% 
  summarytools::tb() %>% 
  mutate(across(c('mean','sd','pct.valid'), ~round(.,digits = 1))) %>% 
  gt() %>% 
  tab_header(title = "Max height", subtitle = "All shrub species combined, all site and mgt group types")
Max height
All shrub species combined, all site and mgt group types
timeClass SITE_TYPE variable mean sd min med max n.valid pct.valid
BL WC MAX_HEIGHT_CM 83.5 81.3 3.0 60 535 154 92.8
BL WNC MAX_HEIGHT_CM 166.6 128.4 25.0 120 520 79 92.9
2013 WC MAX_HEIGHT_CM 97.7 99.5 15.0 60 750 247 95.4
2013 WNC MAX_HEIGHT_CM 222.2 138.8 15.0 210 500 95 96.0
2018 WC MAX_HEIGHT_CM 113.4 92.9 0.7 80 810 409 96.2
2018 WNC MAX_HEIGHT_CM 174.7 129.3 10.0 130 610 163 97.0
2018 WK MAX_HEIGHT_CM 81.7 44.3 15.0 75 250 47 97.9 
## all shrub species
csv.all.lc.li.df %>%
  filter(pType == "willow") %>% 
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  filter(!is.na(MAX_HEIGHT_CM)) %>%
  filter(SITE_TYPE == "WC") %>% 
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill=FENCED), color = 'black', outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  scale_fill_viridis(discrete = TRUE) +
  theme_minimal() +
  ylim(0,320) +
  scale_fill_manual(values = c("grey90","grey50")) +
  labs(x = "Year", y = "Max height (cm)", fill = "", title = "Core Winter Range Plots: All Shrub Species", caption = "Line intercept data")

### All shrubs, site type, fenced/unfenced
#### Table of mean max height
csv.all.lc.li.df %>% 
  filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>% 
  group_by(timeClass, SITE_TYPE, FENCED) %>%
  descr(MAX_HEIGHT_CM, stats = "common") %>% 
  summarytools::tb() %>% 
  mutate(across(c('mean','sd','pct.valid'), ~round(.,digits = 1))) %>% 
  gt() %>% 
  tab_header(title = "LI, Mean max height", subtitle = "All shrubs, all site types, fenced/unfenced")
LI, Mean max height
All shrubs, all site types, fenced/unfenced
timeClass SITE_TYPE FENCED variable mean sd min med max n.valid pct.valid
BL WC Unfenced MAX_HEIGHT_CM 95.0 101.1 3.0 60 535 91 91.9
BL WC Fenced MAX_HEIGHT_CM 66.8 31.8 5.0 65 190 63 94.0
BL WNC Unfenced MAX_HEIGHT_CM 166.6 128.4 25.0 120 520 79 92.9
2013 WC Unfenced MAX_HEIGHT_CM 106.2 115.5 15.0 60 750 158 95.8
2013 WC Fenced MAX_HEIGHT_CM 82.6 59.5 15.0 65 300 89 94.7
2013 WNC Unfenced MAX_HEIGHT_CM 222.2 138.8 15.0 210 500 95 96.0
2018 WC Unfenced MAX_HEIGHT_CM 106.0 110.8 0.7 60 810 210 95.0
2018 WC Fenced MAX_HEIGHT_CM 121.1 68.5 10.0 120 290 199 97.5
2018 WNC Unfenced MAX_HEIGHT_CM 174.7 129.3 10.0 130 610 163 97.0
2018 WK Unfenced MAX_HEIGHT_CM 56.2 44.8 15.0 50 250 24 96.0
2018 WK Fenced MAX_HEIGHT_CM 108.3 24.0 65.0 110 155 23 100.0 
### All shrubs, site type, fenced/unfenced, burned/unburned
## Boxplot
csv.all.lc.li.df %>% 
  filter(pType == "willow") %>%
  mutate(yr = as.character(yr)) %>%
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill = FENCED), outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  ylim(0,300) +
  labs(x="", y= "Height (cm)", title = "Maximum shrub height", caption = "WC, All shrub species line intercept plots") +
  scale_fill_manual(values = c("grey90","grey50")) +
  theme_minimal() + 
  # theme(axis.text.x=element_text(angle = 45, hjust = 1)) #+
  facet_wrap(BURNED~SITE_TYPE)

# ggsave("./output/figures_exported/WCWNCWK_LI_shrubHt_FenBur_boxplot.png", width = 6.5, height = 4.875)

Willow Species

csv.wil.lc.li.df <- csv.all.lc.li.df %>% 
  filter(str_detect(SPECIES_CODE, "^SA")) %>% 
  mutate(FENCED = fct_rev(FENCED))

## Boxplot
csv.wil.lc.li.df %>% 
  filter(pType == "willow") %>%
  mutate(yr = as.character(yr)) %>%
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill = FENCED), outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  ylim(0,600) +
  labs(x="", y= "Height (cm)", title = "Maximum willow height", caption = "WC, All willow species line intercept plots") +
  scale_fill_manual(values = c("grey90","grey50")) +
  theme_minimal() + 
  facet_wrap(~SITE_TYPE)

### All willow species, site type, fenced/unfenced
#### Table of mean max height
csv.wil.lc.li.df %>% 
  filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>% 
  group_by(timeClass, SITE_TYPE, FENCED) %>%
  descr(MAX_HEIGHT_CM, stats = "common") %>% 
  summarytools::tb() %>% 
  mutate(across(c('mean','sd','pct.valid'), ~round(.,digits = 1))) %>% 
  gt() %>% 
  tab_header(title = "Maximum height", subtitle = "All willows combined, all site types, fenced/unfenced")
Maximum height
All willows combined, all site types, fenced/unfenced
timeClass SITE_TYPE FENCED variable mean sd min med max n.valid pct.valid
BL WC Fenced MAX_HEIGHT_CM 72.7 36.0 10 70.0 190 38 97.4
BL WC Unfenced MAX_HEIGHT_CM 97.3 98.3 3 60.0 440 66 97.1
BL WNC Unfenced MAX_HEIGHT_CM 170.0 127.4 25 120.0 520 74 100.0
2013 WC Fenced MAX_HEIGHT_CM 109.2 62.4 25 97.5 300 52 100.0
2013 WC Unfenced MAX_HEIGHT_CM 129.3 109.3 20 80.0 470 72 100.0
2013 WNC Unfenced MAX_HEIGHT_CM 221.2 139.4 15 215.0 500 82 100.0
2018 WC Fenced MAX_HEIGHT_CM 151.3 59.2 30 145.0 290 135 99.3
2018 WC Unfenced MAX_HEIGHT_CM 130.5 102.4 20 95.0 450 99 100.0
2018 WNC Unfenced MAX_HEIGHT_CM 175.3 119.9 10 135.0 610 139 100.0
2018 WK Fenced MAX_HEIGHT_CM 108.3 24.0 65 110.0 155 23 100.0
2018 WK Unfenced MAX_HEIGHT_CM 50.8 15.5 30 52.5 75 18 100.0 
### All Willow Species, site type, fenced/unfenced, burned/unburned
## Boxplot
csv.wil.lc.li.df %>% 
  filter(pType == "willow") %>%
  mutate(yr = as.character(yr)) %>%
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill = FENCED), outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  ylim(0,600) +
  labs(x="", y= "Height (cm)", title = "Maximum shrub height", caption = "WC, All willow species line intercept plots") +
  scale_fill_manual(values = c("grey90","grey50")) +
  theme_minimal() + 
  facet_wrap(BURNED~SITE_TYPE)

#### Table of mean max height
csv.wil.lc.li.df %>% 
  filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>% 
  group_by(timeClass, SITE_TYPE, FENCED, BURNED) %>%
  descr(MAX_HEIGHT_CM, stats = "common") %>% 
  summarytools::tb() %>% 
  mutate(across(c('mean','sd','pct.valid'), ~round(.,digits = 1))) %>% 
  gt() %>% 
  tab_header(title = "Maximum height", subtitle = "All willow, all site types, fenced/unfenced, burned/unburned")
Maximum height
All willow, all site types, fenced/unfenced, burned/unburned
timeClass SITE_TYPE FENCED BURNED variable mean sd min med max n.valid pct.valid
BL WC Fenced Burned MAX_HEIGHT_CM 46.9 31.6 10 45.0 110 8 100.0
BL WC Fenced Unburned MAX_HEIGHT_CM 79.6 34.4 30 70.0 190 30 96.8
BL WC Unfenced Burned MAX_HEIGHT_CM 108.2 95.1 40 76.0 275 5 100.0
BL WC Unfenced Unburned MAX_HEIGHT_CM 96.4 99.3 3 50.0 440 61 96.8
BL WNC Unfenced Unburned MAX_HEIGHT_CM 170.0 127.4 25 120.0 520 74 100.0
2013 WC Fenced Burned MAX_HEIGHT_CM 66.7 28.1 25 60.0 110 18 100.0
2013 WC Fenced Unburned MAX_HEIGHT_CM 131.8 64.1 40 112.5 300 34 100.0
2013 WC Unfenced Burned MAX_HEIGHT_CM 45.0 21.2 30 45.0 60 2 100.0
2013 WC Unfenced Unburned MAX_HEIGHT_CM 131.7 109.8 20 80.0 470 70 100.0
2013 WNC Unfenced Unburned MAX_HEIGHT_CM 221.2 139.4 15 215.0 500 82 100.0
2018 WC Fenced Burned MAX_HEIGHT_CM 134.6 37.8 50 142.5 225 66 100.0
2018 WC Fenced Unburned MAX_HEIGHT_CM 167.2 70.8 30 175.0 290 69 98.6
2018 WC Unfenced Burned MAX_HEIGHT_CM 32.5 6.5 25 32.5 40 4 100.0
2018 WC Unfenced Unburned MAX_HEIGHT_CM 134.6 102.5 20 95.0 450 95 100.0
2018 WNC Unfenced Unburned MAX_HEIGHT_CM 175.3 119.9 10 135.0 610 139 100.0
2018 WK Fenced Unburned MAX_HEIGHT_CM 108.3 24.0 65 110.0 155 23 100.0
2018 WK Unfenced Unburned MAX_HEIGHT_CM 50.8 15.5 30 52.5 75 18 100.0 
csv.all.lc.li.df %>%
  filter(yr !=2009 & yr != 2015 & yr !=2017) %>% 
  filter(SITE_TYPE != "WK") %>%
  filter(SPECIES_CODE == "SAMO" | SPECIES_CODE == "SAGE" | SPECIES_CODE == "SAPL" ) %>% 
  filter(!is.na(MAX_HEIGHT_CM)) %>% 
  group_by(SPECIES_CODE, timeClass, FENCED) %>% 
  descr(MAX_HEIGHT_CM, stats = "common") %>% 
  tb() %>% 
  select(-variable) %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  gt() %>% 
  tab_header(title = "Maximum height: common willow species comparisons", subtitle = "Combined WC and WNC; SAPL, SAGE, and SAMO")
Maximum height: common willow species comparisons
Combined WC and WNC; SAPL, SAGE, and SAMO
SPECIES_CODE timeClass FENCED mean sd min med max n.valid pct.valid
SAGE BL Unfenced 210.6 152.1 3 245.0 455 21 100
SAGE 2013 Unfenced 227.2 148.3 30 250.0 490 27 100
SAGE 2013 Fenced 123.0 35.1 85 115.0 180 5 100
SAGE 2018 Unfenced 157.6 129.0 30 85.0 410 25 100
SAGE 2018 Fenced 156.0 65.6 60 155.0 270 15 100
SAMO BL Unfenced 115.1 111.0 10 70.0 440 83 100
SAMO BL Fenced 71.1 36.8 10 70.0 190 35 100
SAMO 2013 Unfenced 159.3 127.9 15 100.0 500 77 100
SAMO 2013 Fenced 113.3 74.1 30 95.0 300 23 100
SAMO 2018 Unfenced 161.6 115.9 10 130.0 610 132 100
SAMO 2018 Fenced 155.7 71.4 40 147.5 280 42 100
SAPL BL Unfenced 164.7 117.6 25 150.0 520 19 100
SAPL BL Fenced 82.5 17.7 70 82.5 95 2 100
SAPL 2013 Unfenced 154.0 121.5 40 102.5 450 34 100
SAPL 2013 Fenced 108.8 66.8 25 97.5 270 16 100
SAPL 2018 Unfenced 124.0 95.3 20 95.0 380 36 100
SAPL 2018 Fenced 147.3 68.0 30 140.0 290 28 100 
## salix monitcola
csv.all.lc.li.df %>%
  filter(pType == "willow") %>% 
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  filter(!is.na(MAX_HEIGHT_CM)) %>%
  filter(SPECIES_CODE =="SAMO") %>% 
  filter(SITE_TYPE == "WC") %>% 
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill=FENCED), color = 'black', outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  theme_minimal() +
  ylim(0,320) +
  scale_fill_manual(values = c("grey90","grey50")) +
  labs(x = "Year", y = "Max height (cm)", fill = "", title = "Core Winter Range Plots: SAMO only", caption = "Line intercept data")

## salix geyeriana
csv.all.lc.li.df %>%
  filter(pType == "willow") %>% 
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  filter(!is.na(MAX_HEIGHT_CM)) %>%
  filter(SPECIES_CODE =="SAGE") %>% 
  filter(SITE_TYPE == "WC") %>% 
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill=FENCED), color = 'black', outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  theme_minimal() +
  ylim(0,320) +
  scale_fill_manual(values = c("grey90","grey50")) +
  labs(x = "Year", y = "Max height (cm)", fill = "", title = "Core Winter Range Plots: SAGE only", caption = "Line intercept data")

## Salix planifolia
csv.all.lc.li.df %>%
  filter(pType == "willow") %>% 
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  filter(!is.na(MAX_HEIGHT_CM)) %>%
  filter(SPECIES_CODE =="SAPL") %>% 
  filter(SITE_TYPE == "WC") %>% 
  ggplot(aes(timeClass, MAX_HEIGHT_CM)) +
  geom_boxplot(aes(fill=FENCED), color = 'black', outlier.shape = NA) +
  geom_hline(aes(yintercept = 110), color = "red", lty = "dashed", size = 1) +
  theme_minimal() +
  ylim(0,320) +
  scale_fill_manual(values = c("grey90","grey50")) +
  labs(x = "Year", y = "Max height (cm)", fill = "", title = "Core Winter Range Plots: SAPL only", caption = "Line intercept data")

csv.all.lc.li.df %>%
  filter(pType == "willow") %>% 
  filter(SPECIES_CODE != "SAXX") %>% 
  filter(yr == 2008 | yr == 2013 | yr == 2018) %>%
  filter(!is.na(MAX_HEIGHT_CM)) %>%
  filter(FENCED == "Unfenced" | FENCED == "Fenced") %>% 
  group_by(yr, SPECIES_CODE, FENCED) %>% 
  summarise(mean.ht = mean(MAX_HEIGHT_CM, na.rm = TRUE)) %>% 
  pivot_wider(names_from = yr, values_from = c(mean.ht)) %>% 
  gt::gt() %>% 
  fmt_number(
    columns = vars('2008','2013','2018'),
    decimals = 1,
    use_seps = FALSE
  ) %>% 
  tab_header("Mean height by shrub species")
Mean height by shrub species
FENCED 2008 2013 2018
BEGL
Fenced 41.7 130.0 165.0
Unfenced NA 126.7 108.8
BEOC
Unfenced 322.5 348.0 176.9
Fenced NA 45.0 82.5
DAFR
Unfenced 52.1 49.3 51.8
Fenced 60.0 42.9 52.4
SAGE
Unfenced 210.6 227.2 153.7
Fenced NA 123.0 141.5
SALUC
Unfenced 80.0 NA NA
SAMO
Unfenced 115.1 159.3 157.8
Fenced 71.1 113.3 149.9
SAPL
Unfenced 164.7 154.0 107.2
Fenced 82.5 108.8 144.8
ALIN
Unfenced NA 330.0 345.3
Fenced NA NA 170.0
LOIN
Unfenced NA 112.5 122.5
Fenced NA 40.0 93.3
POTR
Unfenced NA 126.9 NA
RIIN
Unfenced NA 78.6 59.1
ROWO
Unfenced NA 45.9 43.6
Fenced NA NA 33.2
RUDE
Unfenced NA 50.0 NA
SABE
Unfenced NA 239.6 137.8
Fenced NA NA 131.1
SADR
Unfenced NA 231.7 253.1
Fenced NA NA 114.2
SAPE
Fenced NA 86.2 147.3
Unfenced NA NA 63.3
JUCO
Unfenced NA NA 40.0
RUID
Unfenced NA NA 35.0
SABR
Unfenced NA NA 40.0
SAER
Fenced NA NA 156.0

Modeling and Posterior Description - Combined Willow Species

# clean names
csv.wil.lc.li.df <- csv.wil.lc.li.df %>% 
  clean_names() %>% 
  filter(!is.na(max_height_cm)) %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018")

## run with gamma
wc.wnc.stmod_li_TCxF <- stan_glmer(max_height_cm ~ time_class*fenced +  (1 | site_id), data = csv.wil.lc.li.df,
                      family=Gamma(link="log"),
                      iter = 10000,
                      )
## Priors for model 'wc.wnc.stmod_li_TCxF' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [5.71,5.00,5.40,...])
## 
## Auxiliary (shape)
##  ~ exponential(rate = 1)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details
#### PP check plot
ppc.plot(wc.wnc.stmod_li_TCxF) +
  labs(caption = "WC+WNC: willow_ht~TCxF")

#### Describe posterior distributions
fun.desc.post1.gt.rope(wc.wnc.stmod_li_TCxF) %>% 
  tab_header(title = "Combined WC and WNC", subtitle = "Willow Height Posterior Description - Line Intercept")
Combined WC and WNC
Willow Height Posterior Description - Line Intercept
Parameter Median MAD CI CI_low CI_high pd ps ROPE_CI ROPE_low ROPE_high ROPE_Percentage Rhat ESS
(Intercept) 4.12 0.14 0.90 3.90 4.35 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 2.75K
time_class2013 0.46 0.10 0.90 0.29 0.63 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 7.83K
time_class2018 0.82 0.09 0.90 0.66 0.97 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 7.63K
fencedUnfenced 0.56 0.16 0.90 0.31 0.84 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 2.60K
time_class2013:fencedUnfenced −0.19 0.12 0.90 −0.38 0.00 0.95 0.78 90.00 −0.10 0.10 0.19 1.00 7.75K
time_class2018:fencedUnfenced −0.60 0.11 0.90 −0.78 −0.43 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 7.51K
shape 4.87 0.25 0.90 4.44 5.27 1.00 1.00 90.00 −0.10 0.10 0.00 1.00 15.91K

Probability of Direction

p_direction(wc.wnc.stmod_li_TCxF) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(caption = "Willow height, Combined Core and Non-core Winter Range")

p_direction(wc.wnc.stmod_li_TCxF) %>% 
  as_tibble() %>% 
  gt() %>% 
  tab_header(title = "Probability of Direction", subtitle = "WC & WNC Combined Willows")
Probability of Direction
WC & WNC Combined Willows
Parameter pd Effects Component
(Intercept) 1.0000 fixed conditional
time_class2013 1.0000 fixed conditional
time_class2018 1.0000 fixed conditional
fencedUnfenced 0.9995 fixed conditional
time_class2013:fencedUnfenced 0.9508 fixed conditional
time_class2018:fencedUnfenced 1.0000 fixed conditional
shape 1.0000 fixed conditional

Contrasts

## contrasts
emmeans(wc.wnc.stmod_li_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  # plot() + # on the lm predict scale
  plot(type = "response") +
  theme_minimal() +
  geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
  labs(subtitle = "Point estimate displayed: median 
HPD interval probability: 0.95", x = "Estimate", y = "Contrast", caption = "WC & WNC Max Height, Combined Willow, Results are given on the response scale.")

# table of contrasts. Results are given on the response scale.
emmeans(wc.wnc.stmod_li_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE, type="response") %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the response scale.")
Pairwise contrasts
Results are given on the response scale.
contrast fenced ratio lower.HPD upper.HPD
2013 / BL Fenced 1.59 1.28 1.92
2018 / BL Fenced 2.26 1.88 2.69
2018 / 2013 Fenced 1.43 1.21 1.65
2013 / BL Unfenced 1.31 1.16 1.45
2018 / BL Unfenced 1.24 1.11 1.37
2018 / 2013 Unfenced 0.95 0.86 1.04 
## contrasts: 
emmeans(wc.wnc.stmod_li_TCxF, ~ time_class | 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 = "WC & WNC Max Height, Combined Willow, Results are given on the log (not response) scale.")

# table of contrasts. Results are given on the log (not the response) scale.
emmeans(wc.wnc.stmod_li_TCxF, ~ time_class | fenced) %>% 
  pairs(reverse = TRUE) %>% 
  as_tibble() %>% 
  mutate(across(where(is.numeric),round, 2)) %>% 
  gt() %>% 
  tab_header(title = "Pairwise contrasts", subtitle = "Results are given on the log (not the response) scale.")
Pairwise contrasts
Results are given on the log (not the response) scale.
contrast fenced estimate lower.HPD upper.HPD
2013 - BL Fenced 0.46 0.25 0.66
2018 - BL Fenced 0.82 0.63 0.99
2018 - 2013 Fenced 0.36 0.21 0.51
2013 - BL Unfenced 0.27 0.16 0.38
2018 - BL Unfenced 0.21 0.11 0.31
2018 - 2013 Unfenced -0.06 -0.15 0.04

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 March 10

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/.
Cribari-Neto, Francisco, and Achim Zeileis. 2009. “Beta Regression in r.”
Douma, Jacob C., and James T. Weedon. 2019. “Analysing Continuous Proportions in Ecology and Evolution: A Practical Introduction to Beta and Dirichlet Regression.” Methods in Ecology and Evolution 10 (9): 1412–30. https://doi.org/10.1111/2041-210X.13234.
Ferrari, Silvia, and Francisco Cribari-Neto. 2004. “Beta Regression for Modelling Rates and Proportions.” Journal of Applied Statistics 31 (7): 799–815.
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.”
Lenth, R. 2020. “Estimated Marginal Means, Aka Least-Squares Means. R Package Version 1.4.6. Https://CRAN.r-Project.org/Package=emmeans.”
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.
———. 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.
———. 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.
———. 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.”