R: Diversity-Interactions Model Fitting

DI {DImodels}

R Documentation

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)

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.

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.

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

Author(s)

Rafael A. Moral, John Connolly and Caroline Brophy

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

 

## 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)

[Package DImodels version 1.3.2 Index]