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

Introduction

This document contains code and results from frequentist statistical analyses of Elk Vegetation Management Plan (EVMP) data collected through the 2018 sampling season for a subset of data as requested by RMNP staff. Analyses are provided as a supplement to the NPS NRR Report “Monitoring of Vegetation Response to Elk Population and Habitat Management in Rocky Mountain National Park - Analysis of Elk Vegetation Management Plan monitoring Data: 2008–2018.” The frequentist analyses here differ in philosophical foundation and interpretation than the Bayesian analyses presented in the report and are presented for comparison to past analyses.

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. See other Appendices for details on data processing used to transform raw files to the derived data ingested for analyses here.

Per guidance from RMNP, for willow plots, only macroplot data are analyzed here. Based on data structure, a repeated measures Analysis of Variance (ANOVA) approach is used. Following methods used in Zeigenfuss and Johnson (2015), the minimum non-zero proportion of shrub cover was substituted for areas with no shrub cover and (1 - this minimum) was used for sites with a proportion of 1 so that a log of proportion of cover coul be tested statistically (Zeigenfuss and Johnson 2015).

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

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

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

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

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

### calculate the mean and max height for all the willows
willow.ht <- willow.ht %>% 
  group_by(time_class,site_type, fenced, burned, site_id) %>% 
  summarise(ht_mean = mean(plant_ht_cm, na.rm=TRUE), ht_max = max(plant_ht_cm, na.rm=TRUE)) %>% 
  mutate(burned = case_when(is.na(burned) ~ "N",
                            TRUE ~ burned)) %>% 
  ungroup()


# parse out the WK plots
willow.ht.wk <- willow.ht %>%
  filter(site_type == "WK")


## 
## filter to 5 yr sampling interval, chr to factor conversions
willow.ht <- willow.ht %>% 
  filter(time_class == "BL" | time_class == "2013" | time_class == "2018") %>%
  mutate(time_class = as.factor(time_class)) %>% 
  mutate(time_class = fct_relevel(time_class, "2018", "2013", "BL")) %>% 
  mutate(fenced = as.factor(fenced))

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

## factor reverse: fenced
willow.ht <- willow.ht %>% 
  mutate(fenced = fct_rev(fenced))
## create separate df for WC, WNC, combined WC & WNC, and KV
willow.ht.wc <- willow.ht %>% 
  filter(site_type == "WC") 

# non core
willow.ht.nc <- willow.ht %>% 
  # distinct(site_type)
  filter(site_type == "WNC") 

## combined winter range (WC and WNC)
willow.ht.wcwnc <- willow.ht %>% 
  filter(site_type == "WNC" | site_type == "WC")

# WK defined above

Results

Willow Height

Combined core and Non-core range

# time class only
aov_model_wcwnc <- aov(ht_mean ~ time_class + Error(1/site_id), data = willow.ht.wcwnc)
summary(aov_model_wcwnc)
## 
## Error: Within
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## time_class   2   49318   24659    2.76 0.0651 .
## Residuals  275 2457341    8936                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# time x fenced x location
aov_model_wcwnc.tcXfXsite_type <- aov(ht_mean ~ time_class*fenced*site_type + Error(1/site_id), data = willow.ht.wcwnc)

broom::tidy(aov_model_wcwnc.tcXfXsite_type) %>% 
  mutate(across(where(is.numeric),round,2)) %>% 
  gt() %>% 
  tab_header(title = "ht_mean ~ time_class*fenced*site_type + Error(1/site_id)")
ht_mean ~ time_class*fenced*site_type + Error(1/site_id)
stratum term df sumsq meansq statistic p.value
Within time_class 2 49318.36 24659.18 3.36 0.04
Within fenced 1 56195.44 56195.44 7.65 0.01
Within site_type 1 373791.01 373791.01 50.89 0.00
Within time_class:fenced 2 50881.55 25440.77 3.46 0.03
Within time_class:site_type 2 553.81 276.90 0.04 0.96
Within Residuals 269 1975919.69 7345.43 NA NA 
# add burned
aov_model_wcwnc.tcXfXbXsite_type <- aov(ht_mean ~ time_class*site_type*fenced*burned + Error(1/site_id), data = willow.ht.wcwnc)

broom::tidy(aov_model_wcwnc.tcXfXbXsite_type) %>% 
  mutate(across(where(is.numeric),round,2)) %>% 
  gt() %>% 
  tab_header(title = "plant_ht_cm ~ time_class*site_type*fenced*burned + Error(1/site_id)")
plant_ht_cm ~ time_class*site_type*fenced*burned + Error(1/site_id)
stratum term df sumsq meansq statistic p.value
Within time_class 2 49318.36 24659.18 3.54 0.03
Within site_type 1 424969.12 424969.12 61.05 0.00
Within fenced 1 5017.34 5017.34 0.72 0.40
Within burned 1 117511.68 117511.68 16.88 0.00
Within time_class:site_type 2 7154.26 3577.13 0.51 0.60
Within time_class:fenced 2 47366.86 23683.43 3.40 0.03
Within time_class:burned 2 24103.04 12051.52 1.73 0.18
Within fenced:burned 1 17.94 17.94 0.00 0.96
Within time_class:fenced:burned 2 442.06 221.03 0.03 0.97
Within Residuals 263 1830759.19 6961.06 NA NA
ht_mean ~ time_class
Effect DFn DFd F p p<.05 ges
time_class 2 275 2.76 0.065 0.02
ht_mean ~ time_class+fenced+burned
Effect DFn DFd F p p<.05 ges
time_class 2 273 3.674 2.70e-02 * 0.026
fenced 1 273 0.428 5.13e-01 0.002
burned 1 273 30.178 9.03e-08 * 0.100
ht_mean ~ time_class*fenced*burned
Effect DFn DFd F p p<.05 ges
time_class 2 266 3.76 0.03 * 0.03
fenced 1 266 0.49 0.49 0.00
burned 1 266 31.19 0.00 * 0.10
time_class:fenced 2 266 3.82 0.02 * 0.03
time_class:burned 2 266 1.32 0.27 0.01
fenced:burned 1 266 2.24 0.14 0.01
time_class:fenced:burned 2 266 0.00 1.00 0.00

Model summary, combined core and non-core winter range:
The repeated-measures ANOVA (formula: ht_mean ~ time_class * site_type * fenced * burned + Error(1/site_id)) suggests that:

  • The main effect of time_class is significant and small (F(2, 263) = 3.54, p = 0.030; Eta2 (partial) = 0.03, 90% CI [1.34e-03, 0.06])
  • The main effect of site_type is significant and large (F(1, 263) = 61.05, p < .001; Eta2 (partial) = 0.19, 90% CI [0.12, 0.26])
  • The main effect of fenced is not significant and very small (F(1, 263) = 0.72, p = 0.397; Eta2 (partial) = 2.73e-03, 90% CI [0.00, 0.02])
  • The main effect of burned is significant and medium (F(1, 263) = 16.88, p < .001; Eta2 (partial) = 0.06, 90% CI [0.02, 0.11])
  • The interaction between time_class and site_type is not significant and very small (F(2, 263) = 0.51, p = 0.599; Eta2 (partial) = 3.89e-03, 90% CI [0.00, 0.02])
  • The interaction between time_class and fenced is significant and small (F(2, 263) = 3.40, p = 0.035; Eta2 (partial) = 0.03, 90% CI [9.64e-04, 0.06])
  • The interaction between time_class and burned is not significant and small (F(2, 263) = 1.73, p = 0.179; Eta2 (partial) = 0.01, 90% CI [0.00, 0.04])
  • The interaction between fenced and burned is not significant and very small (F(1, 263) = 2.58e-03, p = 0.960; Eta2 (partial) = 9.80e-06, 90% CI [0.00, 0.00])
  • The interaction between time_class, fenced and burned is not significant and very small (F(2, 263) = 0.03, p = 0.969; Eta2 (partial) = 2.41e-04, 90% CI [0.00, 0.00])

Effect sizes were labelled following Field’s (2013) recommendations.

# GT table: location: WC and WNC only
aov_model_wcwnc %>% 
  broom::tidy() %>% 
  mutate(across(where(is.numeric), round,2)) %>% 
  gt() %>% 
  tab_header(title="plant_ht_cm ~ time_class*location*fenced*burned + Error(1/site_id)", subtitle = "WC and WNC")
plant_ht_cm ~ time_class*location*fenced*burned + Error(1/site_id)
WC and WNC
stratum term df sumsq meansq statistic p.value
Within time_class 2 49318.36 24659.18 2.76 0.07
Within Residuals 275 2457341.50 8935.79 NA NA

Core range

stratum term df sumsq meansq statistic p.value
Within time_class 2 65553.866 32776.933 6.4713998 0.001941084
Within fenced 1 4928.339 4928.339 0.9730395 0.325277455
Within time_class:fenced 2 43983.834 21991.917 4.3420319 0.014431970
Within Residuals 176 891420.775 5064.891 NA NA

Model summary, mean plot willow height, core winter range:
The repeated-measures ANOVA (formula: ht_mean ~ time_class * fenced + Error(1/site_id)) suggests that:

  • The main effect of time_class is significant and medium (F(2, 176) = 6.47, p = 0.002; Eta2 (partial) = 0.07, 90% CI [0.02, 0.13])
  • The main effect of fenced is not significant and very small (F(1, 176) = 0.97, p = 0.325; Eta2 (partial) = 5.50e-03, 90% CI [0.00, 0.04])
  • The interaction between time_class and fenced is significant and small (F(2, 176) = 4.34, p = 0.014; Eta2 (partial) = 0.05, 90% CI [5.39e-03, 0.10])

Effect sizes were labelled following Field’s (2013) recommendations.

broom::tidy(aov_model_wc) %>% 
  mutate(across(.cols = where(is.numeric),.fns = round, 2)) %>% 
  gt()
stratum term df sumsq meansq statistic p.value
Within time_class 2 65553.87 32776.93 6.47 0.00
Within fenced 1 4928.34 4928.34 0.97 0.33
Within time_class:fenced 2 43983.83 21991.92 4.34 0.01
Within Residuals 176 891420.77 5064.89 NA NA 
## save as table (.csv)
broom::tidy(aov_model_wc) %>% 
  mutate(across(.cols = where(is.numeric),.fns = round, 3)) %>% 
  write_csv("output/tables/aov_wc_ht_TC_fenced.csv")
## aov_model_wc <- aov(plant_ht_cm ~ time_class*fenced + Error(1/site_id), data = willow.ht.wc)
willow.ht.stype.nest <- willow.ht %>% 
  nest_by(site_type)

### create AOV models

# plant_ht_cm ~ time_class + Error(1/site_id)
# not enough factor levels for WK
mods.aov.willHt.tc <- willow.ht.stype.nest %>%
  # filter(site_type == "WC") %>% 
  filter(site_type != "WK") %>%
  mutate(mod = list(aov(ht_mean ~ time_class + Error(1/site_id), data = data)))

# plant_ht_cm ~ time_class*fenced + Error(1/site_id)
# not enough factor levels for WK OR WNC
mods.aov.willHt.tc.fen <- willow.ht.stype.nest %>%
  filter(site_type == "WC") %>% 
  mutate(mod = list(aov(ht_mean ~ time_class*fenced + Error(1/site_id), data = data))) %>% 
  mutate(mod.max = list(aov(ht_max ~ time_class*fenced + Error(1/site_id), data = data)))


## WNC
mods.aov.willHt.wnc <- willow.ht.stype.nest %>%
  filter(site_type == "WNC") %>% 
  mutate(mod = list(aov(ht_mean ~ time_class + Error(1/site_id), data = data)))%>% 
  mutate(mod.max = list(aov(ht_max ~ time_class + Error(1/site_id), data = data)))


# GT table: location: WC only
aov.willHt.tc.fen <- mods.aov.willHt.tc.fen %>% 
  mutate(aov_summary = list(summary(mod))) %>% 
  mutate(aov_tidy = list(broom::tidy(mod))) %>% 
  unnest(aov_tidy) %>%
  mutate(across(where(is.numeric), round, 3)) %>% 
  select(site_type,term,df,sumsq,meansq,statistic,p.value) %>% 
  gt::gt() %>% 
  gt::tab_header(title = "aov: plant_ht_cm ~ time_class*fenced + Error(1/site_id)", subtitle = "WC")

aov.willHt.tc.fen
aov: plant_ht_cm ~ time_class*fenced + Error(1/site_id)
WC
term df sumsq meansq statistic p.value
WC
time_class 2 65553.866 32776.933 6.471 0.002
fenced 1 4928.339 4928.339 0.973 0.325
time_class:fenced 2 43983.834 21991.917 4.342 0.014
Residuals 176 891420.775 5064.891 NA NA

Non-core range

aov_model_wnc <- aov(ht_mean ~ time_class + Error(1/site_id), data = willow.ht.nc)

Model summary, non-core winter range:
The repeated-measures ANOVA (formula: ht_mean ~ time_class + Error(1/site_id)) suggests that:

  • The main effect of time_class is not significant and small (F(2, 93) = 0.63, p = 0.533; Eta2 (partial) = 0.01, 90% CI [0.00, 0.06])

Effect sizes were labelled following Field’s (2013) recommendations.

# GT table: location: WNC only
aov.willHt.tc.fen.wnc <- mods.aov.willHt.wnc %>% 
  mutate(aov_summary = list(summary(mod))) %>% 
  mutate(aov_tidy = list(broom::tidy(mod))) %>% 
  unnest(aov_tidy) %>%
  mutate(across(where(is.numeric), round, 3)) %>% 
  select(site_type,term,df,sumsq,meansq,statistic,p.value) %>% 
  gt::gt() %>% 
  gt::tab_header(title = "aov: plant_ht_cm ~ time_class + Error(1/site_id)", subtitle = "WNC. All unfenced")

aov.willHt.tc.fen.wnc
aov: plant_ht_cm ~ time_class + Error(1/site_id)
WNC. All unfenced
term df sumsq meansq statistic p.value
WNC
time_class 2 14754.75 7377.374 0.633 0.533
Residuals 93 1084498.91 11661.279 NA NA

Kawuneeche Valley

2016-2018 only

# model for wk
# willow.ht.wk %>% 
#   distinct(time_class)
aov_model_wk <- aov(ht_mean ~ time_class + Error(1/site_id), data = willow.ht.wk)

Model summary, non-core winter range:
The repeated-measures ANOVA (formula: ht_mean ~ time_class + Error(1/site_id)) suggests that:

  • The main effect of time_class is not significant and very small (F(2, 20) = 0.01, p = 0.986; Eta2 (partial) = 1.37e-03, 90% CI [0.00, 0.00])

Effect sizes were labelled following Field’s (2013) recommendations.

aov_model_wk %>% 
  broom::tidy() %>% 
  mutate(across(where(is.numeric), round,2)) %>% 
  gt() %>% 
  tab_header(title="ht_mean ~ time_class + Error(1/site_id)", subtitle = "Kawuneeche Valley")
ht_mean ~ time_class + Error(1/site_id)
Kawuneeche Valley
stratum term df sumsq meansq statistic p.value
Within time_class 2 22.17 11.09 0.01 0.99
Within Residuals 20 16118.89 805.94 NA NA

Mean versus maximum willow height

willow.ht %>%
  filter(site_type != "WK") %>% 
  ggplot(aes(ht_mean, ht_max)) +
  geom_point(aes(color=fenced, shape = fenced)) +
  geom_smooth(aes(color = fenced), se = FALSE, ) +
  geom_abline(intercept = 0, slope = 1, lty = "dashed", size=1, color="gray80") +
  scale_color_manual(values = c("blue","red")) +
  theme_minimal() +
  labs(x="Mean willow height (cm)", y = "Mean willow height (cm)", shape = "", color = "") +
  facet_wrap(~site_type)

Willow Cover

### Combined Core and Noncore Winter Range  

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

canopy.all.willow.sum <- read_csv("./output/exported_data/willow_macroplot_cover_BLto2018.csv") %>% 
  janitor::clean_names()
### create subsets of data for modeling
# pull out WK, note different years of data
cov.wk <- canopy.all.willow.sum %>% 
  filter(site_type == "WK") %>%
  mutate(time_class = as.character(time_class)) %>% 
  filter(time_class == "2016" | time_class == "2017" | time_class == "2018") %>% 
  mutate(time_class = as.factor(time_class))

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

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

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

Summary table: time class and fenced

#### Summary Tables
#### summary tables: time class and fenced
cov.wc.wnc %>% 
  group_by(time_class, fenced) %>% 
  summarytools::descr(cover_allwillow,stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "Combined WC and WNC macroplot data")
Combined WC and WNC macroplot data
time_class fenced variable mean sd min med max n.valid pct.valid
BL Unfenced cover_allwillow 0.33 0.27 0.00 0.27 1 53 100
BL Fenced cover_allwillow 0.16 0.23 0.00 0.10 1 19 100
2013 Unfenced cover_allwillow 0.29 0.29 0.00 0.24 1 64 100
2013 Fenced cover_allwillow 0.28 0.29 0.00 0.21 1 30 100
2018 Unfenced cover_allwillow 0.43 0.35 0.00 0.39 1 69 100
2018 Fenced cover_allwillow 0.54 0.34 0.05 0.50 1 30 100

Summary table: time class,fenced,burned

## time class, fenced and burned
cov.wc.wnc %>% 
  group_by(time_class, fenced, burned) %>% 
  summarytools::descr(cover_allwillow, stats = "common") %>% 
  summarytools::tb() %>%
  mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>% 
  gt() %>% 
  tab_header(title = "Combined WC and WNC macroplot data")
Combined WC and WNC macroplot data
time_class fenced burned variable mean sd min med max n.valid pct.valid
BL Unfenced Unburned cover_allwillow 0.33 0.27 0.00 0.27 1.00 53 100
BL Fenced Unburned cover_allwillow 0.16 0.23 0.00 0.10 1.00 19 100
2013 Unfenced Unburned cover_allwillow 0.33 0.29 0.00 0.26 1.00 54 100
2013 Unfenced Burned cover_allwillow 0.10 0.22 0.01 0.02 0.72 10 100
2013 Fenced Unburned cover_allwillow 0.41 0.29 0.02 0.40 1.00 17 100
2013 Fenced Burned cover_allwillow 0.12 0.21 0.00 0.04 0.79 13 100
2018 Unfenced Unburned cover_allwillow 0.47 0.34 0.00 0.41 1.00 61 100
2018 Unfenced Burned cover_allwillow 0.15 0.34 0.00 0.03 1.00 8 100
2018 Fenced Unburned cover_allwillow 0.66 0.33 0.14 0.64 1.00 17 100
2018 Fenced Burned cover_allwillow 0.39 0.29 0.05 0.26 1.00 13 100

Combined core and Non-core range

# ### AOV: cov ~ time class * condition
# cull to just bl, 13, 18 and WC WNC
canopy.all.willow.sum.5yr.wcwnc <- canopy.all.willow.sum %>%
  filter(site_type != "WK") %>% 
  filter(time_class %in% c("BL", "2013","2018")) %>% 
  mutate(time_class = fct_drop(time_class))
# nest by TC
willow.cov.tc.nest <- canopy.all.willow.sum.5yr.wcwnc %>%
  nest_by(time_class)

### create AOV models
mods.aov.willcov.tcxZcond2 <- willow.cov.tc.nest %>%
  mutate(mod = list(aov(cover_transformed ~ z_cond2 + Error(1/site_id), data = data)))
## nest by z cond & TC
willow.cov.zc.wcwnc.nest <- canopy.all.willow.sum.5yr.wcwnc %>% 
  nest_by(z_cond2)

willow.cov.zc.wcwnc.nest <- willow.cov.zc.wcwnc.nest %>% 
  mutate(mod = list(aov(cover_transformed ~ time_class + Error(1/site_id), data = data)))
# GT table
willow.cov.zc.wcwnc.nest %>% 
  mutate(aov_summary = list(summary(mod))) %>% 
  mutate(aov_tidy = list(broom::tidy(mod))) %>% 
  unnest(aov_tidy) %>%
  mutate(across(where(is.numeric), round, 3)) %>% 
  select(z_cond2,term,df,sumsq,meansq,statistic,p.value) %>% 
  gt::gt() %>% 
  gt::tab_header(title = "aov: cover (transformed) ~ condition + Error(1/site_id)", subtitle = "analysis on transformed scale")
aov: cover (transformed) ~ condition + Error(1/site_id)
analysis on transformed scale
term df sumsq meansq statistic p.value
WC-BF
time_class 2 71.862 35.931 11.251 0.000
Residuals 30 95.811 3.194 NA NA
WC-BG
time_class 2 0.361 0.181 0.041 0.960
Residuals 16 71.135 4.446 NA NA
WC-UF
time_class 2 36.998 18.499 7.158 0.002
Residuals 43 111.128 2.584 NA NA
WC-UG
time_class 2 8.453 4.227 1.097 0.339
Residuals 73 281.261 3.853 NA NA
WNC-UG
time_class 2 13.411 6.706 2.254 0.111
Residuals 88 261.813 2.975 NA NA 
# write csv
willow.cov.zc.wcwnc.nest %>%
  mutate(aov_summary = list(summary(mod))) %>%
  mutate(aov_tidy = list(broom::tidy(mod))) %>%
  unnest(aov_tidy) %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  select(z_cond2,term,df,sumsq,meansq,statistic,p.value) %>%
  write_csv("output/tables/willow_cover/frequentist/aov_willcov_zcond2.csv")
# ### Table for manuscript
# extract model info, tidy
wcov.pval <- willow.cov.zc.wcwnc.nest %>%
  mutate(aov_summary = list(summary(mod))) %>%
  mutate(aov_tidy = list(broom::tidy(mod))) %>%
  unnest(aov_tidy) %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  select(z_cond2,term,df,sumsq,meansq,statistic,p.value) %>% 
  filter(term == "time_class") 


## calc empirical mean of cover then join in results from freq aov
# cover transformed ~ 
canopy.all.willow.sum.5yr.wcwnc %>% 
  group_by(time_class, z_cond2) %>% 
  summarise(cov.mean = mean(cover_allwillow, na.rm=TRUE)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = time_class,
              names_glue = "Cover-{time_class}",
              values_from = cov.mean) %>% 
  left_join(.,wcov.pval, by="z_cond2") %>% 
  mutate(across(where(is.numeric), round, 2)) %>% 
  gt()
z_cond2 Cover-BL Cover-2013 Cover-2018 term df sumsq meansq statistic p.value
WC-BF 2.89 12.00 38.63 time_class 2 71.86 35.93 11.25 0.00
WC-BG 2.82 10.11 15.14 time_class 2 0.36 0.18 0.04 0.96
WC-UF 24.04 41.07 65.54 time_class 2 37.00 18.50 7.16 0.00
WC-UG 26.80 22.69 37.08 time_class 2 8.45 4.23 1.10 0.34
WNC-UG 36.99 43.00 55.57 time_class 2 13.41 6.71 2.25 0.11 
# wcov.pval %>%
#   write_csv("output/tables/willow_cover/frequentist/aov_willcov_pval_zcond2_wcwnc.csv")
#### Frequentist AOV. Single model: WNC ~ TC*zcond2
aov_model_wcov_wcwnc <- aov(cover_transformed ~ time_class*z_cond2 + Error(1/site_id), data = cov.wc.wnc)

Model summary:
The repeated-measures ANOVA (formula: cover_transformed ~ time_class * z_cond2 + Error(1/site_id)) suggests that:

  • The main effect of time_class is significant and medium (F(2, 250) = 10.53, p < .001; Eta2 (partial) = 0.08, 90% CI [0.03, 0.13])
  • The main effect of z_cond2 is significant and large (F(4, 250) = 19.40, p < .001; Eta2 (partial) = 0.24, 90% CI [0.16, 0.30])
  • The interaction between time_class and z_cond2 is not significant and small (F(8, 250) = 1.87, p = 0.065; Eta2 (partial) = 0.06, 90% CI [0.00, 0.12])

Effect sizes were labelled following Field’s (2013) recommendations.

broom::tidy(aov_model_wcov_wcwnc) %>% 
  mutate(across(.cols = where(is.numeric),.fns = round, 3)) %>% 
  gt() %>% 
  tab_header(title = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)", subtitle = "Combined WC and WNC")
cover_transformed ~ time_class + z_cond2 + Error(1/site_id)
Combined WC and WNC
stratum term df sumsq meansq statistic p.value
Within time_class 2 69.203 34.601 10.534 0.000
Within z_cond2 4 254.848 63.712 19.397 0.000
Within time_class:z_cond2 8 49.106 6.138 1.869 0.065
Within Residuals 250 821.147 3.285 NA NA 
## save as table (.csv)
# broom::tidy(aov_model_wc) %>% 
#   mutate(across(.cols = where(is.numeric),.fns = round, 3)) %>% 
#   write_csv("output/tables/aov_wc_ht_TC_fenced.csv")
#_Estimated marginal means_
emmeans (aov_model_wcwnc,  ~ time_class) %>% 
  broom::tidy() %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  gt()
time_class estimate std.error df statistic p.value
BL 118.103 10.145 275 11.642 0
2013 136.237 9.707 275 14.035 0
2018 151.716 9.737 275 15.581 0 
emmeans(aov_model_wcov_wcwnc,  ~ z_cond2|time_class) %>% 
  pairs() %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept=0), color = "grey40", size=1) +
  labs(x = "Contrast", y = "", caption = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)") 

emmeans(aov_model_wcov_wcwnc,  ~ z_cond2|time_class) %>% 
  pairs() %>% 
  as_tibble() %>%
  dplyr::arrange(time_class, contrast) %>% 
  mutate(across(where(is.numeric), round,2)) %>%
  filter(p.value <0.1) %>% 
  gt() %>%
  tab_header(title = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)", subtitle = "Only contrasts with p <0.1 displayed. P value adjustment: tukey method for comparing estimates")
cover_transformed ~ time_class + z_cond2 + Error(1/site_id)
Only contrasts with p <0.1 displayed. P value adjustment: tukey method for comparing estimates
contrast time_class estimate SE df t.ratio p.value
(WC-BF) - (WC-UF) BL -2.66 0.86 250 -3.08 0.02
(WC-BF) - (WC-UG) BL -2.64 0.80 250 -3.31 0.01
(WC-BF) - (WNC-UG) BL -3.41 0.76 250 -4.51 0.00
(WC-BF) - (WC-UF) 2013 -2.68 0.67 250 -4.01 0.00
(WC-BF) - (WNC-UG) 2013 -2.87 0.61 250 -4.68 0.00
(WC-BG) - (WC-UF) 2013 -2.82 0.72 250 -3.90 0.00
(WC-BG) - (WNC-UG) 2013 -3.01 0.67 250 -4.48 0.00
(WC-UG) - (WNC-UG) 2013 -1.41 0.49 250 -2.87 0.04
(WC-BF) - (WC-BG) 2018 2.52 0.81 250 3.10 0.02
(WC-BG) - (WC-UF) 2018 -3.91 0.78 250 -5.03 0.00
(WC-BG) - (WC-UG) 2018 -2.10 0.72 250 -2.91 0.03
(WC-BG) - (WNC-UG) 2018 -3.33 0.72 250 -4.64 0.00
(WC-UF) - (WC-UG) 2018 1.80 0.55 250 3.26 0.01
(WC-UG) - (WNC-UG) 2018 -1.22 0.46 250 -2.63 0.07 
cov.wc.wnc %>% 
  ggplot() +
  geom_boxplot(aes(x=z_cond2, y = cover_allwillow), fill="grey80") +
  facet_wrap(~time_class) +
  theme_minimal() +
  labs(x="Condition", y = "Mean cover", caption = "WC & WNC plots, all willow species")

cov.wc.wnc %>% 
  group_by(time_class, z_cond2) %>% 
  descr(cover_allwillow, stats = "common") %>% 
  tb() %>% 
  mutate(across(where(is.numeric), round,2)) %>%
  gt() %>% 
  tab_header(title="Summary statistics", subtitle = "WC & WNC plots, all willow species")
Summary statistics
WC & WNC plots, all willow species
time_class z_cond2 variable mean sd min med max n.valid pct.valid
BL WC-BF cover_allwillow 0.03 0.03 0.00 0.02 0.09 7 100
BL WC-BG cover_allwillow 0.03 NA 0.03 0.03 0.03 1 100
BL WC-UF cover_allwillow 0.24 0.26 0.01 0.17 1.00 12 100
BL WC-UG cover_allwillow 0.27 0.23 0.00 0.25 1.00 20 100
BL WNC-UG cover_allwillow 0.37 0.28 0.03 0.32 1.00 32 100
2013 WC-BF cover_allwillow 0.12 0.21 0.00 0.04 0.79 13 100
2013 WC-BG cover_allwillow 0.10 0.22 0.01 0.02 0.72 10 100
2013 WC-UF cover_allwillow 0.41 0.29 0.02 0.40 1.00 17 100
2013 WC-UG cover_allwillow 0.23 0.22 0.00 0.18 1.00 27 100
2013 WNC-UG cover_allwillow 0.43 0.32 0.02 0.37 1.00 27 100
2018 WC-BF cover_allwillow 0.39 0.29 0.05 0.26 1.00 13 100
2018 WC-BG cover_allwillow 0.15 0.34 0.00 0.03 1.00 8 100
2018 WC-UF cover_allwillow 0.66 0.33 0.14 0.64 1.00 17 100
2018 WC-UG cover_allwillow 0.37 0.32 0.00 0.29 1.00 29 100
2018 WNC-UG cover_allwillow 0.56 0.33 0.01 0.53 1.00 32 100 
emmeans(aov_model_wcov_wcwnc,  ~ time_class) %>% 
  pairs(reverse=TRUE) %>%
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept=0), color = "grey40", size=1) +
  labs(x = "Contrast", y = "", caption = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)")

emmeans(aov_model_wcov_wcwnc,  ~ time_class) %>% 
  pairs(reverse=TRUE) %>% 
  as_tibble() %>% 
  # dplyr::arrange(p.value) %>% 
  mutate(across(where(is.numeric), round,2)) %>% 
  gt() %>% 
  tab_header(title = "cover_transformed ~ time_class + Error(1/site_id)", subtitle = "P value adjustment: tukey method for comparing estimates")
cover_transformed ~ time_class + Error(1/site_id)
P value adjustment: tukey method for comparing estimates
contrast estimate SE df t.ratio p.value
2013 - BL 0.37 0.46 250 0.81 0.7
2018 - BL 1.50 0.46 250 3.23 0.0
2018 - 2013 1.12 0.29 250 3.88 0.0 
emmeans(aov_model_wcov_wcwnc,  ~ z_cond2) %>% 
  pairs(reverse=TRUE) %>% 
  plot() +
  theme_minimal() +
  geom_vline(aes(xintercept=0), color = "grey40", size=1) +
  labs(x = "Contrast", y = "", caption = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)")

emmeans(aov_model_wcov_wcwnc,  ~ z_cond2) %>% 
  pairs(reverse=TRUE) %>% 
  as_tibble() %>% 
  dplyr::arrange(p.value) %>% 
  mutate(across(where(is.numeric), round,2)) %>% 
  gt() %>% 
  tab_header(title = "cover_transformed ~ time_class + z_cond2 + Error(1/site_id)", subtitle = "P value adjustment: tukey method for comparing estimates")
cover_transformed ~ time_class + z_cond2 + Error(1/site_id)
P value adjustment: tukey method for comparing estimates
contrast estimate SE df t.ratio p.value
(WNC-UG) - (WC-BF) 2.36 0.38 250 6.20 0.00
(WC-UF) - (WC-BF) 2.24 0.43 250 5.25 0.00
(WNC-UG) - (WC-BG) 3.05 0.70 250 4.38 0.00
(WC-UF) - (WC-BG) 2.93 0.72 250 4.06 0.00
(WNC-UG) - (WC-UG) 1.14 0.28 250 4.00 0.00
(WC-UG) - (WC-BF) 1.22 0.39 250 3.13 0.02
(WC-UG) - (WC-UF) -1.02 0.34 250 -2.97 0.03
(WC-UG) - (WC-BG) 1.91 0.70 250 2.73 0.05
(WC-BG) - (WC-BF) -0.69 0.75 250 -0.92 0.89
(WNC-UG) - (WC-UF) 0.12 0.33 250 0.36 1.00

Core range

Model summary, core winter range:
The repeated-measures ANOVA (formula: cover_transformed ~ time_class * z_cond2 + Error(1/site_id)) suggests that:

  • The main effect of time_class is significant and medium (F(2, 162) = 9.48, p < .001; Eta2 (partial) = 0.10, 90% CI [0.04, 0.18])
  • The main effect of z_cond2 is significant and large (F(3, 162) = 15.58, p < .001; Eta2 (partial) = 0.22, 90% CI [0.13, 0.31])
  • The interaction between time_class and z_cond2 is not significant and medium (F(6, 162) = 2.12, p = 0.054; Eta2 (partial) = 0.07, 90% CI [0.00, 0.12])

Effect sizes were labelled following Field’s (2013) recommendations.

Non-core range

Model summary, non-core winter range:
The repeated-measures ANOVA (formula: cover_transformed ~ time_class + Error(1/site_id)) suggests that:

  • The main effect of time_class is not significant and small (F(2, 88) = 2.25, p = 0.111; Eta2 (partial) = 0.05, 90% CI [0.00, 0.13])

Effect sizes were labelled following Field’s (2013) recommendations.

# _Separate analyses by condition class._
### create subsets of data for modeling
# pull out WK, note different years of data
cov.wk <- canopy.all.willow.sum %>%
  filter(site_type == "WK") %>%
  mutate(time_class = as.character(time_class)) %>%
  filter(time_class == "2016" | time_class == "2017" | time_class == "2018") %>%
  mutate(time_class = as.factor(time_class))

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

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

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

Kawuneeche Valley

# model for wk
aov_model_wcov_wk <- aov(cover_transformed ~ time_class + Error(1/site_id), data = cov.wk)

Model summary, Kawuneeche Valley:
The repeated-measures ANOVA (formula: cover_transformed ~ time_class + Error(1/site_id)) suggests that:

  • The main effect of time_class is not significant and small (F(2, 20) = 0.12, p = 0.887; Eta2 (partial) = 0.01, 90% CI [0.00, 0.08])

Effect sizes were labelled following Field’s (2013) recommendations.

aov_model_wcov_wk %>% 
  report::report() %>% 
  as.data.frame() %>% 
  gt()
Group Parameter Sum_Squares df Mean_Square F p Eta2_partial Eta2_partial_CI_low Eta2_partial_CI_high
Within time_class 0.5346485 2 0.2673242 0.1208116 0.8868427 0.01193694 0 0.0802567
Within Residuals 44.2547459 20 2.2127373 NA NA NA NA NA 
aov_model_wcov_wk %>% 
  report::report_model() %>% 
  as.data.frame() %>% 
  gt()
.
repeated-measures ANOVA (formula: cover_transformed ~ time_class + Error(1/site_id))

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 18

report::report_table(sessionInfo()) %>% 
  as_tibble() %>% 
  gt() %>% 
  tab_header(title="Package Information")
Package Information
Package Version Reference
dplyr 1.0.4 Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2021). dplyr: A Grammar of Data Manipulation. R package version 1.0.4. https://CRAN.R-project.org/package=dplyr
DT 0.17 Yihui Xie, Joe Cheng and Xianying Tan (2021). DT: A Wrapper of the JavaScript Library 'DataTables'. R package version 0.17. https://CRAN.R-project.org/package=DT
emmeans 1.5.4 Russell V. Lenth (2021). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.5.4. 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.3 H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
ggpubr 0.4.0 Alboukadel Kassambara (2020). ggpubr: 'ggplot2' Based Publication Ready Plots. R package version 0.4.0. https://CRAN.R-project.org/package=ggpubr
glue 1.4.2 Jim Hester (2020). glue: Interpreted String Literals. R package version 1.4.2. https://CRAN.R-project.org/package=glue
gt 0.2.2 Richard Iannone, Joe Cheng and Barret Schloerke (2020). gt: Easily Create Presentation-Ready Display Tables. R package version 0.2.2. 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.31 Yihui Xie (2021). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.31.
lubridate 1.7.9.2 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/.
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/.
readr 1.4.0 Hadley Wickham and Jim Hester (2020). readr: Read Rectangular Text Data. R package version 1.4.0. https://CRAN.R-project.org/package=readr
rstatix 0.7.0 Alboukadel Kassambara (2021). rstatix: Pipe-Friendly Framework for Basic Statistical Tests. R package version 0.7.0. https://CRAN.R-project.org/package=rstatix
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 0.9.8 Dominic Comtois (2020). summarytools: Tools to Quickly and Neatly Summarize Data. R package version 0.9.8. https://CRAN.R-project.org/package=summarytools
tibble 3.0.6 Kirill Müller and Hadley Wickham (2021). tibble: Simple Data Frames. R package version 3.0.6. https://CRAN.R-project.org/package=tibble
tidyr 1.1.2 Hadley Wickham (2020). tidyr: Tidy Messy Data. R package version 1.1.2. https://CRAN.R-project.org/package=tidyr
tidyverse 1.3.0 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.5.1 Simon Garnier (2018). viridis: Default Color Maps from 'matplotlib'. R package version 0.5.1. https://CRAN.R-project.org/package=viridis
viridisLite 0.3.0 Simon Garnier (2018). viridisLite: Default Color Maps from 'matplotlib' (Lite Version). R package version 0.3.0. https://CRAN.R-project.org/package=viridisLite

References

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.