sim4 | R Documentation |
The sim4
dataset was simulated. There is a covariate treatment and six species that vary in proportions (p1 - p6
). It is assumed that species 1 and 2 come from functional group 1, species 3 and 4 from functional group 2 and species 5 and 6 from functional group 3. The response was simulated assuming that there were species identity effects, separate pairwise interaction effects and a covariate effect.
data(sim4)
A data frame with 141 observations on the following nine variables:
richness
A numeric vector identifying the number of species in the initial composition.
treatment
A covariate taking values 50, 150 or 250.
p1
A numeric vector indicating the initial proportion of species 1.
p2
A numeric vector indicating the initial proportion of species 2.
p3
A numeric vector indicating the initial proportion of species 3.
p4
A numeric vector indicating the initial proportion of species 4.
p5
A numeric vector indicating the initial proportion of species 5.
p6
A numeric vector indicating the initial proportion of species 6.
response
A numeric vector giving the simulated response variable.
What are Diversity-Interactions (DI) 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 strongly recommend that users read the short introduction to Diversity-Interactions models (available at: DImodels
). Further information on Diversity-Interactions models is also available in Kirwan et al 2009 and Connolly et al 2013.
Parameter values for the simulation
DI models take the general form of:
y = Identities + Interactions + Structures + \epsilon
where y is a community-level response, the Identities are the effects of species identities and enter the model as individual species proportions at the beginning of the time period, the Interactions are the interactions among the species proportions, while Structures include other experimental structures such as blocks, treatments or density.
The dataset sim4
was simulated with:
identity effects for the six species with values = 25, 16, 18, 20, 10, 12
a covariate effect = 0.03
all 15 pairwise interaction effects with values: 30, 27, 20, 15, 10, 9, 14, 18, 36, 17, 26, 32, 9, 21, 16 (for pairs of species 1-2, 1-3, 1-4, 1-5, 1-6, 2-3, 2-4, ... , 5-6 respectively).
theta = 1 (where \theta
is a non-linear parameter included as a power on each pipj
product within interaction variables, see Connolly et al 2013 for details)
\epsilon
assumed normally distributed with mean 0 and standard deviation 2.
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.
####################################
## Code to simulate the sim4 dataset
## Simulate dataset sim4 with 6 species, three functional groups and three levels of a covariate
## The species 1-2 are FG1, species 3-4 are FG2 and species 5-6 are FG3.
## Assume ID effects and the full pairwise interaction model, with a covariate.
## Set up proportions
data("design_b")
sim4a <- design_b
# Replicate the design for three values of a covariate
sim4b <- sim4a[rep(seq_len(nrow(sim4a)), times = 3), ]
sim4c <- data.frame(treatment = rep(c(50, 150, 250), each = 47))
sim4 <- data.frame(richness = sim4b[,1], sim4c, sim4b[,2:7])
row.names(sim4) <- NULL
## To simulate the response, first create a matrix of predictors that includes p1-p6, the treatment
## and all pairwise interaction variables
X <- model.matrix(~ p1 + p2 + p3 + p4 + p5 + p6 + treatment + (p1 + p2 + p3 + p4 + p5 + p6)^2 -1,
data = sim4)
## Create a vector of 'known' parameter values for simulating the response.
## The first six are the p1-p6 parameters, and the second set of one is the treatment parameter
## and the third set of 15 are the interaction parameters.
sim4_coeff <- c(25,16,18,20,10,12, 0.03, 30,27,20,15,10,9,14,18,36,17,26,32,9,21,16)
## Create response and add normally distributed error
sim4$response <- as.numeric(X %*% sim4_coeff)
set.seed(34261)
r <- rnorm(n = 141, mean = 0, sd = 2)
sim4$response <- round(sim4$response + r, digits = 3)
###########################
## Analyse the sim4 dataset
## Load the sim4 data
data(sim4)
## View the first few entries
head(sim4)
## Explore the variables in sim4
str(sim4)
## Check characteristics of sim4
hist(sim4$response)
summary(sim4$response)
plot(sim4$richness, sim4$response)
plot(sim4$richness[sim4$treatment==50], sim4$response[sim4$treatment==50], ylim=c(0,40))
plot(sim4$richness[sim4$treatment==150], sim4$response[sim4$treatment==150], ylim=c(0,40))
plot(sim4$richness[sim4$treatment==250], sim4$response[sim4$treatment==250], ylim=c(0,40))
plot(sim4$p1, sim4$response)
plot(sim4$p2, sim4$response)
plot(sim4$p3, sim4$response)
plot(sim4$p4, sim4$response)
plot(sim4$p5, sim4$response)
plot(sim4$p6, sim4$response)
## What model fits best? Selection using F-test
auto1 <- autoDI(y = "response", prop = 3:8, treat = "treatment",
FG = c("FG1","FG1","FG2","FG2","FG3","FG3"), data = sim4, selection = "Ftest")
summary(auto1)
## Ignore functional groups (will replace FG model with ADD model in Step 1 selection)
auto2 <- autoDI(y = "response", prop = 3:8, treat = "treatment", data = sim4, selection = "Ftest")
summary(auto2)
## Fit the functional group model using DI and the FG tag
m1 <- DI(y = "response", prop = 3:8, treat = "treatment",
FG = c("FG1","FG1","FG2","FG2","FG3","FG3"), DImodel = "FG", data = sim4)
summary(m1)
## Fit the additive species model using DI and the ADD tag
m2 <- DI(y = "response", prop = 3:8, treat = "treatment", DImodel = "ADD", data = sim4)
summary(m2)
## Fit the full pairwise model using DI and the FULL tag
m3 <- DI(y = "response", prop = 3:8, treat = "treatment", DImodel = "FULL", data = sim4)
summary(m3)
plot(m3)
## Check goodness-of-fit using a half-normal plot with a simulated envelope
library(hnp)
hnp(m3)
## Create the functional group and additive species interaction variables,
## and store in a new data frame called sim4a
newlist <- DI_data(prop = 3:8, FG = c("FG1","FG1","FG2","FG2","FG3","FG3"),
data = sim4, what = c("FG", "ADD"))
sim4a <- data.frame(sim4, newlist$FG, newlist$ADD)
## Fit the functional group model using DI and custom_formula (equivalent to m1)
m4 <- DI(custom_formula = response ~ 0 + p1 + p2 + p3 + p4 + p5 + p6 + bfg_FG1_FG2
+ bfg_FG1_FG3 + bfg_FG2_FG3 + wfg_FG1 + wfg_FG2 + wfg_FG3 + treatment, data = sim4a)
summary(m4)
## Fit the additive species model using DI and custom_formula (equivalent to m2)
m5 <- DI(custom_formula = response ~ 0 + p1 + p2 + p3 + p4 + p5 + p6 + p1_add
+ p2_add + p3_add + p4_add + p5_add + p6_add + treatment, data = sim4a)
summary(m5)
## Fit the full pairwise model using DI and custom_formula (equivalent to m3)
m6 <- DI(custom_formula = response ~ 0 + p1 + p2 + p3 + p4 + p5 + p6
+ (p1 + p2 + p3 + p4 + p5 + p6)^2 + treatment, data = sim4a)
summary(m6)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.