R: The simulated multivariate "simMV" dataset

simMV {DImodelsMulti}

R Documentation

The simulated multivariate "simMV" dataset

Description

The simMV dataset was simulated from a multivariate DI model. It contains 192 plots comprising of six species that vary in proportions (p1 - p6). Each plot was replicated once for testing a two-level factor treat, included at levels 0 and 1, resulting in a total of 384 plots. There are four simulated responses (Y1 - Y4) recorded in a wide data format (one column per response). The data was simulated assuming that there were existing covariances between the responses, an additive treatment effect, and both species identity and species interaction effects were present.

Usage

data("simMV")

Format

A data frame with 384 observations on the following twelve variables.

plot: A numeric vector identifying each unique experimental unit
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
treat: A two-level factor indicating whether a treatment was applied (1) or not (0)
Y1: A numeric vector giving the simulated response for ecosystem function 1
Y2: A numeric vector giving the simulated response for ecosystem function 2
Y3: A numeric vector giving the simulated response for ecosystem function 3
Y4: A numeric vector giving the simulated response for ecosystem function 4

Details

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. The multivariate DI model was developed by Dooley et al., 2015.

Parameter values for the simulation

Multivariate DI models take the general form of:

{y}_{km} = {Identities}_{km} + {Interactions}_{km} + {Structures}_{k} + {\epsilon}_{km}

where y are the community-level responses, the Identities are the effects of species identities for each response 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 simMV was simulated with:

identity effects for the six species for response:
- Y1 = 6.9, -0.3, 6.6, 1.7, -0.8, 4.3
- Y2 = 4.9, 3.6, 4.4, 2.3, 4.3, 6.6
- Y3 = -0.3, 4.6, 1.2, 6.8, 1.4, 6.9
- Y4 = 4.1, 4.2, -0.5, 4.9, 6.7, -0.9, -1.0
an average pairwise interaction effect for response:
- Y1 = 1.8
- Y2 = 6.8
- Y3 = 1.4
- Y4 = 0.3
a teatment effect for response:
- Y1 = 3.5
- Y2 = -0.3
- Y3 = 5.5
- Y4 = -1.0
\epsilon assumed to have a multivaraite normal distribution with mean 0 and Sigma:

[1,] 3.87 -0.17 -0.23 0.31

[2,] -0.17 1.34 -0.11 0.49

[3,] -0.23 -0.11 2.95 -0.36

[4,] 0.31 0.49 -0.36 1.83

References

Dooley, A., Isbell, F., Kirwan, L., Connolly, J., Finn, J.A. and Brophy, C., 2015.
Testing the effects of diversity on ecosystem multifunctionality using a multivariate model.
Ecology Letters, 18(11), pp.1242-1251.

Finn, J.A., Kirwan, L., Connolly, J., Sebastia, M.T., Helgadottir, A., Baadshaug, O.H., Belanger, G., Black, A., Brophy, C., Collins, R.P. and Cop, J., 2013.
Ecosystem function enhanced by combining four functional types of plant species in intensively managed grassland mixtures: a 3-year continental-scale field experiment.
Journal of Applied Ecology, 50(2), pp.365-375 .

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

Examples

###################################################################################################
###################################################################################################



## Modelling Example

# For a more thorough example of the workflow of this package, please see vignette
# DImulti_workflow using the following code:

# vignette("DImulti_workflow")

# We use head() to view the dataset
head(simMV)


# We can use the function DImulti() to fit a multivariate DI model, with an intercept for "treat"
# for each ecosystem function. The dataset is wide, so the Y columns are all entered through 'y'
# and the first index of eco_func is "NA". We fit the average interaction structure and use "ML"
# so that we can perform model comparisons with varying fixed effects
MVmodel <- DImulti(y = paste0("Y", 1:4), eco_func = c("NA", "UN"), unit_IDs = 1,
                   prop = paste0("p", 1:6), data = simMV, DImodel = "AV", extra_fixed = ~ treat,
                   method = "ML")
print(MVmodel)

# We can adjust the previous model to now cross "treat" with each other ID effect. We also specify
# different values of theta for each ecosystem function and use the "REML" method to get unbiased
# estimates.

MVmodel_theta <- DImulti(y = paste0("Y", 1:4), eco_func = c("NA", "UN"), unit_IDs = 1,
                   prop = paste0("p", 1:6), data = simMV, DImodel = "AV", extra_fixed = ~ 1:treat,
                   theta = c(1, 0.5, 0.8, 0.6), method = "REML")

summary(MVmodel_theta)


# We can now use any S3 method compatible with a gls object, for example, predict()

predict(MVmodel_theta, newdata = simMV[which(simMV$plot == 1), ])


##################################################################################################
#
##################################################################################################
## Code to simulate data
#
#

set.seed(412)

props <- data.frame(plot = integer(),
                    p1 = integer(),
                    p2 = integer(),
                    p3 = integer(),
                    p4 = integer(),
                    p5 = integer(),
                    p6 = integer())

index <- 1 #row number

#Monocultures
for(i in 1:6) #6 species
{
  for(j in 1:2) #2 technical reps
  {
    props[index, i+1] <- 1
    index <- index + 1
  }
}


#Equal Mixtures
for(rich in sort(rep(2:6, 3))) #3 reps at each richness level
{
  sp <- sample(1:6, rich) #randomly pick species from pool

  for(j in 1:2) #2 technical reps
  {
    for(i in sp)
    {
      props[index, i+1] <- 1/rich #equal proportions
    }
    index <- index + 1
  }
}


#Unequal Mixtures
for(rich in sort(rep(c(2, 3, 4, 5, 6), 15))) #15 reps at each richness level
{
  sp <- sample(1:6, rich, replace = TRUE) #randomly pick species from pool

  for(j in 1:2) #2 technical reps
  {
    for(i in 1:6)
    {
      props[index, i+1] <- sum(sp==i)/rich #equal proportions
    }
    index <- index + 1
  }
}


props[is.na(props)] <- 0

mySimData <- props

mySimData$treat <- 0
mySimDataDupe <- mySimData
mySimDataDupe$treat <- 1
mySimData <- rbind(mySimData, mySimDataDupe)

mySimData$plot <- 1:nrow(mySimData)

mySimData$Y1 <- NA
mySimData$Y2 <- NA
mySimData$Y3 <- NA
mySimData$Y4 <- NA

ADDs <- DImodels::DI_data(prop=2:7, what=c("ADD"), data=mySimData)
mySimData <- cbind(mySimData, ADDs)

E_AV <- DImodels::DI_data(prop=2:7, what=c("E", "AV"), data=mySimData)
mySimData <- cbind(mySimData, E_AV)

n <- 4 #Number of Ys
p <- qr.Q(qr(matrix(stats::rnorm(n^2), n))) #Principal Components (make sure it's positive definite)
S <- crossprod(p, p*(n:1)) #Sigma
m <- stats::runif(n, -0.25, 1.5)

#runif(8, -1, 7) #decide on betas randomly

for(i in 1:nrow(mySimData))
{
  #Within subject error
  error <- MASS::mvrnorm(n=1, mu=m, Sigma=S)
  mySimData$Y1[i] <- 6.9*mySimData$p1[i] + -0.3*mySimData$p2[i] + 6.6*mySimData$p3[i] +
  1.7*mySimData$p4[i] + -0.8*mySimData$p5[i] + 4.3*mySimData$p6[i] + 3.5*mySimData$treat[i] +
  1.8*mySimData$AV[i] + error[1]
  mySimData$Y2[i] <- 4.9*mySimData$p1[i] + 3.6*mySimData$p2[i] + 4.4*mySimData$p3[i] +
  2.3*mySimData$p4[i] + 4.3*mySimData$p5[i] + 6.6*mySimData$p6[i] + -0.3*mySimData$treat[i] +
  6.8*mySimData$AV[i] + error[2]
  mySimData$Y3[i] <- -0.3*mySimData$p1[i] + 4.6*mySimData$p2[i] + 1.2*mySimData$p3[i] +
  6.8*mySimData$p4[i] + 1.4*mySimData$p5[i] + 6.9*mySimData$p6[i] + 5.5*mySimData$treat[i] +
  1.4*mySimData$AV[i] + error[3]
  mySimData$Y4[i] <- 4.1*mySimData$p1[i] + 4.2*mySimData$p2[i] + -0.5*mySimData$p3[i] +
  4.9*mySimData$p4[i] + 6.7*mySimData$p5[i] + -0.9*mySimData$p6[i] + -1.0*mySimData$treat[i] +
  0.3*mySimData$AV[i] + error[4]
}

mySimData$treat <- as.factor(mySimData$treat)

[Package DImodelsMulti version 1.1.1 Index]