Appendix C. EVMP Upland Statistical Analyses
- Introduction
- Methods
- Results
- Individual species
- Upland Shrub Cover, Purshia tridentata - Combined Winter Range
- Upland Shrub Cover, Purshia tridentata - Core Winter Range
- Upland Shrub Cover, Purshia tridentata - Noncore Winter Range
- Upland Shrub Cover, Ribes cereum - Combined Winter Range
- Upland Shrub Cover, Ribes cereum - Core Winter Range
- Upland Shrub Cover, Ribes cereum - Noncore Winter Range
- Additional plots
- Session info
Code author: E. Gage
mycologica@gmail.com
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")# purrr plots
tmp1 <- function(df){
dcanopy.all.willow.sum %>%
group_nest(SITE_TYPE) %>%
mutate(plots = map2(.y = SITE_TYPE, .x = data, ~{ggplot(data = .x) +
geom_boxplot(aes(y = cover.allwillow, x = timeClass), fill="Gray80") +
geom_hline(aes(xintercept=0.31)) +
# ggtitle(paste0("Shrub cover, site: ", .y)) +
labs(title = paste0("Site type: ", .y), x = "", y = "Willow cover (%)") +
theme_minimal()}))
tmp1 %>%
pull(plots)Introduction
This document contains code for running statistical analyses of of Elk Vegetation Management Plan (EVMP) data collected through the 2018 sampling season. Analyses are limited to data from plots established in upland communities in areas of the elk winter range (core and noncore area, see Zeigenfuss et al. 2011 for details) 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
Background
The code and results presented here build on derived data produced in a separate code document focused on ingesting, compiling, and cleaning raw data provided by RMNP staff (Appendix B). Extensive “munging” of raw files was required, so the decision was made to separate code used to pre-process raw data from code used for the analyses presented here. Likewise, separate code documents were developed for other distinct elements of the larger EVMP (e.g., aspen and upland plots). In general, packages included in the “Tidyverse” ecosystem of R packages were used for data transformation and visualization (Wickham et al. 2019), although specialized packages particular to specific tasks (e.g., the “bayesplot” plot for visualization of Bayesian posterior distributions) were also used (Gabry et al. 2019).
Bayesian 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) 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 (Lemoine 2019). Zero-bounded continuous response variables (e.g., shrub height) were fit with a Gamma distribution, while proportion data (e.g., cover) were modeled using a beta distribution (Ferrari and Cribari-Neto 2004) that 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 we have and a specific parameter value
- \(P(\theta)\): Our view of the possible values of our model parameters before observing any actual data, i.e. the prior
- \(P(y)\): The unconditional probability of the outcome, i.e. the probability after considering all possible parameter. The objective is to estimate the posterior density of parameter \(P(\theta|y,X)\), typically described in terms of a point estimate of the expected value of posterior density (e.g., the median of the distribution) and a measure of the variance. To solve for \(P(\theta|y,X)\) analytically, the data likelihood \(P(y|\theta,X)\), prior \(P(\theta)\), and marginal likelihood of outcome \(P(y)\) are needed, but because there is usually no closed-form analystical solution, Markov chain Monte Carlo methods are used to numerically solve the posterior density by directly generating random draws of parameters. The general steps in a Bayesian analysis followed included: Specify a joint distribution for the outcome(s) and all the unknowns, which proportional to a posterior distribution of the unknowns conditional on the observed data
Use Markov Chain Monte Carlo (MCMC) to draw from posterior distribution.
*Evaluation of model fit and revise the model as appropriate.
Shrub height modeling
Bayesian repeated measures analyses of shrub height were fit separately for combined core and non-core winter range, core winter range (WC) and non-core winter range plots (WNC). Separate models were also fit to different groups of species. Shrub height models were fit using weekly informative priors and a gamma distribution with the “stan_glmer” function in the “rstanarm” package (Goodrich et al. 2020)(Brilleman et al. 2018). Gamma regression is commonly used when the response variable is continuous and positive.
For the simplest case a GLM for a continuous outcome is simply a linear model and the likelihood for one observation is a conditionally normal probability density function.
\[\frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{1}{2}
\left(\frac{y - \mu}{\sigma}\right)^2},\] where \(\mu = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) is a linear predictor and \(\sigma\) is the standard deviation of the error in predicting the outcome, \(y\). A linear predictor \(\eta = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) can more generally be related to the conditional mean of the outcome via a link function \(g\). This maps the range of values on which the outcome is defined. For the gamma distribution used in the modeling of willow height, the log link function was used. The likelihood is the product of the likelihood contributions of individual observations.
Bayesian analysis requires specifying prior distributions \(f(\alpha)\) and \(f(\boldsymbol{\beta})\) for the intercept and vector of regression coefficients. Prior distributions were represented by normal distributions with a mean zero and a small standard deviation (scale)). The joint posterior distribution for \(\alpha\) and \(\boldsymbol{\beta}\) is proportional to the product of the priors and the \(N\) likelihood contributions:
\[f\left(\boldsymbol{\beta} | \mathbf{y},\mathbf{X}\right) \propto f\left(\alpha\right) \times \prod_{k=1}^K f\left(\beta_k\right) \times \prod_{i=1}^N {f(y_i|\eta_i)},\]
where \(\mathbf{X}\) is the matrix of predictors and \(\eta\) the linear predictor, i.e., the posterior distribution that stan_glmer draws from using MCMC.
Shrub cover modeling
Shrub cover data were modeled using the beta distribution as the likelihood for the data, \[ f(y_i | a, b) = \frac{y_i^{(a-1)}(1-y_i)^{(b-1)}}{B(a,b)} \] where \(B(\cdot)\) is the beta function. The shape parameters for the distribution are \(a\) and \(b\) and enter into the model according to the following transformations, \[ a = \mu\cdot\phi \\ b = (1-\mu)\cdot\phi \]
If \(g_1(\cdot)\) is some link function, the specification of the shape parameters, \(\mu = g_1^{-1}(\mathbf{X}\boldsymbol{\beta})\), where \(\boldsymbol{X}\) is a \(N\times K\) dimensional matrix of predictors, and \(\boldsymbol{\beta}\), is a \(K\) dimensional vector of parameters associated with each predictor. \(\phi\) is a scalar parameter. With the shape parameter values included, the likelihood takes the form: \[ f(y_i | \mu, \phi) = \frac{y_i^{(\mu\phi-1)}(1-y_i)^{((1-\mu)\phi-1)}}{B(\mu\phi,(1-\mu)\phi)} \]
Bayesian analysis requires specifying prior distributions \(f(\boldsymbol{\beta})\) and \(f(\phi)\) for the vector of regression coefficients and \(\phi\). When modeling \(\phi\) with a linear predictor, a Bayesian analysis requires specifying the prior distributions \(f(\boldsymbol{\beta})\) and \(f(\boldsymbol{\gamma})\).
With only a single set of regressors, the posterior distribution of \(\boldsymbol{\beta}\) and \(\phi\) is proportional to the product of the likelihood contributions, the \(K\) priors on the \(\beta_k\) parameters, and \(\phi\), \[ f(\boldsymbol{\beta},\phi|\mathbf{y},\mathbf{X}) \propto \prod_{i=1}^N f(y_i | a, b) \times \prod_{k=1}^K f(\beta_k) \times f(\phi) \]
When using two sets of regressors, the posterior distribution of \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) is proportional to the product of the likelihood contribution, the \(K\) priors on the \(\beta_k\) parameters, and the \(J\) priors on the \(\gamma_j\) parameters,
\[ f(\boldsymbol{\beta},\boldsymbol{\gamma}|\mathbf{y},\mathbf{X}) \propto \prod_{i=1}^N f(y_i | a, b) \times \prod_{k=1}^K f(\beta_k) \times \prod_{j=1}^J f(\gamma_j) \]
### Data import and munging
## read in derived data
upl.allshrubs <- read_csv("./output/exported_data/upland_LI_allShrubsPooled_20200814.csv")
upl.byspp <- read_csv("./output/exported_data/upland_LI_byShrubSpp_20200814.csv")
### upland pooled species
cov.allshrubs <- upl.allshrubs %>%
filter(!is.na(perc_cover)) %>%
mutate(time_class = case_when(yr == 2007 ~ "BL",
TRUE ~ as.character(yr))) %>%
mutate(timeClass = as.factor(time_class)) %>%
mutate(site_type = as.factor(site_type)) %>%
mutate(time_class = fct_relevel(timeClass, "2018", "2013", "BL"))
## reorder the time class
cov.allshrubs <- cov.allshrubs %>%
mutate(time_class = fct_rev(time_class))
### upland by species
cov.byspp <- upl.byspp %>%
filter(!is.na(perc_cover)) %>%
mutate(time_class = case_when(yr == 2007 ~ "BL",
TRUE ~ as.character(yr))) %>%
mutate(time_class = as.factor(time_class)) %>%
mutate(site_type = as.factor(site_type)) %>%
mutate(time_class = fct_relevel(time_class, "2018", "2013", "BL"))
## reorder the time class
cov.byspp <- cov.byspp %>%
mutate(time_class = fct_rev(time_class))
# %>%
# mutate(timeClass = case_when(yr == 2008 ~ "BL",
# yr == 2009 ~ "BL",
# TRUE ~ as.character(yr))) %>%
# mutate(timeClass = as.factor(timeClass)) %>%
# mutate(timeClass = fct_relevel(timeClass, "2018", "2013", "BL"))## create seperate subsets of ac and anc
cov.allshrubs %>%
glimpse()
cov.allshrubs %>%
names()
# cov.allshrubs %>%
# tabyl()# join in site info
## join in the site info
site.info.clean <- read_csv("./data/EVMP_derived/site_info_clean.csv")
site.info.cleannames <- site.info.clean %>%
clean_names() %>%
select(-starts_with('utm')) %>%
select(-c(wilderness,east_west)) %>%
distinct() %>%
mutate(burned = case_when(burned == "Not burned" ~ "Unburned",
is.na(burned) ~ "Unburned",
TRUE ~ burned))
cov.allshrubs <- left_join(cov.allshrubs, site.info.cleannames, by = "site_id")
cov.byspp <- left_join(cov.byspp, site.info.cleannames, by = "site_id")Results
# purrr plots
list.pl.allshrub.wcwnc <- cov.allshrubs %>%
group_nest(site_type) %>%
mutate(plots = map2(.y = site_type, .x = data, ~{ggplot(data = .x) +
geom_boxplot(aes(y = perc_cover, x = time_class), fill="Gray80") +
# geom_hline(aes(xintercept=0.31)) +
# ggtitle(paste0("Shrub cover, site: ", .y)) +
labs(title = paste0("Site type: ", .y), x = "", y = "All shrub cover (%)") +
theme_minimal()}))
list.pl.allshrub.wcwnc %>%
pull(plots)## [[1]]
##
## [[2]]
# purrr plots
list.pl.allshrub.wc.f <- cov.allshrubs %>%
filter(site_type == "core range") %>%
group_nest(fenced) %>%
mutate(plots = map2(.y = fenced, .x = data, ~{ggplot(data = .x) +
geom_boxplot(aes(y = perc_cover, x = time_class), fill="Gray80") +
# geom_hline(aes(xintercept=0.31)) +
# ggtitle(paste0("Shrub cover, site: ", .y)) +
labs(title = paste0("Site type: ", .y), x = "", y = "All shrub cover (%)") +
theme_minimal()}))
list.pl.allshrub.wc.f %>%
pull(plots)Upland Shrub Cover, Pooled Shrub Species - Combined Winter Range
#### summary tables
cov.allshrubs %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Upland cover - site type + time_class")| Upland cover - site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 13.32 | 8.76 | 2.00 | 11.00 | 31.00 | 21 | 100 |
| BL | non-core range | perc_cover | 21.21 | 15.03 | 1.33 | 19.50 | 47.33 | 20 | 100 |
| 2013 | core range | perc_cover | 10.75 | 8.11 | 0.33 | 8.83 | 24.33 | 20 | 100 |
| 2013 | non-core range | perc_cover | 24.32 | 19.75 | 1.67 | 19.00 | 66.33 | 22 | 100 |
| 2018 | core range | perc_cover | 13.61 | 8.04 | 1.67 | 12.67 | 34.33 | 18 | 100 |
| 2018 | non-core range | perc_cover | 20.06 | 15.11 | 0.33 | 18.67 | 56.67 | 22 | 100 |
cov.allshrubs %>%
filter(!is.na(range_type)) %>%
ggplot(aes(time_class, perc_cover)) +
geom_boxplot(fill="grey80") +
geom_point(color="#0095AF") +
theme_minimal() +
labs(x="Time class", y="% cover") +
facet_grid(burned~range_type)
# save plots
ggsave("./output/figures_202107_stats_upland/cover_range_type_boxplot.png", dpi = 300, width = 5.75, height = 3.75)cov.allshrubs %>%
filter(!is.na(valley_full)) %>%
ggplot(aes(time_class, perc_cover)) +
geom_boxplot(fill="grey80") +
geom_point(color="#0095AF") +
theme_minimal() +
labs(x="Time class", y="% cover") +
facet_wrap(~valley_full, ncol=3)
# save plots
ggsave("./output/figures_202107_stats_upland/cover_valley_boxplot.png", dpi = 300, width = 6.75, height = 6.25)Summary by Range Type
cov.allshrubs %>%
group_by(site_type, time_class, burned) %>%
descr(perc_cover, stats = "common") %>%
tb() %>%
mutate(across(where(is.numeric), round, 1)) %>%
datatable(caption = "Upland percent cover -- all shrub species pooled", filter = "bottom")Modeling and Posterior Description
## add cover in 0-1 range
cov.allshrubs <- cov.allshrubs %>%
mutate(cover = perc_cover/100)
stmod_cov.allshr <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.allshrubs,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.upl.allshrubcov <- insight::get_parameters(stmod_cov.allshr)
# tidy
posteriors.upl.allshrubcov.tidy <- posteriors.upl.allshrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.allshr)## Priors for model 'stmod_cov.allshr'
## ------
## 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.32])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.allshr)posteriors.upl.allshrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.upl.allshrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland, all shrub cover - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.allshr) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. all shrubs pooled") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs_ucunc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.allshr$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.65 | 0.9 | -1.92 | -1.37 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.02 | 0.9 | -0.26 | 0.29 | 0.55 | 0.43 | 0.05 |
| time_class2018 | -0.06 | 0.9 | -0.34 | 0.23 | 0.63 | 0.51 | 0.08 |
| (phi) | 10.25 | 0.9 | 7.53 | 13.24 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.18 | 0.9 | 0.16 | 0.21 | 1.00 | 1.00 | 0.00 |
| log-posterior | 47.51 | 0.9 | 34.00 | 60.04 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.allshr$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_allshrubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.allshr$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/PUTR_allshrubs_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.allshr) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_uplcov_ucunc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.uplcov.wcwnc.gt <- rope(stmod_cov.allshr$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.uplcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.50 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.47 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.uplcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_uplcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Pooled Shrub Species - Core Winter Range
#### summary tables
cov.allshrubs %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Upland cover - site type + time_class")| Upland cover - site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 13.32 | 8.76 | 2.00 | 11.00 | 31.00 | 21 | 100 |
| BL | non-core range | perc_cover | 21.21 | 15.03 | 1.33 | 19.50 | 47.33 | 20 | 100 |
| 2013 | core range | perc_cover | 10.75 | 8.11 | 0.33 | 8.83 | 24.33 | 20 | 100 |
| 2013 | non-core range | perc_cover | 24.32 | 19.75 | 1.67 | 19.00 | 66.33 | 22 | 100 |
| 2018 | core range | perc_cover | 13.61 | 8.04 | 1.67 | 12.67 | 34.33 | 18 | 100 |
| 2018 | non-core range | perc_cover | 20.06 | 15.11 | 0.33 | 18.67 | 56.67 | 22 | 100 |
Modeling and Posterior Description
cov.allshrubs.uc <- cov.allshrubs %>%
filter(site_type == "core range") # add
## add cover in 0-1 range
cov.allshrubs.uc <- cov.allshrubs.uc %>%
mutate(cover = perc_cover/100)
stmod_cov.allshr.uc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.allshrubs,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.upl.allshrubcov.uc <- insight::get_parameters(stmod_cov.allshr.uc)
# tidy
posteriors.upl.allshrubcov.tidy <- posteriors.upl.allshrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.allshr.uc)## Priors for model 'stmod_cov.allshr.uc'
## ------
## 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.32])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.allshr.uc)posteriors.upl.allshrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.upl.allshrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland, all shrub cover - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.allshr.uc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. all shrubs pooled") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs_uc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.allshr.uc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.65 | 0.9 | -1.92 | -1.37 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.02 | 0.9 | -0.26 | 0.29 | 0.55 | 0.43 | 0.05 |
| time_class2018 | -0.06 | 0.9 | -0.34 | 0.23 | 0.63 | 0.51 | 0.08 |
| (phi) | 10.25 | 0.9 | 7.53 | 13.24 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.18 | 0.9 | 0.16 | 0.21 | 1.00 | 1.00 | 0.00 |
| log-posterior | 47.51 | 0.9 | 34.00 | 60.04 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.allshr.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_allshrubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.allshr.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/PUTR_allshrubs_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.allshr.uc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_upl_uc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.uplcov.wcwnc.gt <- rope(stmod_cov.allshr.uc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.uplcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.50 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.47 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.uplcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_uplcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Pooled Shrub Species - Noncore Winter Range
Modeling and Posterior Description
cov.allshrubs.unc <- cov.allshrubs %>%
filter(site_type == "non-core range") # add
## add cover in 0-1 range
cov.allshrubs.unc <- cov.allshrubs.unc %>%
mutate(cover = perc_cover/100)
stmod_cov.allshr.unc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.allshrubs,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.upl.allshrubcov.unc <- insight::get_parameters(stmod_cov.allshr.unc)
# tidy
posteriors.upl.allshrubcov.tidy <- posteriors.upl.allshrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.allshr.unc)## Priors for model 'stmod_cov.allshr.unc'
## ------
## 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.32])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.allshr.unc)posteriors.upl.allshrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.upl.allshrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland, all shrub cover - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.allshr.unc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. all shrubs pooled") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs_unc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_stemcnt_upl_allshrubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.allshr.unc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.65 | 0.9 | -1.92 | -1.37 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.02 | 0.9 | -0.26 | 0.29 | 0.55 | 0.43 | 0.05 |
| time_class2018 | -0.06 | 0.9 | -0.34 | 0.23 | 0.63 | 0.51 | 0.08 |
| (phi) | 10.25 | 0.9 | 7.53 | 13.24 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.18 | 0.9 | 0.16 | 0.21 | 1.00 | 1.00 | 0.00 |
| log-posterior | 47.51 | 0.9 | 34.00 | 60.04 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.allshr.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_allshrubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.allshr.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/PUTR_allshrubs_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.allshr.unc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_uplcov_unc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.uplcov.wcwnc.gt <- rope(stmod_cov.allshr.unc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.uplcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.50 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.47 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.uplcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_uplcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
# Point estimate displayed: median
# HPD interval probability: 0.95
pl.emm.upl.allshr.ucunc <- emmeans(stmod_cov.allshr, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Combined Range", x = "Estimate", y = "Contrast")
pl.emm.upl.allshr.uc <- emmeans(stmod_cov.allshr.uc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Core Range", x = "Estimate", y = "Contrast")
pl.emm.upl.allshr.unc <- emmeans(stmod_cov.allshr.unc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Noncore Range", x = "Estimate", y = "Contrast")
pl.emm.upl.allshr.ucunc + pl.emm.upl.allshr.uc + pl.emm.upl.allshr.unc +
patchwork::plot_layout(ncol=3) +
plot_annotation('Upland Shrub Cover', caption = 'All Shrub Species Pooled')
Individual species
cov.byspp %>%
tabyl(shrub_species) %>%
arrange(-n) %>%
gt()| shrub_species | n | percent |
|---|---|---|
| PUTR | 79 | 0.26072607 |
| RICE | 52 | 0.17161716 |
| CHNA | 39 | 0.12871287 |
| ARTR | 35 | 0.11551155 |
| ROWO | 24 | 0.07920792 |
| JUCO | 19 | 0.06270627 |
| PRVI | 13 | 0.04290429 |
| CEFE | 8 | 0.02640264 |
| SYRO | 5 | 0.01650165 |
| ACGL | 3 | 0.00990099 |
| RIXX | 3 | 0.00990099 |
| RUDE | 3 | 0.00990099 |
| DAFR | 2 | 0.00660066 |
| JUSC | 2 | 0.00660066 |
| MARE | 2 | 0.00660066 |
| PHMO | 2 | 0.00660066 |
| SYOR | 2 | 0.00660066 |
| XXXX | 2 | 0.00660066 |
| ALIN | 1 | 0.00330033 |
| AMAL | 1 | 0.00330033 |
| CEXX | 1 | 0.00330033 |
| JAAM | 1 | 0.00330033 |
| PRXX | 1 | 0.00330033 |
| PSME | 1 | 0.00330033 |
| SHCA | 1 | 0.00330033 |
| UNKN | 1 | 0.00330033 |
#### summary
cov.byspp.summary <- cov.byspp %>%
group_by(time_class, site_type, shrub_species) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate(across(where(is.numeric), round, 1)) %>%
mutate(label = glue("{mean} ({sd})"))
# library(ggtext)
## ggtext for ggrich text
cov.byspp.summary %>%
# filter(site_type == "")
filter(n.valid>=3) %>%
filter(shrub_species != "SYRO" & shrub_species != "CEFE"& shrub_species != "PRVI"& shrub_species != "") %>%
ggplot(aes(time_class, shrub_species)) +
geom_tile(aes(fill = mean), color="grey80") +
geom_text(aes(label = label), size=.8) +
theme_minimal() +
scale_fill_gradientn(colors = c("#0095AF","#9ADCBB", "#FCFFDD")) +
facet_wrap(~site_type) +
theme(legend.position = "bottom") +
labs(x="", y="Species", title = "Winter Range", fill="Mean cover")
ggsave("./output/figures_202107_stats_upland/ucunc_tile_mean_sd_cov.png", dpi = 300, width = 6.75, height = 4.75)
cov.byspp.summary %>%
# filter(site_type == "")
filter(n.valid>=3) %>%
filter(shrub_species != "SYRO" & shrub_species != "CEFE"& shrub_species != "PRVI"& shrub_species != "") %>%
ggplot(aes(time_class, shrub_species)) +
geom_tile(aes(fill = mean), color="grey80") +
# geom_text(aes(label = label)) +
theme_minimal() +
scale_fill_gradientn(colors = c("#0095AF","#9ADCBB", "#FCFFDD")) +
facet_wrap(~site_type) +
theme(legend.position = "bottom") +
labs(x="", y="Species", title = "Winter Range", fill="Mean cover")
ggsave("./output/figures_202107_stats_upland/ucunc_tile_mean_sd_cov_nolbl.png", dpi = 300, width = 6.75, height = 4.75)cov.byspp.summary %>%
# filter(site_type == "")
filter(n.valid>=3) %>%
filter(shrub_species != "SYRO" & shrub_species != "CEFE"& shrub_species != "PRVI"& shrub_species != "") %>%
ggplot(aes(time_class, shrub_species)) +
geom_col(aes(fill = mean), color="grey80") +
# geom_text(aes(label = label)) +
theme_minimal() +
scale_fill_gradientn(colors = c("#0095AF","#9ADCBB", "#FCFFDD")) +
facet_wrap(~site_type) +
theme(legend.position = "bottom") +
labs(x="", y="Species", title = "Winter Range", fill="Mean cover")
ggsave("./output/figures_202107_stats_upland/ucunc_bar_mean_sd_cov_nolbl.png", dpi = 300, width = 6.75, height = 4.75)Upland Shrub Cover, Purshia tridentata - Combined Winter Range
# PUTR
cov.putr <- cov.byspp %>%
filter(shrub_species == "PUTR")
#### summary tables
cov.putr %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "PUTR cover ~ site type + time_class")| PUTR cover ~ site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 7.23 | 4.40 | 1.00 | 6.67 | 14.00 | 13 | 100 |
| BL | non-core range | perc_cover | 12.41 | 5.09 | 5.33 | 12.67 | 20.67 | 13 | 100 |
| 2013 | core range | perc_cover | 9.12 | 4.37 | 3.67 | 7.67 | 18.33 | 11 | 100 |
| 2013 | non-core range | perc_cover | 15.31 | 10.12 | 2.00 | 14.00 | 30.67 | 15 | 100 |
| 2018 | core range | perc_cover | 8.95 | 6.30 | 1.67 | 6.00 | 20.33 | 13 | 100 |
| 2018 | non-core range | perc_cover | 13.02 | 8.26 | 2.00 | 13.83 | 25.00 | 14 | 100 |
Modeling and Posterior Description
## add cover in 0-1 range
cov.putr <- cov.putr %>%
mutate(cover = perc_cover/100)
stmod_cov.putr <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.putr,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.putrubcov <- insight::get_parameters(stmod_cov.putr)
# tidy
posteriors.cov.putrubcov.tidy <- posteriors.cov.putrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.putr)## Priors for model 'stmod_cov.putr'
## ------
## 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.29,5.24])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.putr)posteriors.cov.putrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.putrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland, PUTR cover - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.putr) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. PUTR") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_PUTR_putrubs_unc.png", dpi = 300, width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.putr$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.98 | 0.9 | -2.27 | -1.69 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.16 | 0.9 | -0.21 | 0.54 | 0.76 | 0.69 | 0.27 |
| time_class2018 | 0.03 | 0.9 | -0.34 | 0.42 | 0.55 | 0.47 | 0.12 |
| (phi) | 9.33 | 0.9 | 6.63 | 12.00 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.13 | 0.9 | 0.10 | 0.16 | 1.00 | 1.00 | 0.00 |
| log-posterior | 26.11 | 0.9 | 17.15 | 35.34 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.putr$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_putr_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.putr$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description - PUTR", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/PUTR_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.putr) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_putr_cov_unc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.PUTRcov.wcwnc.gt <- rope(stmod_cov.putr$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.PUTRcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.30 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.36 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.PUTRcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_putr_cov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Purshia tridentata - Core Winter Range
#### summary tables
cov.putr %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Upland cover - site type + time_class")| Upland cover - site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 7.23 | 4.40 | 1.00 | 6.67 | 14.00 | 13 | 100 |
| BL | non-core range | perc_cover | 12.41 | 5.09 | 5.33 | 12.67 | 20.67 | 13 | 100 |
| 2013 | core range | perc_cover | 9.12 | 4.37 | 3.67 | 7.67 | 18.33 | 11 | 100 |
| 2013 | non-core range | perc_cover | 15.31 | 10.12 | 2.00 | 14.00 | 30.67 | 15 | 100 |
| 2018 | core range | perc_cover | 8.95 | 6.30 | 1.67 | 6.00 | 20.33 | 13 | 100 |
| 2018 | non-core range | perc_cover | 13.02 | 8.26 | 2.00 | 13.83 | 25.00 | 14 | 100 |
Modeling and Posterior Description
cov.putr.uc <- cov.putr %>%
filter(site_type == "core range") # add
## add cover in 0-1 range
cov.putr.uc <- cov.putr.uc %>%
mutate(cover = perc_cover/100)
stmod_cov.putr.uc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.putr,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.putrubcov.uc <- insight::get_parameters(stmod_cov.putr.uc)
# tidy
posteriors.cov.putrubcov.tidy <- posteriors.cov.putrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.putr.uc)## Priors for model 'stmod_cov.putr.uc'
## ------
## 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.29,5.24])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.putr.uc)posteriors.cov.putrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.putrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland,PUTR - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.putr.uc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. PUTR") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_PUTR_putrubs_uc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_PUTR_putrubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.putr.uc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.98 | 0.9 | -2.27 | -1.69 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.16 | 0.9 | -0.21 | 0.54 | 0.76 | 0.69 | 0.27 |
| time_class2018 | 0.03 | 0.9 | -0.34 | 0.42 | 0.55 | 0.47 | 0.12 |
| (phi) | 9.33 | 0.9 | 6.63 | 12.00 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.13 | 0.9 | 0.10 | 0.16 | 1.00 | 1.00 | 0.00 |
| log-posterior | 26.11 | 0.9 | 17.15 | 35.34 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.putr.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.putr.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/RICE_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.putr.uc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_upl_uc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.PUTRcov.wcwnc.gt <- rope(stmod_cov.putr.uc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.PUTRcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.30 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.36 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.PUTRcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_PUTRcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Purshia tridentata - Noncore Winter Range
Modeling and Posterior Description
cov.putr.unc <- cov.putr %>%
filter(site_type == "non-core range") # add
## add cover in 0-1 range
cov.putr.unc <- cov.putr.unc %>%
mutate(cover = perc_cover/100)
stmod_cov.putr.unc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.putr,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.putrubcov.unc <- insight::get_parameters(stmod_cov.putr.unc)
# tidy
posteriors.cov.putrubcov.tidy <- posteriors.cov.putrubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.putr.unc)## Priors for model 'stmod_cov.putr.unc'
## ------
## 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.29,5.24])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.putr.unc)posteriors.cov.putrubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.putrubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland,PUTR - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.putr.unc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. PUTR") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_PUTR_putrubs_unc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_PUTR_putrubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.putr.unc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -1.98 | 0.9 | -2.27 | -1.69 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.16 | 0.9 | -0.21 | 0.54 | 0.76 | 0.69 | 0.27 |
| time_class2018 | 0.03 | 0.9 | -0.34 | 0.42 | 0.55 | 0.47 | 0.12 |
| (phi) | 9.33 | 0.9 | 6.63 | 12.00 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.13 | 0.9 | 0.10 | 0.16 | 1.00 | 1.00 | 0.00 |
| log-posterior | 26.11 | 0.9 | 17.15 | 35.34 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.putr.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xRICE_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.putr.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/RICE_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.putr.unc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_PUTRcov_unc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.PUTRcov.wcwnc.gt <- rope(stmod_cov.putr.unc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.PUTRcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.30 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.36 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.PUTRcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_RICEcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Ribes cereum - Combined Winter Range
# PUTR
cov.rice <- cov.byspp %>%
filter(shrub_species == "RICE")
#### summary tables
cov.rice %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Upland cover ~ site type + time_class")| Upland cover ~ site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 7.19 | 6.28 | 0.33 | 10.00 | 15.67 | 9 | 100 |
| BL | non-core range | perc_cover | 3.96 | 2.97 | 1.00 | 3.50 | 8.67 | 8 | 100 |
| 2013 | core range | perc_cover | 3.67 | 3.94 | 0.33 | 2.00 | 13.33 | 9 | 100 |
| 2013 | non-core range | perc_cover | 6.03 | 6.85 | 0.33 | 3.00 | 20.00 | 10 | 100 |
| 2018 | core range | perc_cover | 2.38 | 3.54 | 0.33 | 1.33 | 11.00 | 8 | 100 |
| 2018 | non-core range | perc_cover | 3.67 | 4.52 | 0.33 | 1.33 | 12.33 | 8 | 100 |
Modeling and Posterior Description
## add cover in 0-1 range
cov.rice <- cov.rice %>%
mutate(cover = perc_cover/100)
stmod_cov.rice <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.rice,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.riceubcov <- insight::get_parameters(stmod_cov.rice)
# tidy
posteriors.cov.riceubcov.tidy <- posteriors.cov.riceubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.rice)## Priors for model 'stmod_cov.rice'
## ------
## 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.14,5.36])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.rice)posteriors.cov.riceubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.riceubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland,RICE - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.rice) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. RICE") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_PUTR_riceubs_unc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_PUTR_riceubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.rice$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -2.37 | 0.9 | -2.81 | -1.91 | 1.00 | 1.00 | 1.00 |
| time_class2013 | -0.06 | 0.9 | -0.60 | 0.47 | 0.57 | 0.51 | 0.23 |
| time_class2018 | -0.32 | 0.9 | -0.88 | 0.25 | 0.83 | 0.79 | 0.53 |
| (phi) | 6.35 | 0.9 | 4.02 | 8.95 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.08 | 0.9 | 0.05 | 0.11 | 1.00 | 0.96 | 0.00 |
| log-posterior | 41.59 | 0.9 | 34.56 | 48.62 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.rice$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xPUTR_riceubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.rice$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/rice_cov_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.rice) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_RICEcov_ucunc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.RICEcov.wcwnc.gt <- rope(stmod_cov.rice$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.RICEcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.27 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.17 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.90 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.RICEcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_RICEcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Ribes cereum - Core Winter Range
#### summary tables
cov.rice %>%
group_by(time_class, site_type) %>%
summarytools::descr(perc_cover, stats = "common") %>%
summarytools::tb() %>%
mutate_if(.predicate = is.numeric,.funs = round, digits = 2) %>%
gt() %>%
tab_header(title = "Upland cover - site type + time_class")| Upland cover - site type + time_class | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| time_class | site_type | variable | mean | sd | min | med | max | n.valid | pct.valid |
| BL | core range | perc_cover | 7.19 | 6.28 | 0.33 | 10.00 | 15.67 | 9 | 100 |
| BL | non-core range | perc_cover | 3.96 | 2.97 | 1.00 | 3.50 | 8.67 | 8 | 100 |
| 2013 | core range | perc_cover | 3.67 | 3.94 | 0.33 | 2.00 | 13.33 | 9 | 100 |
| 2013 | non-core range | perc_cover | 6.03 | 6.85 | 0.33 | 3.00 | 20.00 | 10 | 100 |
| 2018 | core range | perc_cover | 2.38 | 3.54 | 0.33 | 1.33 | 11.00 | 8 | 100 |
| 2018 | non-core range | perc_cover | 3.67 | 4.52 | 0.33 | 1.33 | 12.33 | 8 | 100 |
Modeling and Posterior Description
cov.rice.uc <- cov.rice %>%
filter(site_type == "core range") # add
## add cover in 0-1 range
cov.rice.uc <- cov.rice.uc %>%
mutate(cover = perc_cover/100)
stmod_cov.rice.uc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.rice,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.riceubcov.uc <- insight::get_parameters(stmod_cov.rice.uc)
# tidy
posteriors.cov.riceubcov.tidy <- posteriors.cov.riceubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.rice.uc)## Priors for model 'stmod_cov.rice.uc'
## ------
## 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.14,5.36])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.rice.uc)posteriors.cov.riceubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.riceubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland,RICE - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.rice.uc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. RICE") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_RICE_riceubs_uc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_RICE_riceubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.rice.uc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -2.37 | 0.9 | -2.81 | -1.91 | 1.00 | 1.00 | 1.00 |
| time_class2013 | -0.06 | 0.9 | -0.60 | 0.47 | 0.57 | 0.51 | 0.23 |
| time_class2018 | -0.32 | 0.9 | -0.88 | 0.25 | 0.83 | 0.79 | 0.53 |
| (phi) | 6.35 | 0.9 | 4.02 | 8.95 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.08 | 0.9 | 0.05 | 0.11 | 1.00 | 0.96 | 0.00 |
| log-posterior | 41.59 | 0.9 | 34.56 | 48.62 | 1.00 | 1.00 | 1.00 |
fun.sexit.gt(stmod_cov.rice.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/RICE_riceubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.rice.uc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/RICE_riceubs_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.rice.uc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_upl_uc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.RICEcov.wcwnc.gt <- rope(stmod_cov.rice.uc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.RICEcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.27 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.17 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.90 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.RICEcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_RICEcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Upland Shrub Cover, Ribes cereum - Noncore Winter Range
Modeling and Posterior Description
cov.rice.unc <- cov.rice %>%
filter(site_type == "non-core range") # add
## add cover in 0-1 range
cov.rice.unc <- cov.rice.unc %>%
mutate(cover = perc_cover/100)
stmod_cov.rice.unc <- stan_glmer(cover ~ time_class + (1 | site_id),
data = cov.rice.unc,
family = mgcv::betar,
prior_aux = exponential(2),
iter = 8000,
seed = 1234)# extract posteriors
posteriors.cov.riceubcov.unc <- insight::get_parameters(stmod_cov.rice.unc)
# tidy
posteriors.cov.riceubcov.tidy <- posteriors.cov.riceubcov %>%
pivot_longer(cols = starts_with("time"), names_to = "parameter", values_to = "values")Prior summary: cover ~ time_class + (1 | site_id)
prior_summary(stmod_cov.rice.unc)## Priors for model 'stmod_cov.rice.unc'
## ------
## 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.04,5.31])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsfun.model.summary(stmod_cov.rice.unc)posteriors.cov.riceubcov.tidy %>%
ggplot() +
geom_density(aes(fill = parameter, x = values),color="grey", alpha = .25) +
# scale_fill_viridis(discrete = TRUE) +
scale_fill_manual(values = colfunc2(ncol(posteriors.cov.riceubcov.tidy))) +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland,RICE - Combined Winter Range")
#### PP check plot
pp_check(stmod_cov.rice.unc) +
labs(x = "Value", y = "Density", caption = "upland cover ppc plot. RICE") +
theme_minimal()
# save plots
ggsave("./output/figures_202107_stats_upland/ppc_RICE_riceubs_unc.png", dpi = 300, width = 4.75, height = 3.75)
# ggsave("./output/figures_202107_stats_upland/ppc_RICE_riceubs.pdf", width = 4.75, height = 3.75)fun.sexit.gt(stmod_cov.rice.unc$stanfit)| Effect Existence and Significance Testing | |||||||
|---|---|---|---|---|---|---|---|
| Parameter | Median | CI | CI_low | CI_high | Direction | Significance | Large |
| (Intercept) | -2.11 | 0.9 | -2.85 | -1.42 | 1.00 | 1.00 | 1.00 |
| time_class2013 | 0.07 | 0.9 | -0.78 | 0.86 | 0.56 | 0.52 | 0.33 |
| time_class2018 | -0.18 | 0.9 | -1.03 | 0.71 | 0.63 | 0.60 | 0.41 |
| (phi) | 3.63 | 0.9 | 1.77 | 5.73 | 1.00 | 1.00 | 1.00 |
| mean_PPD | 0.11 | 0.9 | 0.05 | 0.18 | 1.00 | 0.97 | 0.00 |
| log-posterior | 10.00 | 0.9 | 4.27 | 14.94 | 0.99 | 0.99 | 0.99 |
fun.sexit.gt(stmod_cov.rice.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/xRICE_riceubs_stmod_cover1_posteriors.rtf")## save
fun.desc.post1.gt(stmod_cov.rice.unc$stanfit) %>%
tab_header(title = "Upland Shrub Cover Posterior Description", subtitle = "Combined Core & Noncore Winter Range") %>%
gt::gtsave(filename = "./output/tables/upland_cover/RICE_riceubs_stmod_cover1_posteriors.rtf")Effect Description
Probability of Direction
# Visualize the pd
#### PD - AC
p_direction(stmod_cov.rice.unc) %>% plot() +
theme_minimal() +
scale_fill_manual(values = colfunc3(2)) +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(caption = "Upland Shrub Cover Combined Winter Range") +
scale_fill_manual(values = colfunc3(2))
# save plots
ggsave("./output/figures_202107_stats_upland/pd_RICEcov_unc.png", dpi = 300, width = 4.75, height = 3.75)Region of Practical Equivalence (ROPE)
rope.RICEcov.wcwnc.gt <- rope(stmod_cov.rice.unc$stanfit, ci=0.9) %>%
gt() %>%
tab_header(title = "Percent in ROPE") %>%
fmt_number(
columns = 5,
decimals = 2,
suffixing = TRUE
)
rope.RICEcov.wcwnc.gt| Percent in ROPE | ||||||
|---|---|---|---|---|---|---|
| Parameter | CI | ROPE_low | ROPE_high | ROPE_Percentage | Effects | Component |
| (Intercept) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| time_class2013 | 0.9 | -0.1 | 0.1 | 0.18 | fixed | conditional |
| time_class2018 | 0.9 | -0.1 | 0.1 | 0.16 | fixed | conditional |
| (phi) | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
| mean_PPD | 0.9 | -0.1 | 0.1 | 0.43 | fixed | conditional |
| log-posterior | 0.9 | -0.1 | 0.1 | 0.00 | fixed | conditional |
# rope.RICEcov.wcwnc.gt %>%
# gt::gtsave(filename = "./output/tables/upland_cover/rope_RICEcov_wcwnc.rtf")Contrasts
EMM summaries with type = “response,” the tests and confidence intervals are done before back-transforming. The ratios estimated are ratios of geometric means. 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.
Additional plots
# putr
pl.emm.cov.putr.ucunc <- emmeans(stmod_cov.putr, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Combined Range", x = "", y = "Contrast")
pl.emm.cov.putr.uc <- emmeans(stmod_cov.putr.uc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Core Range", x = "Estimate", y = "")
pl.emm.cov.putr.unc <- emmeans(stmod_cov.putr.unc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Noncore Range", x = "", y = "")
pl.patch.putr <- pl.emm.cov.putr.ucunc + pl.emm.cov.putr.uc + pl.emm.cov.putr.unc +
patchwork::plot_layout(ncol=3) +
plot_annotation(title = expression(paste(italic("Pushia tridentata"))))
pl.patch.putr
ggsave("./output/figures_202107_stats_upland/uplcov_3panel_putr.png", dpi = 300, width = 7.75, height = 2.75)Ribes cereum
# Point estimate displayed: median
# HPD interval probability: 0.95
pl.emm.cov.rice.ucunc <- emmeans(stmod_cov.rice, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Combined Range", x = "", y = "Contrast")
pl.emm.cov.rice.uc <- emmeans(stmod_cov.rice.uc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Core Range", x = "Estimate", y = "")
pl.emm.cov.rice.unc <- emmeans(stmod_cov.rice.unc, ~ time_class) %>%
pairs() %>%
plot() +
theme_minimal() +
geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
labs(subtitle = "Noncore Range", x = "", y = "")
pl.patch.rice <- pl.emm.cov.rice.ucunc + pl.emm.cov.rice.uc + pl.emm.cov.rice.unc +
patchwork::plot_layout(ncol=3) +
plot_annotation(title = expression(paste(italic("Ribes cereum"))))
pl.patch.rice
# pl.patch.putr + pl.patch.rice + plot_layout(ncol=1)
ggsave("./output/figures_202107_stats_upland/uplcov_3panel_rice.png", dpi = 300, width = 7.75, height = 3.75)# rice
# pl.emm.cov.rice.ucunc <- emmeans(stmod_cov.rice, ~ time_class) %>%
# pairs() %>%
# plot() +
# theme_minimal() +
# geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
# labs(title= "Combined Range", subtitle = expression(paste(italic("Ribes cereum"))), x = "", y = "Contrast")
#
# pl.emm.cov.rice.uc <- emmeans(stmod_cov.rice.uc, ~ time_class) %>%
# pairs() %>%
# plot() +
# theme_minimal() +
# geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
# labs(title= "Core Range", subtitle = expression(paste(italic("Ribes cereum"))), x = "Estimate", y = "")
#
# pl.emm.cov.rice.unc <- emmeans(stmod_cov.rice.unc, ~ time_class) %>%
# pairs() %>%
# plot() +
# theme_minimal() +
# geom_vline(aes(xintercept = 0), color= "black", size=1, lty="dashed") +
# labs(title= "Noncore Range", subtitle = expression(paste(italic("Ribes cereum"))), x = "", y = "")
#
#
# pl.emm.cov.rice.ucunc + pl.emm.cov.rice.uc + pl.emm.cov.rice.unc +
# patchwork::plot_layout(ncol=3) +
# plot_annotation(title = 'Upland Shrub Cover', subtitle = expression(paste(italic("RICE"))))
#
# ggsave("./output/figures_202107_stats_upland/uplcov_3panel_rice.png", dpi = 300, width = 7.75, height = 3.75)Combined range
(pl.emm.upl.allshr.ucunc + labs(x="", subtitle="All shrubs")) + (pl.emm.cov.putr.ucunc + labs(x="Estimate",subtitle="PUTR")) + (pl.emm.cov.rice.ucunc + labs(subtitle="RICE")) +
patchwork::plot_layout(ncol=3) +
plot_annotation(title = 'Combined Range')
ggsave("./output/figures_202107_stats_upland/uplcov_3panel_ucunc_allshr_putr_rice.png", dpi = 300, width = 7.75, height = 3.75)Core range
(pl.emm.upl.allshr.uc + labs(subtitle="All shrubs")) + (pl.emm.cov.putr.uc + labs(subtitle="PUTR"))+ (pl.emm.cov.rice.uc + labs(subtitle="RICE")) +
patchwork::plot_layout(ncol=3) +
plot_annotation(title = 'Core Range')
ggsave("./output/figures_202107_stats_upland/uplcov_3panel_uc_allshr_putr_rice.png", dpi = 300, width = 7.75, height = 1.75)Noncore range
(pl.emm.upl.allshr.unc + labs(subtitle="All shrubs")) + (pl.emm.cov.putr.unc + labs(subtitle="PUTR")) + (pl.emm.cov.rice.unc + labs(subtitle="RICE")) +
patchwork::plot_layout(ncol=3) +
plot_annotation(title = 'Noncore Range')
ggsave("./output/figures_202107_stats_upland/uplcov_3panel_unc_allshr_putr_rice.png", dpi = 300, width = 7.75, height = 3.75)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)
Last updated: 2021 August 28
| Attached packages and versions | |
|---|---|
| package | loadedversion |
| bayestestR | 0.10.5 |
| dataMaid | 1.4.0 |
| dplyr | 1.0.7 |
| DT | 0.18 |
| emmeans | 1.6.2-1 |
| forcats | 0.5.1 |
| fs | 1.5.0 |
| ggplot2 | 3.3.5 |
| ggstance | 0.3.5 |
| ggtext | 0.1.1 |
| glue | 1.4.2 |
| gt | 0.3.1 |
| janitor | 2.1.0 |
| knitr | 1.33 |
| lubridate | 1.7.10 |
| mapview | 2.10.0 |
| patchwork | 1.1.1 |
| purrr | 0.3.4 |
| Rcpp | 1.0.7 |
| readr | 2.0.1 |
| readxl | 1.3.1 |
| rlang | 0.4.11 |
| rstanarm | 2.21.1 |
| sf | 1.0-2 |
| stringr | 1.4.0 |
| summarytools | 1.0.0 |
| tibble | 3.1.3 |
| tidyr | 1.1.3 |
| tidyverse | 1.3.1 |
| viridis | 0.6.1 |
| viridisLite | 0.4.0 |
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.”
Lemoine, Nathan P. 2019. “Moving Beyond Noninformative Priors: Why and How to Choose Weakly Informative Priors in Bayesian Analyses.” Oikos 128 (7): 912–28. https://doi.org/10.1111/oik.05985.
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.”