The recommended workflow of DImulti() analyses"

In DImodelsMulti: Fit Multivariate Diversity-Interactions Models with Repeated Measures

```{css styling, echo=FALSE} span.R { font-family: Courier New; }

```r
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

options(crayon.enabled = TRUE)

ansi_aware_handler <- function(x, options)
{
  paste0(
    "<pre class=\"r-output\"><code>",
    fansi::sgr_to_html(x = x, warn = FALSE, term.cap = "256"),
    "</code></pre>"
  )
}

old_hooks <- fansi::set_knit_hooks(knitr::knit_hooks,
                                 which = c("output", "message", "error", "warning"))
knitr::knit_hooks$set(
  output = ansi_aware_handler,
  message = ansi_aware_handler,
  warning = ansi_aware_handler,
  error = ansi_aware_handler
)

library(DImodelsMulti)
library(dplyr)
library(nlme)
library(ggplot2)

This vignette aims to give a rough outline for the process to be followed when examining multivariate and/or repeated measures BEF relationship study data using the DImulti() function from the DImodelsMulti R package.

We will use the dataset simMVRM, which is included in the package, to illustrate the workflow of the package.

Examining the data

We always start by looking at our data using View() or head(), to make sure it is as expected and includes all required information.

head(simMVRM)

We can then look at some summarising metrics, such as the mean of each ecosystem function, separated by time as below, which can help us know what to expect in the analysis, in this case we see that ecosystem functioning improves at time point 2 for the average multifunctionality index, MFindex, but the function Y2 performs better at time 1.

simMVRM_group <-  dplyr::summarise(dplyr::group_by(simMVRM, time),
                                   Y1 = mean(Y1),
                                   Y2 = mean(Y2),
                                   Y3 = mean(Y3),
                                   MFindex = mean(Y1 + Y2 + Y3))
simMVRM_group

We can also produce plots to the same effect.

hist(simMVRM[which(simMVRM$time == 1), ]$Y1, main = "Y1 at time 1", xlab = "Y1")

Fitting the first model

We use the main function of the DImodelsMulti R package, DImulti(), to fit a repeated measures Diversity-Interactions model to the data.

It is recommended to begin with the simplest model options and increase complexity as model selection/comparisons permit. The simplest model available to be fit by the package is the structural model, "STR", which only includes an intercept for each time and ecosystem function. $$ y_{kmt} = \beta_{0kt} + \epsilon_{kmt} $$

modelSTR <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "STR", method = "ML")

As our dataset is multivariate and in a wide format, we specify our response columns through the parameter y, "NA" through the first index of the vector eco_func, and select the unstructured ("UN") autocorrelation structure for our ecosystem functions through the second index of eco_func. The data also contains multiple time points, so we pass the column holding the time point indicator through the parameter time, along with the chosen autocorrelation structure, in this case we use compound symmetry, "CS". To facilitate grouping the recorded responses, a column index for the unique identifier of the experimental units is passed through unit_IDs. We choose to use the maximum likelihood ("ML") estimation method as we intend to compare models with differing fixed effects.

We can use print() or summary() to view information on the model in our console.

print(modelSTR)

From this evaluation, we can see that each of our coefficients are significant at an alpha significance level of 0.05 ($\alpha = 0.05$). We can also see the variance covariance matrices for our repeated measures and multiple ecosystem functions, along with the combined matrix, from which we can see the estimated strength and direction of the covarying relationships between response types.

Fitting a Diversity-Interactions model

The next model to be fit is the simplest Diversity-Interactions (DI) model available in the package, the identity ("ID") structure. $$ y_{kmt} = \sum^{S}{i=1}{\beta{ikt} p_{im}} + \epsilon_{kmt} $$ The DI modelling framework assumes that the main driver behind changes in ecosystem functioning is the initial relative abundance/proportion of the species.To reflect this, the intercept from the structural model we fit previously is replaced by the initial proportion of the species in the study, each with their own $\beta$ coefficient. These fixed effects form a simplex space. The only code that we need to change for this model is the DImodel tag, from "STR" to "ID".

modelID <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "ID", method = "ML")

We can now view this model, as before. In this case, we simply extract the t-table from the model summary.

summary(modelID)$tTable

The models within the DI modelling framework can be easily compared as they are hierarchical in nature, i.e. they are nested within one another. We can compare nested models using a likelihood ratio test, through anova(), or information criteria, such as AIC() or BIC() (the correct versions of which, AICc() or BICc() should be used in the cased of models with a large number of terms).

{width=100%} Image Credit: Kirwan et al. 2009

We choose to compare these models using a likelihood ratio test, where the null hypothesis states that the likelihoods of the two models do not significantly differ from one another, i.e. the extra terms in our larger model are not worth the added complexity.

anova(modelSTR, modelID)

As the p-value of the test is less than our chosen $\alpha$ value of 0.05, we reject the null hypothesis and proceed with the more complex model. We can confirm this selection by choosing the model with the lower AIC and BIC values, which are also printed by the anova() call.

Fitting and comparing interactions structures

We continue to increase the complexity of the models by fitting the next interaction structure in the nesting series, the average interaction model, "AV", a reparameterisation of the evenness model, "E". Again, the only code change required is the change of the DImodel tag. We then compare the two models.

modelAV <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "AV", method = "ML")

coef(modelAV)
anova(modelID, modelAV)

Once again, as our p-value < 0.05, we select the more complex model and continue increasing complexity.

Estimating theta

Now that we are fitting interactions, we may elect to estimate the non-linear parameter $\theta$, see vignette onTheta for a more in depth look at this process. We include the parameter estimate_theta with argument TRUE to have DImulti() automatically fit and compare different values of $\theta$, or if a user has a priori information on the nature of the term, they may pass values through the parameter theta. We can then test the two models against one another using information criteria, to see if the theta values significantly differ from 1.

modelAV_theta <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "AV", method = "ML",
                    estimate_theta = TRUE)

thetaVals <- modelAV_theta$theta
thetaVals

AICc(modelAV) 
AICc(modelAV_theta)

In this case, the inclusion of theta does improve the model.

The next models in the series are not nested within one another, so either can be used next but they cannot be compared to one another using a likelihood ratio test. As this dataset was simulated without any functional groupings, i.e. any groupings would be arbitrary, we choose to fit the additive interaction structure, "ADD".

modelADD <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "ADD", method = "ML",
                    theta = thetaVals)

modelADD$coefficients
anova(modelAV_theta, modelADD)

In this case, the p-value is greater than our chosen alpha level of 0.05, therefore we fail to reject the null hypothesis that the added complexity of the model does not significantly increase the model likelihood, and proceed in our analysis with the simpler model, modelAV.

As we do not need to increase the complexity of the model through the interaction structures any further, we turn to other options. It is at this stage that one could introduce treatment effects or environmental variables. As the simMVRM dataset does not contain any such information, example code is included below, where model modelAV_treat1 includes an intercept for treat for each level of ecosystem function and time point, while modelAV_treat2 crosses treat with each other fixed effect term.

modelAV_treat1 <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "ADD", method = "ML",
                    theta = thetaVals, extra_fixed = ~treat)

modelAV_treat2 <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "ADD", method = "ML",
                    theta = thetaVals, extra_fixed = ~1:treat)

The model could also alternatively be simplified here, by adding ID grouping, i.e. introducing functional redundancy within species ID effects. As simMVRM is simulated with no ID grouping in mind, any groups selected would be arbitrary and thus not recommended, however, example code for doing so is included below. These groupings do not affect the interaction term(s).

modelAV_ID <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "AV", method = "ML",
                    theta = thetaVals, ID = c("Group1", "Group1", "Group2", "Group2"))

summary(modelAV_ID)$tTable
anova(modelAV_ID, modelAV)

As the p-value is lower than 0.05, we reject the null hypothesis and continue with the more complex model, without any ID groupings.

If any errors are encountered in the process of fitting these models, please see vignette commonErrors for a guide.

Final model

As we have concluded our model selection process, the final step before interpretation is to refit the chosen model design using the method "REML", for unbiased estimates.

modelFinal <- DImulti(y = c("Y1", "Y2", "Y3"), eco_func = c("NA", "UN"), time = c("time", "CS"),
                    unit_IDs = 1, prop = 2:5, data = simMVRM, DImodel = "AV", method = "REML",
                    theta = thetaVals)

summary(modelFinal)

Interpretation

The final and most important step in the workflow of any analysis is the interpretation of the final model.

We first look at the fixed effect coefficients, their significance and direction. For example, we can see that the identity effect of species p1 has a significant, negative effect on ecosystem function Y1 at time point 1, but a significant positive effect on ecosystem function Y3 at the same time point. We can also see that the interaction effect AV has a significant positive effect on all ecosystem functions at both time points. As this interaction term is maxmimised in an ecosystem design where all four species are represented evenly, it may be beneficial to include species p1 despite its negative effect on Y1, the presence of other species, which have a positive effect on this ecosystem function could also 'make up' for this fault of species p1 (assuming that our goal is the maximisation of all ecosystem functions).

summary(modelFinal)$tTable

Next we could examine the variance covariance matrices of the model. These can be extracted from the model by either using the function getVarCov() from the R package nlme or by using the $ operator (the components of a model in R can be viewed using View(model)).

nlme::getVarCov(modelFinal)

modelFinal$vcov

nlme::getVarCov(modelFinal)

An example of an interpretation that we can make from these variance covariance matrices, is that the ecosystem functions Y1 and Y3 are negatively correlated, therefore would likely result in trade-offs when trying to estimated one over the other.

As multivariate repeated measures DI models can contain a large number of fixed effects, they can be difficult to interpret from. Our suggestions to circumvent this issue is to first predict from the model for community designs of interest, using the predict() function, as seen below, see vignette prediction for further details on this.

comms <- simMVRM[c(1, 4, 7, 10, 21), ]
print(comms)


commPred <- predict(modelFinal, newdata = comms)
commPred

Then, with these predictions, we can employ summarising metrics, such as a multifunctionality index, to summarise our findings while maintaining the benefits of having modelled each ecosystem function separately but simultaneously, alla Suter et al 2021. We could also use these predictions to generate a plot/graph to display our findings, such as a histogram, grouped bar chart, or ternary diagram.

ggplot2::ggplot(commPred, ggplot2::aes(fill=Ytype, y=Yvalue, x=plot)) +
  ggplot2::geom_bar(position="dodge", stat="identity")

A new R package, DImodelsVis, is in development, which will make plotting DI models fit using the R packages DImodels or DImodelsMulti a simple process.

References

Suter, M., Huguenin-Elie, O. and Lüscher, A., 2021. Multispecies for multifunctions: combining four complementary species enhances multifunctionality of sown grassland. Scientific Reports, 11(1), p.3835.

Kirwan, L., Connolly, J., Finn, J.A., Brophy, C., Lüscher, A., Nyfeler, D. and Sebastià, M.T., 2009. Diversity-interaction modeling: estimating contributions of species identities and interactions to ecosystem function. Ecology, 90(8), pp.2032-2038.

DImodelsMulti documentation built on May 29, 2024, 2:15 a.m.