Code author:
E. Gage
mycologica@gmail.com
https://github.com/egage/EVMP>
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.”
They are provided with important caveats and reservations. The frequentist analyses here fundamentally differ in philosophical foundation, basic inference, and interpretation than the Bayesian analyses presented in the report. No attempt is made here made to reconcile differences. Caution in interpreting results presented here is recommended, particularly in regards to concepts like p-values, confidence intervals, and statistical significance. These concepts are very commonly misunderstood (often by incorrectly assuming a Bayesian interpretation), in part why a Bayesian approach was utilized for the main analyses in the document.
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 wrangling used to transform raw files to the derived data ingested for analyses here.
Per guidance from RMNP, only willow macroplot data are analyzed here. Based on data structure, a repeated measures Analysis of Variance (ANOVA) approach is used.
No significant outliers. This was checked by visualizing the data
Normality of the dependent variable. This was evaluated using the Shapiro-Wilk normality test and visual inspection of QQ plots.
Equal variance between groups. This was checked using the Mauchly’s test of sphericity
Per methods used in Zeigenfuss and Johnson, the minimum non-zero proportion of shrub cover is substituted for areas with no shrub cover and (1 - this minimum) is used for sites with a proportion of 1 so that a log of proportion of cover can 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")
willow <- willow %>%
select(-c(contains("UTM")))
## 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 to 5 yr sampling interval, chr to factor conversions
willow.ht <- willow.ht %>%
filter(timeClass == "BL" | timeClass == "2013" | timeClass == "2018") %>%
mutate(timeClass = as.factor(timeClass)) %>%
mutate(timeClass = fct_relevel(timeClass, "2018", "2013", "BL")) %>%
mutate(FENCED = as.factor(FENCED))
## clean names
willow.ht <- willow.ht %>%
janitor::clean_names()
## filter out just the willow species
willow.ht <- willow.ht %>%
filter(str_detect(species_code, pattern = "SA"))
## filter out 0
willow.ht <- willow.ht %>%
filter(plant_ht_cm > 0)
## reorder the time class
willow.ht <- willow.ht %>%
mutate(time_class = fct_rev(time_class))
## factor reverse: fenced
willow.ht <- willow.ht %>%
mutate(fenced = fct_rev(fenced))## create separate df for WC, WNC, combined WC & WNC, and KV
willow.ht.wc <- willow.ht %>%
filter(site_type == "WC")
# non core
willow.ht.nc <- willow.ht %>%
# distinct(site_type)
filter(site_type == "WNC")
# combined WC and WNC
willow.ht.wc.wnc <- willow.ht %>%
# distinct(site_type)
filter(site_type == "WNC" | site_type == "WC")
# kw
willow.ht.kw <- willow.ht %>%
filter(site_type == "WK") Willow Macroplot Data
## # A tibble: 4 x 7
## stratum term df sumsq meansq statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Within time_class 2 656403. 328202. 74.9 0
## 2 Within fenced 1 286047. 286047. 65.3 0
## 3 Within time_class:fenced 2 464675. 232337. 53.0 0
## 4 Within Residuals 2494 10930705. 4383. NA NA## 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(plant_ht_cm ~ 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(plant_ht_cm ~ time_class*fenced + 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(plant_ht_cm ~ time_class*fenced + 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)")
# write csv
# 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) %>%
# write_csv("output/aov_willHt_tc_fen.csv")#### Post-hoc tests and multiple pairwise comparisons between groups to identify which groups are different# 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 Tables
#### summary tables: time class and burned
cov.wc.wnc %>%
group_by(time_class, fenced) %>%
summarytools::descr(cover_allwillow,stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Combined WC and WNC macroplot data: willow cover - fenced + time_class")| Combined WC and WNC macroplot data: willow cover - fenced + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | fenced | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | 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 |
## time class, fenced and burned
cov.wc.wnc %>%
group_by(time_class, fenced, burned) %>%
summarytools::descr(cover_allwillow, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Combined WC and WNC macroplot data: willow cover ~ burned + fenced + time_class")| Combined WC and WNC macroplot data: willow cover ~ burned + fenced + time_class | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| time_class | fenced | burned | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | Unfenced | Unburned | cover_allwillow | 0.33 | 0.27 | 0.00 | 0.27 | 1.00 | 52 | 100 |
| BL | Unfenced | Burned | cover_allwillow | 0.03 | NA | 0.03 | 0.03 | 0.03 | 1 | 100 |
| BL | Fenced | Unburned | cover_allwillow | 0.24 | 0.26 | 0.01 | 0.17 | 1.00 | 12 | 100 |
| BL | Fenced | Burned | cover_allwillow | 0.03 | 0.03 | 0.00 | 0.02 | 0.09 | 7 | 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 |
#### Post-hoc tests and multiple pairwise comparisons between groups to identify which groups are different# add transform following method in Zeigenfuss 2015. Treatment of 0 and 1 values in
canopy.all.willow.sum <- canopy.all.willow.sum %>%
mutate(cover_transformed = cover_allwillow/100) %>%
mutate(cover_transformed = case_when(cover_transformed == 1 ~ 1-0.0883573,
cover_transformed == 0 ~ 0.0883573,
TRUE ~ cover_transformed)) %>%
mutate(cover_transformed = log(cover_transformed/(1-cover_transformed)))# cull to just bl, 13, 18
canopy.all.willow.sum.5yr <- canopy.all.willow.sum %>%
filter(time_class %in% c("BL", "2013","2018")) %>%
mutate(time_class = fct_drop(time_class))
# canopy.all.willow.sum.5yr[,1:8]# nest by TC
willow.cov.tc.nest <- canopy.all.willow.sum.5yr %>%
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.nest <- canopy.all.willow.sum.5yr %>%
nest_by(z_cond2)
## remove the KW plots
mods.aov.willcov.z_cond2 <- willow.cov.zc.nest %>%
filter(z_cond2 == "WC-BF" | z_cond2 == "WC-BG" | z_cond2 == "WC-UG" | z_cond2 == "WC-UF" | z_cond2 == "WNC-UG") %>%
mutate(mod = list(aov(cover_transformed ~ time_class + Error(1/site_id), data = data)))# GT table
mods.aov.willcov.z_cond2 %>%
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
mods.aov.willcov.z_cond2 %>%
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")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))| 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.258 | 0 |
| Within | z_cond2 | 4 | 254.848 | 63.712 | 18.888 | 0 |
| Within | Residuals | 258 | 870.254 | 3.373 | NA | NA |
Summary: The repeated-measures ANOVA (formula: plant_ht_cm ~ time_class * fenced + Error(1/site_id)) suggests that:
Effect sizes were labelled following Field’s (2013) recommendations.
canopy.all.willow.sum.wcwnc %>%
group_by(time_class, z_cond2) %>%
nest()## # A tibble: 15 x 3
## # Groups: time_class, z_cond2 [15]
## time_class z_cond2 data
## <fct> <chr> <list>
## 1 2018 WC-UF <tibble [17 x 10]>
## 2 2018 WC-UG <tibble [29 x 10]>
## 3 2018 WC-BF <tibble [13 x 10]>
## 4 2018 WC-BG <tibble [8 x 10]>
## 5 2018 WNC-UG <tibble [32 x 10]>
## 6 2013 WC-UF <tibble [17 x 10]>
## 7 2013 WC-UG <tibble [29 x 10]>
## 8 2013 WC-BF <tibble [13 x 10]>
## 9 2013 WC-BG <tibble [10 x 10]>
## 10 2013 WNC-UG <tibble [32 x 10]>
## 11 BL WC-UF <tibble [12 x 10]>
## 12 BL WC-UG <tibble [20 x 10]>
## 13 BL WC-BF <tibble [8 x 10]>
## 14 BL WC-BG <tibble [6 x 10]>
## 15 BL WNC-UG <tibble [32 x 10]># library(emmeans)
emmeans (aov_model_wcwnc, ~ time_class + z_cond2) %>%
broom::tidy() %>%
mutate(across(where(is.numeric), round, 3)) %>%
gt()| time_class | z_cond2 | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|---|
| BL | WC-BF | -2.502 | 0.350 | 258 | -7.152 | 0.000 |
| 2013 | WC-BF | -2.234 | 0.341 | 258 | -6.552 | 0.000 |
| 2018 | WC-BF | -1.218 | 0.339 | 258 | -3.592 | 0.000 |
| BL | WC-BG | -3.469 | 0.429 | 258 | -8.086 | 0.000 |
| 2013 | WC-BG | -3.201 | 0.391 | 258 | -8.193 | 0.000 |
| 2018 | WC-BG | -2.185 | 0.399 | 258 | -5.475 | 0.000 |
| BL | WC-UF | -0.377 | 0.317 | 258 | -1.188 | 0.236 |
| 2013 | WC-UF | -0.110 | 0.318 | 258 | -0.344 | 0.731 |
| 2018 | WC-UF | 0.907 | 0.316 | 258 | 2.869 | 0.004 |
| BL | WC-UG | -1.504 | 0.288 | 258 | -5.215 | 0.000 |
| 2013 | WC-UG | -1.237 | 0.291 | 258 | -4.243 | 0.000 |
| 2018 | WC-UG | -0.220 | 0.286 | 258 | -0.771 | 0.441 |
| BL | WNC-UG | -0.324 | 0.268 | 258 | -1.207 | 0.229 |
| 2013 | WNC-UG | -0.056 | 0.294 | 258 | -0.191 | 0.849 |
| 2018 | WNC-UG | 0.960 | 0.285 | 258 | 3.372 | 0.001 |
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 08