Diversity Interaction Models Example
Hira Naveed, Laura Byrne, Rishabh Vishwakarma
Introduction to DI models
Diversity-Interactions (DI) models are a modeling framework that estimates the effects of species identity and diversity on ecosystem function. They Predict the diversity–function relationship across different types of communities. DI models separately estimate the contributions of different species when they interact, and are usually an analysis of simplex design.
A DI model has the following 3 main elements:
Fig1: DI Model Explained
How DI models are used in practice
The models are a family of models and not one single model can explain its full functionality. One example of a DI model is displayed below:
Fig2: DI Model Syntax
How biological assumptions are tested in DI models
The below flow diagram shows the hierarchy of the models, where: * The different models are in boxes (models 1 to 8) * The arrows indicate where a model is hierarchical to the one below * The letters show the alternative biological hypotheses A–G that represent different assumptions about the relationship between diversity and ecosystem function
Our Data: simData
This is the simulated dataset that will be used throughout this example
## Community Richness Treatment p1 p2 p3 p4 Response
## 1 1 A 1.0000000 0.0000000 0.0000000 0.0000000 8.238
## 1 1 B 1.0000000 0.0000000 0.0000000 0.0000000 14.156
## 2 1 A 0.0000000 1.0000000 0.0000000 0.0000000 10.305
## 2 1 B 0.0000000 1.0000000 0.0000000 0.0000000 12.926
## 3 1 A 0.0000000 0.0000000 1.0000000 0.0000000 5.321
## 3 1 B 0.0000000 0.0000000 1.0000000 0.0000000 5.070
## 4 1 A 0.0000000 0.0000000 0.0000000 1.0000000 11.104
## 4 1 B 0.0000000 0.0000000 0.0000000 1.0000000 12.364
## 1 1 A 1.0000000 0.0000000 0.0000000 0.0000000 8.533
## 1 1 B 1.0000000 0.0000000 0.0000000 0.0000000 10.771
## 2 1 A 0.0000000 1.0000000 0.0000000 0.0000000 7.059
## 2 1 B 0.0000000 1.0000000 0.0000000 0.0000000 10.385
## 3 1 A 0.0000000 0.0000000 1.0000000 0.0000000 4.869
## 3 1 B 0.0000000 0.0000000 1.0000000 0.0000000 5.300
## 4 1 A 0.0000000 0.0000000 0.0000000 1.0000000 11.547
## 4 1 B 0.0000000 0.0000000 0.0000000 1.0000000 13.380
## 1 1 A 1.0000000 0.0000000 0.0000000 0.0000000 7.333
## 1 1 B 1.0000000 0.0000000 0.0000000 0.0000000 14.294
## 2 1 A 0.0000000 1.0000000 0.0000000 0.0000000 9.999
## 2 1 B 0.0000000 1.0000000 0.0000000 0.0000000 11.052
## 3 1 A 0.0000000 0.0000000 1.0000000 0.0000000 6.199
## 3 1 B 0.0000000 0.0000000 1.0000000 0.0000000 5.153
## 4 1 A 0.0000000 0.0000000 0.0000000 1.0000000 9.810
## 4 1 B 0.0000000 0.0000000 0.0000000 1.0000000 12.337
## 5 2 A 0.5000000 0.5000000 0.0000000 0.0000000 9.520
## 5 2 B 0.5000000 0.5000000 0.0000000 0.0000000 15.110
## 6 2 A 0.5000000 0.0000000 0.5000000 0.0000000 12.430
## 6 2 B 0.5000000 0.0000000 0.5000000 0.0000000 13.073
## 7 2 A 0.5000000 0.0000000 0.0000000 0.5000000 11.838
## 7 2 B 0.5000000 0.0000000 0.0000000 0.5000000 19.896
## 8 2 A 0.0000000 0.5000000 0.5000000 0.0000000 8.889
## 8 2 B 0.0000000 0.5000000 0.5000000 0.0000000 13.837
## 9 2 A 0.0000000 0.5000000 0.0000000 0.5000000 13.222
## 9 2 B 0.0000000 0.5000000 0.0000000 0.5000000 15.243
## 10 2 A 0.0000000 0.0000000 0.5000000 0.5000000 7.944
## 10 2 B 0.0000000 0.0000000 0.5000000 0.5000000 10.219
## 11 2 A 0.7500000 0.2500000 0.0000000 0.0000000 9.578
## 11 2 B 0.7500000 0.2500000 0.0000000 0.0000000 14.049
## 12 2 A 0.7500000 0.0000000 0.2500000 0.0000000 10.223
## 12 2 B 0.7500000 0.0000000 0.2500000 0.0000000 14.882
## 13 2 A 0.7500000 0.0000000 0.0000000 0.2500000 9.399
## 13 2 B 0.7500000 0.0000000 0.0000000 0.2500000 14.326
## 14 2 A 0.0000000 0.7500000 0.2500000 0.0000000 10.298
## 14 2 B 0.0000000 0.7500000 0.2500000 0.0000000 14.202
## 15 2 A 0.0000000 0.7500000 0.0000000 0.2500000 13.218
## 15 2 B 0.0000000 0.7500000 0.0000000 0.2500000 15.794
## 16 2 A 0.0000000 0.0000000 0.7500000 0.2500000 7.097
## 16 2 B 0.0000000 0.0000000 0.7500000 0.2500000 6.768
## 17 2 A 0.2500000 0.7500000 0.0000000 0.0000000 10.562
## 17 2 B 0.2500000 0.7500000 0.0000000 0.0000000 13.537
## 18 2 A 0.2500000 0.0000000 0.7500000 0.0000000 8.267
## 18 2 B 0.2500000 0.0000000 0.7500000 0.0000000 9.246
## 19 2 A 0.2500000 0.0000000 0.0000000 0.7500000 10.753
## 19 2 B 0.2500000 0.0000000 0.0000000 0.7500000 15.096
## 20 2 A 0.0000000 0.2500000 0.7500000 0.0000000 8.090
## 20 2 B 0.0000000 0.2500000 0.7500000 0.0000000 12.049
## 21 2 A 0.0000000 0.2500000 0.0000000 0.7500000 12.540
## 21 2 B 0.0000000 0.2500000 0.0000000 0.7500000 13.277
## 22 2 A 0.0000000 0.0000000 0.2500000 0.7500000 7.732
## 22 2 B 0.0000000 0.0000000 0.2500000 0.7500000 11.748
## 23 3 A 0.3333333 0.3333333 0.3333333 0.0000000 10.773
## 23 3 B 0.3333333 0.3333333 0.3333333 0.0000000 17.253
## 24 3 A 0.3333333 0.3333333 0.0000000 0.3333333 14.617
## 24 3 B 0.3333333 0.3333333 0.0000000 0.3333333 15.906
## 25 3 A 0.3333333 0.0000000 0.3333333 0.3333333 10.549
## 25 3 B 0.3333333 0.0000000 0.3333333 0.3333333 15.615
## 26 3 A 0.0000000 0.3333333 0.3333333 0.3333333 11.830
## 26 3 B 0.0000000 0.3333333 0.3333333 0.3333333 13.671
## 23 3 A 0.3333333 0.3333333 0.3333333 0.0000000 10.015
## 23 3 B 0.3333333 0.3333333 0.3333333 0.0000000 13.024
## 24 3 A 0.3333333 0.3333333 0.0000000 0.3333333 12.008
## 24 3 B 0.3333333 0.3333333 0.0000000 0.3333333 15.743
## 25 3 A 0.3333333 0.0000000 0.3333333 0.3333333 11.188
## 25 3 B 0.3333333 0.0000000 0.3333333 0.3333333 13.901
## 26 3 A 0.0000000 0.3333333 0.3333333 0.3333333 10.448
## 26 3 B 0.0000000 0.3333333 0.3333333 0.3333333 13.247
## 27 4 A 0.2500000 0.2500000 0.2500000 0.2500000 12.514
## 27 4 B 0.2500000 0.2500000 0.2500000 0.2500000 14.258
## 27 4 A 0.2500000 0.2500000 0.2500000 0.2500000 11.337
## 27 4 B 0.2500000 0.2500000 0.2500000 0.2500000 14.447
## 27 4 A 0.2500000 0.2500000 0.2500000 0.2500000 12.433
## 27 4 B 0.2500000 0.2500000 0.2500000 0.2500000 15.643
## 27 4 A 0.2500000 0.2500000 0.2500000 0.2500000 13.168
## 27 4 B 0.2500000 0.2500000 0.2500000 0.2500000 14.190
## 28 4 A 0.7000000 0.1000000 0.1000000 0.1000000 10.282
## 28 4 B 0.7000000 0.1000000 0.1000000 0.1000000 13.791
## 29 4 A 0.1000000 0.7000000 0.1000000 0.1000000 13.382
## 29 4 B 0.1000000 0.7000000 0.1000000 0.1000000 14.862
## 30 4 A 0.1000000 0.1000000 0.7000000 0.1000000 8.632
## 30 4 B 0.1000000 0.1000000 0.7000000 0.1000000 10.353
## 31 4 A 0.1000000 0.1000000 0.1000000 0.7000000 12.636
## 31 4 B 0.1000000 0.1000000 0.1000000 0.7000000 14.606
## 28 4 A 0.7000000 0.1000000 0.1000000 0.1000000 8.604
## 28 4 B 0.7000000 0.1000000 0.1000000 0.1000000 13.118
## 29 4 A 0.1000000 0.7000000 0.1000000 0.1000000 13.250
## 29 4 B 0.1000000 0.7000000 0.1000000 0.1000000 13.144
## 30 4 A 0.1000000 0.1000000 0.7000000 0.1000000 7.779
## 30 4 B 0.1000000 0.1000000 0.7000000 0.1000000 8.779
## 31 4 A 0.1000000 0.1000000 0.1000000 0.7000000 12.122
## 31 4 B 0.1000000 0.1000000 0.1000000 0.7000000 13.478Explore the Data
A quick summary:
summary(simData)## Community Richness Treatment p1
## Min. : 1.00 Min. :1.0 Length:100 Min. :0.0000
## 1st Qu.: 5.00 1st Qu.:2.0 Class :character 1st Qu.:0.0000
## Median :17.50 Median :2.0 Mode :character Median :0.1000
## Mean :16.26 Mean :2.4 Mean :0.2500
## 3rd Qu.:26.00 3rd Qu.:3.0 3rd Qu.:0.3333
## Max. :31.00 Max. :4.0 Max. :1.0000
## p2 p3 p4 Response
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. : 4.869
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 9.752
## Median :0.1000 Median :0.1000 Median :0.1000 Median :12.028
## Mean :0.2500 Mean :0.2500 Mean :0.2500 Mean :11.581
## 3rd Qu.:0.3333 3rd Qu.:0.3333 3rd Qu.:0.3333 3rd Qu.:13.701
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :19.896
The highest yielding observation:
simData[which(simData$Response == max(simData$Response)),]## Community Richness Treatment p1 p2 p3 p4 Response
## 30 7 2 B 0.5 0 0 0.5 19.896The lowest yielding observation:
simData[which(simData$Response == min(simData$Response)),]## Community Richness Treatment p1 p2 p3 p4 Response
## 13 3 1 A 0 0 1 0 4.869
A histogram of Response:
hist(simData$Response, xlab="Yield", main="")
#a <- ggplot(data=simulated)+
# geom_point(aes(x=p1, y=Response, colour=Treatment))+
# scale_colour_manual(values = c('#228B22','#8A228A'))
#b <- ggplot(data=simulated)+
# geom_point(aes(x=p2, y=Response, colour=Treatment))+
# scale_colour_manual(values = c('#228B22','#8A228A'))
#c <- ggplot(data=simulated)+
# geom_point(aes(x=p3, y=Response, colour=Treatment))+
# scale_colour_manual(values = c('#228B22','#8A228A'))
#d <- ggplot(data=simulated)+
# geom_point(aes(x=p4, y=Response, colour=Treatment))+
# scale_colour_manual(values = c('#228B22','#8A228A'))
#ggarrange(a,b,c,d, nrow=2, ncol=2, common.legend=TRUE, legend='bottom')
Boxplots of yield depending on treatment, community type, and richness:
for(i in 1:100)
{
simData$CType[i] <- if(simData$Richness[i] == 1) 'Monoculture'
else 'Mixture'
}
par(mfrow=c(1,2))
boxplot(Response~CType, data=simData[which(simData$Treatment=='A'),], col="#66CC00", main="Treatment A", ylab="Yield", xlab="")
boxplot(Response~CType, data=simData[which(simData$Treatment=='B'),], col="#66CC00", main="Treatment B", ylab="Yield", xlab="")
boxplot(Response~Richness, data=simData[which(simData$Treatment=='A'),], col="#66CC00", main="Treatment A", ylab="Yield", xlab="Richness")
boxplot(Response~Richness, data=simData[which(simData$Treatment=='B'),], col="#66CC00", main="Treatment B", ylab="Yield", xlab="Richness")
We see that Mixtures perform better than Monocultures, and that Treatment B outperforms Treatment A in most cases.
Data Preparation
Define S as the number of species in the data:
S=4Now calculate E, which is the Evenness term. Species evenness refers to how close in numbers each species in an environment is:
simData$E = ((2*S) / (S-1)) *
(simData$p1 * simData$p2 + simData$p1 * simData$p3 +
simData$p1 * simData$p4 + simData$p2 * simData$p3 +
simData$p2 * simData$p4 + simData$p3 * simData$p4)
Calculate species specific interaction terms:
simData$p1INT = simData$p1 * (1-simData$p1)
simData$p2INT = simData$p2 * (1-simData$p2)
simData$p3INT = simData$p3 * (1-simData$p3)
simData$p4INT = simData$p4 * (1-simData$p4)
Finally calculate the functional group proportions:
simData$Grass = simData$p1 + simData$p2
simData$Legume = simData$p3 + simData$p4Declare factors/Code categorical variables to numeric
for(i in 1:100)
{
simData$N[i] <- if(simData$Treatment[i] == 'A') 0 else 1
}
simData$Treatment <- as.factor(simData$Treatment)
Building the Models
Null Model - Richness
modelNullRich <- glm(Response ~ Richness, data = simData)Null Model - Treatment
modelNullTreat <- glm(Response ~ Treatment, data = simData)
Identity Effects Model
modelID <- glm(Response ~ 0 + p1 + p2 + p3 + p4, data = simData)Identity Effects Model - Treatment
modelIDTreat <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1:Treatment + p2:Treatment + p3:Treatment + p4:Treatment, data = simData)
Pairwise Interaction Model
modelPair <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1:p2 + p1:p3 + p1:p4 + p2:p3 + p2:p4 + p3:p4, data = simData)Pairwise Interaction Model - Treatment
modelPairTreat <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1:p2 + p1:p3 + p1:p4 + p2:p3 + p2:p4 + p3:p4 +
p1:Treatment + p2:Treatment + p3:Treatment + p4:Treatment +
p1:p2:Treatment + p1:p2:Treatment + p1:p3:Treatment + p1:p4:Treatment +
p2:p3:Treatment + p2:p4:Treatment + p3:p4:Treatment, data = simData)
Evenness Model
modelEven <- glm(Response ~ 0 + p1 + p2 + p3 + p4 + E, data = simData)Evenness Model - Treatment
modelEvenTreat <- glm(Response ~ 0 + p1 + p2 + p3 + p4 + E +
p1:Treatment + p2:Treatment + p3:Treatment + p4:Treatment + E:Treatment,
data = simData)
Average Interaction Model
modelAvg <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1INT + p2INT + p3INT + p4INT, data = simData)Average Interaction Model - Treatment
modelAvgTreat <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1INT + p2INT + p3INT + p4INT +
p1:Treatment + p2:Treatment + p3:Treatment + p4:Treatment +
p1INT:Treatment + p2INT:Treatment + p3INT:Treatment + p4INT:Treatment,
data = simData)
Functional Group Interaction Model
modelFG <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1:p2 + p3:p4 + Grass:Legume, data = simData)Functional Group Interaction Model - Treatment
modelFGTreat <- glm(Response ~ 0 + p1 + p2 + p3 + p4 +
p1:p2 + p3:p4 + Grass:Legume +
p1:Treatment + p2:Treatment + p3:Treatment + p4:Treatment +
Grass:Treatment + Legume:Treatment + Grass:Legume:Treatment, data = simData)
Functional Redundancy within Grass Group Model
modelGrass <- glm(Response ~ 0 + Grass + p3 + p4 +
Grass:p3 + Grass:p4 + p3:p4, data = simData)Functional Redundancy within Legume Group Model
modelLegume <- glm(Response ~ 0 + p1 + p2 + Legume +
p1:p2 + Legume:p1 + Legume:p2, data = simData)
F-Tests
Fig 2: Hierarchy of models
H0: The additional terms are all equal to 0
HA: At least one of the additional terms does not equal to 0
Do identity effects differ among species?
anova(modelNullRich, modelID, test = "F")$"Pr(>F)"[2]## [1] 1.548395e-06Yes
Do identity effects differ among species when Treatment is involved?
anova(modelNullTreat, modelIDTreat, test = "F")$"Pr(>F)"[2]## [1] 1.781743e-11Also yes
Does including treatment improve the identity effects model?
anova(modelID, modelIDTreat, test = "F")$"Pr(>F)"[2]## [1] 1.706387e-09Yes
Do species interact to produce diversity effects?
anova(modelIDTreat, modelPairTreat, test = "F")$"Pr(>F)"[2]## [1] 7.024828e-17Yes
Is including treatment still a good idea?
anova(modelPair, modelPairTreat, test = "F")$"Pr(>F)"[2]## [1] 1.044425e-17Also yes
Is a single interaction term sufficient?
anova(modelEvenTreat, modelPairTreat, test = "F")$"Pr(>F)"[2]## [1] 1.25545e-05No
Are additive interaction contributions sufficient?
anova(modelAvgTreat, modelPairTreat, test = "F")$"Pr(>F)"[2]## [1] 3.522333e-07No
Can patterns in separate pairwise interactions be described by functional groupings?
anova(modelFGTreat, modelPairTreat, test = "F")$"Pr(>F)"[2]## [1] 0.0594457(In the case of functional groupings being used, we would also do the following tests)
Are the two grasses functionally redundant?
anova(modelGrass, modelPair, test = "F")$"Pr(>F)"[2]## [1] 0.5383073Yes
Are the two legumes functionally redundant?
anova(modelLegume, modelPair, test = "F")$"Pr(>F)"[2]## [1] 7.046768e-09No
Can also use AIC and other diagnostic measures to decide
AIC(modelPair)## [1] 431.2386AIC(modelPairTreat)## [1] 327.2577AIC(modelEvenTreat)## [1] 355.8766AIC(modelAvgTreat)## [1] 363.2128AIC(modelFGTreat)## [1] 329.3279
## We decide that the best model is The Full Pairwise Interaction model crossed with Treatment
Using autoDI() Instead
autoDI(y="Response", prop=c("p1", "p2", "p3", "p4"), treat = "Treatment", FG = c("G", "G", "L", "L"), data = simData, selection= "Ftest")## One or more rows have species proportions that sum to approximately 1, but not exactly 1. This is typically a rounding issue, and has been corrected internally prior to analysis.
##
## --------------------------------------------------------------------------------
## Step 1: Investigating the diversity effect
##
## Selection using F tests
## Description
## DI Model 1 Structural 'STR' DImodel with treatment
## DI Model 2 Species identity 'ID' DImodel with treatment
## DI Model 3 Average interactions 'AV' DImodel with treatment
## DI Model 4 Functional group effects 'FG' DImodel with treatment
## DI Model 5 Separate pairwise interactions 'FULL' DImodel with treatment
##
## DI_model treat estimate_theta Resid. Df Resid. SSq Resid. MSq
## DI Model 1 STR 'Treatment' FALSE 98 677.7935 6.9163
## DI Model 2 ID 'Treatment' FALSE 95 385.8030 4.0611
## DI Model 3 AV 'Treatment' FALSE 94 202.2147 2.1512
## DI Model 4 FG 'Treatment' FALSE 92 158.7482 1.7255
## DI Model 5 FULL 'Treatment' FALSE 89 157.1654 1.7659
## Df SSq F Pr(>F)
## DI Model 1
## DI Model 2 3 291.9905 55.1164 <0.0001
## DI Model 3 1 183.5883 103.9628 <0.0001
## DI Model 4 2 43.4666 12.3072 <0.0001
## DI Model 5 3 1.5828 0.2988 0.8262
##
## Selected model: Functional group effects 'FG' DImodel with treatment
## Formula: Response = p1 + p2 + p3 + p4 + FG_bfg_G_L + FG_wfg_G + FG_wfg_L + TreatmentA
##
## --------------------------------------------------------------------------------
## Step 2: Investigating the treatment effect
##
## Selection using F tests
## Description
## DI Model 1 Functional group effects 'FG' DImodel
## DI Model 2 Functional group effects 'FG' DImodel with treatment
##
## DI_model treat estimate_theta Resid. Df Resid. SSq Resid. MSq
## DI Model 1 FG none FALSE 93 352.1945 3.7870
## DI Model 2 FG 'Treatment' FALSE 92 158.7482 1.7255
## Df SSq F Pr(>F)
## DI Model 1
## DI Model 2 1 193.4464 112.1088 <0.0001
##
## Selected model: Functional group effects 'FG' DImodel with treatment
## Formula: Response = p1 + p2 + p3 + p4 + FG_bfg_G_L + FG_wfg_G + FG_wfg_L + TreatmentA
##
## --------------------------------------------------------------------------------
## Step 3: Investigating whether theta is equal to 1 or not
##
## Theta estimate: 1.095606
## Selection using F tests
## Description
## DI Model 1 Functional group effects 'FG' DImodel with treatment
## DI Model 2 Functional group effects 'FG' DImodel with treatment, estimating theta
##
## DI_model treat estimate_theta Resid. Df Resid. SSq Resid. MSq
## DI Model 1 FG 'Treatment' FALSE 92 158.7482 1.7255
## DI Model 2 FG 'Treatment' TRUE 91 157.0827 1.7262
## Df SSq F Pr(>F)
## DI Model 1
## DI Model 2 1 1.6654 0.9648 0.3286
##
## Selected model: Functional group effects 'FG' DImodel with treatment
##
## --------------------------------------------------------------------------------
## Step 4: Comparing the final selected model with the reference (community) model
## 'community' is a factor with 31 levels, one for each unique set of proportions.
##
## model Resid. Df Resid. SSq Resid. MSq Df SSq F Pr(>F)
## DI Model 1 Selected 92 158.7482 1.7255
## DI Model 2 Reference 68 124.5743 1.8320 24 34.1739 0.7773 0.7507
##
## --------------------------------------------------------------------------------
## autoDI is limited in terms of model selection. Exercise caution when choosing your final model.
## --------------------------------------------------------------------------------##
## Call: glm(formula = new_fmla, family = family, data = new_data)
##
## Coefficients:
## p1 p2 p3 p4 FG_bfg_G_L FG_wfg_G
## 11.697 12.067 6.767 12.783 14.419 6.581
## FG_wfg_L TreatmentA
## 1.089 -2.782
##
## Degrees of Freedom: 100 Total (i.e. Null); 92 Residual
## Null Deviance: 14280
## Residual Deviance: 158.7 AIC: 348References
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.
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.