Code
::install_github('diazrenata/soar') remotes
This is a modified subset of Appendix S3 from the article “Maintenance of community function through compensation breaks down over time in a desert rodent community” by Renata Diaz and S. K. Morgan Ernest, now published in Ecology (Diaz and Ernest 2022).
Compensation refers to the degree to which the remaining species on kangaroo rat removal plots absorb resources made available via kangaroo rat removal (Figure 1). We fit a generalized least squares (of the form compensation ~ timeperiod; note that “timeperiod” is coded as “oera” throughout) using the gls
function from the R package nlme
(Pinheiro, Bates, and R Core Team 2023). Because values from monthly censuses within each time period are subject to temporal autocorrelation, we included a continuous autoregressive temporal autocorrelation structure of order 1 (using the CORCAR1
function). We compared this model to models fit without the autocorrelation structure and without the time period term using AIC. The model with both the time period term and the autocorrelation structure was the best-fitting model via AIC (Table 1), and we used this model to calculate estimates and contrasts using the package emmeans
(Lenth et al. 2023) (Table 2, Table 3).
The following code downloads the data and prepares the compensation
data frame for analysis.
If needed, install the soar
package:
::install_github('diazrenata/soar') remotes
Data can be downloaded directly from the Portal data repository:
<- get_plot_totals(currency = "biomass")
plotl
<- list_plot_types() %>% filter(plot_type == "EE") plot_types
For interpretability, translating the era and treatment “names” as RMD coded them for analysis to the corresponding dates:
<- data.frame(
oera_df oera = c("a_pre_pb", "b_pre_reorg", "c_post_reorg"),
`Timeperiod` = c("1988-1997", "1997-2010", "2010-2020")
)
<- data.frame(oplottype = c("CC", "EE"),
oplot_df `Treatment` = c("Control", "Exclosure"))
<- data.frame(
contrasts_df contrast = c(
"a_pre_pb - b_pre_reorg",
"a_pre_pb - c_post_reorg",
"b_pre_reorg - c_post_reorg"
),Comparison = c(
"1988-1997 - 1997-2010",
"1988-1997 - 2010-2020",
"1997-2010 - 2010-2020"
) )
Because there are 5 exclosure plots and 4 control plots in these data, we remove 1 exclosure plot to achieve a balanced design. From the 5 possible exclosures to remove, we randomly select 1 using the seed 1977 (the year the Portal Project was initiated).
<- plot_types %>%
plot_types filter(plot_type == "EE")
set.seed(1977)
<- sample(plot_types$plot, 1, F) # results in removing plot 19
remove_plot
<- plotl %>%
plotl filter(plot != remove_plot)
Finally, take treatment-level means and calculate the compensation variable:
# Treatment-level means:
<- plots_to_treatment_means(plotl)
treatl
# Format column types
<- treatl %>%
treatl mutate(censusdate = as.Date(censusdate),
oera = ordered(oera),
oplottype = ordered(oplottype))
<- get_compensation(treatl) compensation
The following code fits the GLS models:
<-
comp_mean_gls gls(smgran_comp ~ oera,
correlation = corCAR1(form = ~ period),
data = compensation)
<-
comp_mean_gls_notime gls(smgran_comp ~ 1,
correlation = corCAR1(form = ~ period),
data = compensation)
<-
comp_mean_gls_noautoc gls(smgran_comp ~ oera, data = compensation)
<- gls(smgran_comp ~ 1, data = compensation) comp_mean_null
Model comparison via AIC:
<- data.frame(
compensation_comparison `Model specification` = c(
"intercept + timeperiod + autocorrelation",
"intercept + autocorrelation",
"intercept + timeperiod",
"intercept"
),AIC = c(
AIC(comp_mean_gls),
AIC(comp_mean_gls_notime),
AIC(comp_mean_gls_noautoc),
AIC(comp_mean_null)
) )
Calculate estimates:
<- emmeans(comp_mean_gls, specs = ~ oera)
comp_mean_gls_emmeans
<- oera_df %>%
compensation_estimates left_join(as.data.frame(comp_mean_gls_emmeans)) %>%
select(-oera)
Calculate contrasts:
<-contrasts_df %>%
compensation_contrasts left_join(as.data.frame(pairs(comp_mean_gls_emmeans))) %>%
mutate(p.value = round(p.value, digits = 4)) %>%
select(-contrast)
::kable(compensation_estimates) knitr
::kable(compensation_contrasts) knitr