DI: Diversity-Interactions Model Fitting

Description

This function will fit a wide range of Diversity-Interactions (DI) models, one at a time. It provides some assisted automated ways to fit DI models, and includes the flexibility to extend DI models in several directions.

Usage

DI(y, prop, DImodel, custom_formula, data,
   block, density, treat, ID, FG, extra_formula,
   estimate_theta = FALSE, theta = 1)

Value

A model object of class "glm" including the components detailed in glm, plus the following:

DIcall: the call of the DI function
DIinternalcall: an internal call made within the DI model fitting process

Arguments

The minimum required arguments to use DI are either:

Argument DImodel with data, y and prop, or
Argument custom_formula with data.

The DImodel argument allows fitting of DI models via a range of 'tag' options that determine the form of the species interactions terms (the tags, described below, are STR, ID, AV, FG, ADD and FULL) and extra terms can be added to the model using the extra_formula argument. Using the argument custom_formula requires full specification of the model to be fitted using standard lm or glm syntax.

y

The column name of the response vector, which must be in quotes, for example, y = "yield".

block

The name of the block variable (if present), which must be in quotes, for example, block = "block". If no blocking variable, omit this argument.

density

The name of the density variable (if present), which must be in quotes, for example, density = "density". If no density variable, omit this argument.

prop

A vector of s column names identifying the species proportions in each community in the dataset. For example, if the species proportions columns are labelled p1 to p4, then prop = c("p1","p2","p3","p4"). Alternatively, the column numbers can be specified, for example, prop = 4:7, where species proportions are in the 4th to 7th columns.

treat

The name of a column in the dataset containing the value of a treatment factor or covariate. The treatment name must be included in quotes, for example, treat = "nitrogen". If the treatment is a factor, the variable must already be specified as a factor prior to using DI.

When used in conjunction with DImodel, the treatment will be included in the model as an additive factor or covariate, for example, specifying treat = nitrogen, DImodel = ID will fit the model p1 + p2 + ... + ps + nitrogen. Additional treatments, or interactions between the treatment and other model terms can be included via the extra_formula argument.
The treat argument is defunct when using the custom_formula argument, and any treatment must be included directly in the custom_formula argument.

ID

This argument takes a text list (of length s) dsecirbing groupings for the identity effects of the species. For example, if there are four species and you wish to group the identity effects all four species into a single term: ID could be ID = c("ID1","ID1","ID1","ID1"), where "ID1" is the name of the ID group. Similarly if the we wish to have two identity effect groups where identity effect of species 1 and 3, and species 2 and 4 are grouped together: ID could be ID = c("ID1","ID2","ID1","ID2"), where "ID1" and "ID2" are the names of the ID groups. These ideas expand to any number of species and any number or combination of groups. Finally, the ID groups do not have to be named "ID1" and "ID2", the user can specify any name for the groups.

If the ID argument is not specified, each species will be assumed to have a separate identity effect.
Specify an grouping for the ID does not affect the interaction terms. The interactions are still calculated using the individual species proportions.
The ID argument is defunct when using the custom_formula argument, since species identity effects must be included directly in the custom_formula argument.

FG

If species are classified by g functional groups, this argument takes a text list (of length s) of the functional group to which each species belongs. For example, for four grassland species with two grasses and two legumes: FG could be FG = c("G","G","L","L"), where G stands for grass and L stands for legume.

The FG argument is required if DImodel = "FG" is specified.
The FG argument is defunct when using the custom_formula argument, since species interactions must be included directly in the custom_formula argument.

DImodel

This argument is chosen (over custom_formula) to fit an automated version of a DI model. The chosen tag determines the form of the species interactions to be included in the model. The tags (or options) are:

STR (no identity or interaction effects, only an intercept is fitted, plus the experiment structural variables block, density and treat, if specified).

Each of the following includes the species proportions as specified in prop, the interaction variables according to the tag, plus block, density and treat if specified.

ID (no interaction terms),
AV (a single average species pairwise interaction variable),
FG (functional group interaction variables, the FG argument must be specified to use this option),
ADD (the additive species interaction variables),
FULL (all pairwise interaction variables).

The DImodel tag should appear in quotes, for example, DImodel = "STR".

extra_formula

In conjunction with DImodel, additional terms can be added using extra_formula. A ~ must be included before specifying the terms. For example, if DImodel = "AV" has been specified, adding extra_formula = ~ I(AV**2) will add a quadratic average pairwise interaction variable or extra_formula = ~ treatment:AV will add an interaction between the average pairwise species interaction variable and the treatment. Any variable included directly in extra_formula must already be contained in the dataset (interaction variables can be created using the function DI_data, if required).

custom_formula

To specify your own DI model, write your own model formula using the custom_formula argument. The standard notation from lm and glm is used here, for example, custom_formula = yield ~ 0 + p1:treatment + p2:treatment + p3:treatment + p4:treatment will fit the DI model with the identity effects each crossed with treatment, where treatment is a variable in the dataset. The custom_formula argument is recommended when the DImodel and extra_formula arguments are not sufficient. Any variable included directly in custom_formula must be already contained in the dataset (interaction variables can be created using the function DI_data, if required).

data

Specify the dataset, for example, data = Switzerland. The dataset name should not appear in quotes.

estimate_theta

By default, theta (the power parameter on all \(pi*pj\) components of each interaction variable in the model) is set equal to one. Specify estimate_theta = TRUE to include the estimation of \(\theta\) in the specified model. The estimate_theta argument can only be used in conjunction with the DImodel argument; if the custom_formula is used, then theta estimation is not available and the default estimate_theta = FALSE must be used.

theta

Users may specify a value of theta different than 1 to fit the DI model. Note that if estimate_theta = TRUE, then theta will be estimated via maximum profile log-likelihood and the value specified for theta will be overridden.

Author

Rafael A. Moral, John Connolly and Caroline Brophy

Details

What are Diversity-Interactions models?

Diversity-Interactions (DI) models (Kirwan et al 2009) are a set of tools for analysing and interpreting data from experiments that explore the effects of species diversity on community-level responses. We recommend that users of the DImodels package read the short introduction to DI models (available at: DImodels). Further information on DI models is available in Kirwan et al 2009 and Connolly et al 2013.

Checks on data prior to using DI.

Before using DI, check that the species proportions for each row in your dataset sum to one. See the 'Examples' section for code to do this. An error message will be generated if the proportions don't sum to one.

How does the DI function work?

The DI function provides wide flexibility in the types of Diversity-Interactions (DI) models that can be fitted. There are two ways to fit models in DI: 1) using DImodel, possibly augmented by extra_formula, or 2) using custom_formula. Models are estimated using iteratively reweighted least squares, via the glm package, when the option estimate_theta = FALSE.

Consider the following DI model, for example (in R formula syntax): y ~ p1 + p2 + p3 + treatment + p1:p2 + p1:p3 + p2:p3 + p1:p2:treatment + p1:p3:treatment + p2:p3:treatment

This model can be fitted using DImodel and extra_formula: DI(y = "y", prop = c("p1", "p2", "p3"), treat = "nitrogen", DImodel = "FULL", extra_formula = ~ p1:p2:treatment + p1:p3:treatment + p2:p3:treatment, data = datasetname)

or, by specifying all of the terms in the model using custom_formula: DI(custom_formula = y ~ p1 + p2 + p3 + treatment + p1:p2 + p1:p3 + p2:p3 + p1:p2:treatment + p1:p3:treatment + p2:p3:treatment, data = datasetname)

We recommend to use DImodel where possible, to augment with extra_formula where required, and to only use custom_formula when DImodel plus extra_formula is insufficient.

Including theta in DI models

A non-linear parameter \(\theta\) can be included in DI models as a power on all \(pi*pj\) components of each pairwise interaction variable. For example (in R formula syntax): y ~ p1 + p2 + p3 + (p1:p2)^theta + (p1:p3)^theta + (p2:p3)^theta for the full pairwise interaction model. Including \(\theta\) alters the contribution of the interaction term to the response (Connolly et al 2013).

By default, the value of \(\theta\) is 1. By specifying estimate_theta = TRUE within DI, a value of \(\theta\) will be estimated using profile likelihood over the space \(\theta\) = 0.01 to 1.5. The option estimate_theta = TRUE can only be used with DImodel, it is not available when using custom_formula.

As a general guideline to testing if \(\theta\) is required, we recommend:

1) finding the best form of the species interaction terms assuming theta = 1, and then,

2) testing if theta differs from 1.

If no species interaction terms are needed, then there is no need to do any testing for theta.

References

Connolly J, T Bell, T Bolger, C Brophy, T Carnus, JA Finn, L Kirwan, F Isbell, J Levine, A Lüscher, V Picasso, C Roscher, MT Sebastia, M Suter and A Weigelt (2013) An improved model to predict the effects of changing biodiversity levels on ecosystem function. Journal of Ecology, 101, 344-355.

Kirwan L, J Connolly, JA Finn, C Brophy, A Lüscher, D Nyfeler and MT Sebastia (2009) Diversity-interaction modelling - estimating contributions of species identities and interactions to ecosystem function. Ecology, 90, 2032-2038.

Examples

Run this code


## Load the Switzerland data
data(Switzerland)

## Check that the proportions sum to 1 (required for DI models)
## p1 to p4 are in the 4th to 7th columns in Switzerland
Switzerlandsums <- rowSums(Switzerland[4:7])
summary(Switzerlandsums)

## Fit the a simple AV DImodel with theta = 1
mod <- DI(y = "yield", prop = c("p1", "p2", "p3", "p4"), 
          DImodel = "AV", data = Switzerland)
summary(mod)

## Fit the same model but group the 4 species identity effect into 2 groups
mod <- DI(y = "yield", prop = c("p1", "p2", "p3", "p4"), 
          ID = c("ID1", "ID1", "ID2", "ID2"), DImodel = "AV", 
          data = Switzerland)
summary(mod)

## Combine the four identity effects into a single term and estimate theta
mod <- DI(y = "yield", prop = c("p1", "p2", "p3", "p4"), 
          ID = c("ID1", "ID1", "ID1", "ID1"), DImodel = "AV", 
          estimate_theta = TRUE, data = Switzerland)
summary(mod)

## Fit the FG DImodel, with factors density and treatment, and with theta = 1
m1 <- DI(y = "yield", density = "density", prop = 4:7, treat = "nitrogen", 
         FG = c("G","G","L","L"), DImodel = "FG", data = Switzerland)
summary(m1)

## Fit the FG DImodel, with factors density and treatment, and a fixed value of theta not equal to 1
m2 <- DI(y = "yield", density = "density", prop = 4:7, treat = "nitrogen", 
         FG = c("G","G","L","L"), DImodel = "FG", data = Switzerland, 
         theta = 0.5)
summary(m2)

## Test if the identity effects interact with nitrogen (and main effect of nitrogen excluded)
m3 <- DI(y = "yield", density = "density", prop = 4:7, FG = c("G", "G", "L", "L"), DImodel = "FG", 
         extra_formula = ~ (p1 + p2 + p3 + p4):nitrogen, data = Switzerland)
summary(m3)

## Fit the average pairwise model and check for a quadratic term using extra_formula.
## Need to create AV variable to be included in extra_formula and put in new dataset Switzerland2.
  AV_variable <- DI_data(prop = c("p1","p2","p3","p4"), data = Switzerland, what = "AV")
  Switzerland2 <- data.frame(Switzerland, "AV" = AV_variable)
  m4 <- DI(y = "yield", density = "density", prop = 4:7, DImodel = "AV", 
           extra_formula = ~ I(AV**2), data = Switzerland2)
  summary(m4)

## Using the custom_formula option to fit some, but not all, of the FG interactions.
## Fit the FG DImodel using custom_formula, with factors density and treatment, and theta = 1.
## Need to create functional group interaction variables for inclusion in custom_formulaand put
##  in new dataset Switzerland3.
  FG_matrix <- DI_data(prop = 4:7, FG = c("G","G","L","L"), data = Switzerland, what = "FG")
  Switzerland3 <- data.frame(Switzerland, FG_matrix)
  m5 <- DI(y = "yield", prop = c("p1","p2","p3","p4"), 
           custom_formula = yield ~ 0 + p1 + p2 + p3 + p4 + bfg_G_L + wfg_G 
           + density + nitrogen,
           data = Switzerland3)
  summary(m5)

Run the code above in your browser using DataLab