The effect of a dominant shrub on annuals plants in a desert system changes along a gradient of precipitation.
Alex Filazzola
May 2017
Abstract
An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.
Hypothesis
We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.
Methods
Load libraries, data, and functions
## load packages
library(tidyverse)
library(effects)
library(lme4)
library(lmerTest)
## Load functions
se <- function(x, ...) sd(x)/sqrt(length(x)) ## standard error
source("functions.r")
source("rsquared.r")
## inverse hyperbolic sine transformation for zero laden data that fits log transformations
##Zhang, M., Fortney, J. C., Tilford, J. M., & Rost, K. M. (2000). An application of the inverse hyperbolic sine transformation—a note. Health Services and Outcomes Research Methodology, 1(2), 165-171.
ihs <- function(x) {
y <- log(x + sqrt(x ^ 2 + 1))
return(y)
}
## Load Aridity gradient
source("continentality.r") ## arid
arid$Year <- as.factor(arid$Year)
## Inverse aridity so larger numbers are more arid, smaler as less arid
arid$aridity <- 11 - arid$aridity
## load data
nutrient <- read.csv("Data//ERG.soilnutrients.csv")
community <- read.csv("Data//ERG.communitydata.csv")
community[is.na(community)] <- 0
census <- read.csv("Data//ERG.phytometer.census.csv")
HOBOdata <- read.csv("Data//ERG.logger.data.csv")
sitevars <- read.csv("Data//ERG.shrub.csv")Weather for growing season
Climate patterns within study
season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)
### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days
season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)## Site frost.days
## 1 Barstow 24.725275
## 2 Cuyama 29.670330
## 3 HeartofMojave 25.824176
## 4 PanocheHills 10.439560
## 5 SheepholeValley 6.043956
## 6 Tecopa 23.626374
## 7 TejonRanch 50.549451season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)## Site frost.days
## 1 Barstow 17.679558
## 2 Cuyama 19.889503
## 3 HeartofMojave 16.574586
## 4 PanocheHills 14.917127
## 5 SheepholeValley 1.785714
## 6 Tecopa 25.966851
## 7 TejonRanch 35.911602## Both years had comparable number of frost days
## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}
season1.frost <- season1 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 10
## 2 Cuyama 10
## 3 HeartofMojave 12
## 4 PanocheHills 5
## 5 SheepholeValley 7
## 6 Tecopa 11
## 7 TejonRanch 14season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 7
## 4 PanocheHills 6
## 5 SheepholeValley 3
## 6 Tecopa 7
## 7 TejonRanch 13## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)## Site frost.days avg.min.temp
## 1 Barstow 12.396694 5.750413
## 2 Cuyama 11.570248 3.352893
## 3 HeartofMojave 16.528926 4.542149
## 4 PanocheHills 2.479339 6.777686
## 5 SheepholeValley 2.479339 9.394167
## 6 Tecopa 11.570248 6.342149
## 7 TejonRanch 33.884298 1.438017season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)## Site frost.days avg.min.temp
## 1 Barstow 11.666667 6.052500
## 2 Cuyama 16.666667 3.624167
## 3 HeartofMojave 8.333333 5.637838
## 4 PanocheHills 8.333333 6.065833
## 5 SheepholeValley 0.000000 9.386916
## 6 Tecopa 20.833333 5.938333
## 7 TejonRanch 28.333333 2.535833season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 4
## 4 PanocheHills 2
## 5 SheepholeValley 2
## 6 Tecopa 4
## 7 TejonRanch 10season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 3
## 4 PanocheHills 3
## 5 SheepholeValley 2
## 6 Tecopa 5
## 7 TejonRanch 10Ambient plant community response to aridity gradient and shrub facilitation
## set colours for plots
obcol <- c("#E69F00","#56B4E9") ## Orange and blue
bgcol <- c("#000000","#707070") ## black and grey
scol <- obcol[1]
ocol <- obcol[2]
### Community First
commArid <- merge(community, arid, by=c("Site","Year"))
commArid[,"Year"] <- as.factor(commArid$Year)
## biomass model
m1 <- lmer(log(Biomass) ~ aridity * Microsite + Year + (1|ID), data=subset(commArid, Biomass>0))
anova(m1, test="Chisq")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## aridity 333.16 333.16 1 216.22 497.4088 < 2.2e-16 ***
## Microsite 71.50 71.50 1 526.28 106.7552 < 2.2e-16 ***
## Year 98.68 98.68 1 681.55 147.3370 < 2.2e-16 ***
## aridity:Microsite 5.62 5.62 1 540.59 8.3933 0.003919 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(residuals(m1))##
## Shapiro-Wilk normality test
##
## data: residuals(m1)
## W = 0.98756, p-value = 7.375e-06r.squared(m1)## Class Family Link Marginal Conditional AIC
## 1 lmerModLmerTest gaussian identity 0.5724485 0.6573578 1914.487## calcualte partial coefficents and confidence interval
ee <- Effect(c("aridity", "Microsite"), m1, xlevels=list(aridity=0:11))
## Extract polynomials
## Plot Biomass
plot1 <- ggplot(data=as.data.frame(ee), aes(x=aridity, y=fit, color=Microsite))+
geom_jitter(data=subset(commArid, Biomass>0), aes(x=aridity, y=log(Biomass), color=Microsite), size=2, width = 0.2, alpha=1) + theme_Publication() + ylab("log-transformed biomass")+
geom_line(lwd=2) +
geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c(scol, ocol)) + scale_fill_manual(values=c(scol, ocol)) +
annotate("text", x=0,y=4, label="a", size=8) +
guides(color=guide_legend(override.aes=list(fill=NA))) + xlim(0,10.5)
## species richness model
m2 <- glmer.nb(Richness ~ poly(aridity,2) * Microsite + Year + (1|ID), data=commArid)
car::Anova(m2, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Richness
## Chisq Df Pr(>Chisq)
## (Intercept) 337.0107 1 < 2e-16 ***
## poly(aridity, 2) 378.6697 2 < 2e-16 ***
## Microsite 4.3008 1 0.03809 *
## Year 2.1686 1 0.14086
## poly(aridity, 2):Microsite 278.0963 2 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## calcualte partial coefficents and confidence interval
ee <- Effect(c("aridity", "Microsite"), m2, xlevels=list(aridity=0:11))
## Plot richness
plot2 <- ggplot(data=as.data.frame(ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() + ylab("species richness")+
geom_jitter(data=subset(commArid, Biomass>0), aes(x=aridity, y=Richness, color=Microsite), size=2, width = 0.2, alpha=1) +
geom_line(lwd=2) +
geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c(scol, ocol)) + scale_fill_manual(values=c(scol, ocol))+
annotate("text", x=0,y=6, label="b", size=8) +
guides(color=guide_legend(override.aes=list(fill=NA))) + xlim(0,10.5)
## native vs non-native
status <- read.csv("Data//ERG.specieslist.csv")
statusComm <- community %>% gather(Species.shorthand, abundance, 13:53)
statusLong <- merge(statusComm , status, by="Species.shorthand")
statusComm <- statusLong %>% group_by(ID, Year, Site, Microsite, status, Rep) %>% summarize(abd= sum(abundance)) %>% data.frame(.) ## plants per plot
statusComm <- merge(statusComm, arid, by=c("Site","Year"))
## native only model
m4 <- glmer.nb( abd ~ poly(aridity,2) * Microsite + Year + (1|ID), data=subset(statusComm, status=="native"), nAGQ=0)
car::Anova(m4, test="Chisq", type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: abd
## Chisq Df Pr(>Chisq)
## (Intercept) 6.4990 1 0.01079 *
## poly(aridity, 2) 71.3878 2 3.150e-16 ***
## Microsite 51.8258 1 6.065e-13 ***
## Year 6.5262 1 0.01063 *
## poly(aridity, 2):Microsite 174.0731 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MuMIn::r.squaredGLMM(m4)## R2m R2c
## delta 0.3226551 0.7416200
## lognormal 0.3503267 0.8052230
## trigamma 0.2768442 0.6363241## calcualte partial coefficents and confidence interval
ee <- Effect(c("aridity", "Microsite"), m4, xlevels=list(aridity=0:11)) %>% data.frame()
## abundance of natives
plot3 <- ggplot(data=data.frame(ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() + ylab("native plant abundance")+
geom_jitter(data=subset(statusComm, status=="native" & abd>0), aes(x=aridity, y=abd, color=Microsite), size=2, width = 0.2, alpha=1) +
geom_line(lwd=2) +
geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c(scol, ocol)) + scale_fill_manual(values=c(scol, ocol))+
scale_y_continuous(trans='log2', limits=c(0.2,140)) +annotate("text", x=0,y=65, label="c", size=8) +
guides(color=guide_legend(override.aes=list(fill=NA))) + xlim(0,10.5)
## non-native only model
m5 <- glmer.nb(abd ~ poly(aridity,2) * Microsite + Year + (1|ID), data=subset(statusComm, status=="non.native"), nAGQ=0 )
car::Anova(m5, test="Chisq", type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: abd
## Chisq Df Pr(>Chisq)
## poly(aridity, 2) 872.2272 2 < 2.2e-16 ***
## Microsite 26.2480 1 3.003e-07 ***
## Year 84.7387 1 < 2.2e-16 ***
## poly(aridity, 2):Microsite 7.7281 2 0.02098 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## calcualte partial coefficents and confidence interval
ee <- Effect(c("aridity","Microsite"), m5, xlevels=list(aridity=0:11))
## abundance of non-natives
plot4 <- ggplot(data=as.data.frame(ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() + ylab("non-native abundance")+
geom_jitter(data=subset(statusComm, status=="non.native" & abd>0), aes(x=aridity, y=abd, color=Microsite), size=2, width = 0.2, alpha=1) +
geom_line(lwd=2) +
geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c(scol, ocol)) + scale_fill_manual(values=c(scol, ocol))+
scale_y_continuous(trans='log2', limits=c(0.5,420)) +annotate("text", x=0,y=420, label="d", size=8) +
guides(color=guide_legend(override.aes=list(fill=NA))) + xlim(0,10.5)
require(gridExtra)
grid.arrange(plot1, plot2, plot3, plot4)
Phytometer response to aridity gradient, nutrient addition, and shrub facilitation
census <- read.csv("Data//ERG.phytometer.census.csv")
end <- subset(census, Census=="end")
## Relevant nutrient to have ambient as control
end$Nutrient <- relevel(end$Nutrient, "Ambient")
phyto <- merge(end, arid, by=c("Year","Site"))
phyto$Year <- as.factor(phyto$Year)
phyto$Phacelia[is.na(phyto$Phacelia)] <- 0
## Run as Hurdle models
## create occurrence columns
phyto[,"pha.occ"] <- ifelse(phyto$Phacelia>0, 1,0)
phyto[,"pla.occ"] <- ifelse(phyto$Plantago>0, 1,0)
phyto[,"sal.occ"] <- ifelse(phyto$Salvia>0, 1,0)
## drop panoche 2017 because too cold
phyto <- phyto %>% filter(Site!="PanocheHills" | Year!="2017")
## Compare global models
occLong <- gather(phyto, species, occ, 23:25)
globalOcc <- glmer(occ ~ species:poly(aridity,2) + Microsite * species + species:Nutrient+ as.factor(Year) + (1|ID), data=occLong, family="binomial" , nAGQ=0)
car::Anova(globalOcc, test="Chisq")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: occ
## Chisq Df Pr(>Chisq)
## Microsite 2.5333 1 0.1114636
## species 98.5490 2 < 2.2e-16 ***
## as.factor(Year) 15.0918 1 0.0001024 ***
## species:poly(aridity, 2) 130.5057 6 < 2.2e-16 ***
## species:Microsite 29.0826 2 4.839e-07 ***
## species:Nutrient 1.2971 3 0.7298133
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bioLong <- gather(phyto, species, bio, 13:15)
globalBio <- lmer(log(bio) ~ species:poly(aridity,2) + Microsite * species + species:Nutrient + as.factor(Year) + (1|ID), data=bioLong)
anova(globalBio, test="Chisq")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Microsite 2.440 2.440 1 596.13 1.8376 0.175741
## species 34.356 17.178 2 608.25 12.9366 3.147e-06 ***
## as.factor(Year) 131.155 131.155 1 692.88 98.7733 < 2.2e-16 ***
## species:poly(aridity, 2) 108.260 18.043 6 613.10 13.5885 1.639e-14 ***
## species:Microsite 8.509 4.255 2 580.30 3.2041 0.041314 *
## species:Nutrient 17.572 5.857 3 444.22 4.4112 0.004519 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Phacelia occurrence
m1.occ <- glmer(pha.occ ~ poly(aridity,2) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial" , nAGQ=0)
glmer(pha.occ ~ poly(aridity,2, raw=T) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial" , nAGQ=0) ## produce co-efficients## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: binomial ( logit )
## Formula: pha.occ ~ poly(aridity, 2, raw = T) * Microsite * Nutrient +
## Year + (1 | ID)
## Data: phyto
## AIC BIC logLik deviance df.resid
## 923.4394 988.6695 -447.7197 895.4394 766
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 1.192
## Number of obs: 780, groups: ID, 210
## Fixed Effects:
## (Intercept)
## -0.80082
## poly(aridity, 2, raw = T)1
## 0.14694
## poly(aridity, 2, raw = T)2
## -0.02443
## Micrositeshrub
## 1.79388
## NutrientAdded
## -0.94322
## Year2017
## -0.55462
## poly(aridity, 2, raw = T)1:Micrositeshrub
## -0.22836
## poly(aridity, 2, raw = T)2:Micrositeshrub
## 0.01199
## poly(aridity, 2, raw = T)1:NutrientAdded
## 0.53514
## poly(aridity, 2, raw = T)2:NutrientAdded
## -0.05144
## Micrositeshrub:NutrientAdded
## -0.08353
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## -0.22829
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## 0.02844car::Anova(m1.occ, test="Chisq", type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: pha.occ
## Chisq Df Pr(>Chisq)
## poly(aridity, 2) 20.6793 2 3.233e-05 ***
## Microsite 24.2749 1 8.352e-07 ***
## Nutrient 0.2708 1 0.602797
## Year 8.4866 1 0.003578 **
## poly(aridity, 2):Microsite 2.8242 2 0.243626
## poly(aridity, 2):Nutrient 2.5338 2 0.281697
## Microsite:Nutrient 0.2762 1 0.599170
## poly(aridity, 2):Microsite:Nutrient 0.8215 2 0.663156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Phacelia Biomass
m1.bio <- lmer(log(Phacelia.biomass) ~ poly(aridity,2) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, pha.occ==1))
lmer(log(Phacelia.biomass) ~ poly(aridity,2, raw=T) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, pha.occ==1)) ## produce co-efficients## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: log(Phacelia.biomass) ~ poly(aridity, 2, raw = T) * Microsite *
## Nutrient + as.factor(Year) + (1 | ID)
## Data: subset(phyto, pha.occ == 1)
## REML criterion at convergence: 779.7963
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 0.7009
## Residual 1.2310
## Number of obs: 214, groups: ID, 135
## Fixed Effects:
## (Intercept)
## -1.7827229
## poly(aridity, 2, raw = T)1
## 0.0410249
## poly(aridity, 2, raw = T)2
## -0.0172080
## Micrositeshrub
## 1.1401106
## NutrientAdded
## 0.0555715
## as.factor(Year)2017
## 0.9680020
## poly(aridity, 2, raw = T)1:Micrositeshrub
## -0.0972093
## poly(aridity, 2, raw = T)2:Micrositeshrub
## -0.0005373
## poly(aridity, 2, raw = T)1:NutrientAdded
## 0.3016420
## poly(aridity, 2, raw = T)2:NutrientAdded
## -0.0310988
## Micrositeshrub:NutrientAdded
## 0.1922660
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## -0.5069094
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## 0.0612628anova(m1.bio, test="Chisq")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## poly(aridity, 2) 30.584 15.292 2 185.92 10.0916
## Microsite 9.673 9.673 1 155.32 6.3833
## Nutrient 5.044 5.044 1 128.69 3.3289
## as.factor(Year) 35.208 35.208 1 162.81 23.2347
## poly(aridity, 2):Microsite 2.321 1.161 2 157.69 0.7660
## poly(aridity, 2):Nutrient 0.520 0.260 2 176.85 0.1715
## Microsite:Nutrient 0.122 0.122 1 154.04 0.0807
## poly(aridity, 2):Microsite:Nutrient 2.972 1.486 2 158.12 0.9805
## Pr(>F)
## poly(aridity, 2) 6.906e-05 ***
## Microsite 0.01252 *
## Nutrient 0.07039 .
## as.factor(Year) 3.266e-06 ***
## poly(aridity, 2):Microsite 0.46661
## poly(aridity, 2):Nutrient 0.84255
## Microsite:Nutrient 0.77669
## poly(aridity, 2):Microsite:Nutrient 0.37738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Plantago occurrence
m2.occ <- glmer(pla.occ ~ poly(aridity,2) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial", nAGQ=0)
glmer(pla.occ ~ poly(aridity,2, raw=T) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial", nAGQ=0) ## produce co-efficient## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: binomial ( logit )
## Formula: pla.occ ~ poly(aridity, 2, raw = T) * Microsite * Nutrient +
## Year + (1 | ID)
## Data: phyto
## AIC BIC logLik deviance df.resid
## 792.6133 857.8255 -382.3067 764.6133 765
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 0.8672
## Number of obs: 779, groups: ID, 210
## Fixed Effects:
## (Intercept)
## -0.56441
## poly(aridity, 2, raw = T)1
## -0.28307
## poly(aridity, 2, raw = T)2
## 0.02314
## Micrositeshrub
## -3.62791
## NutrientAdded
## -0.82673
## Year2017
## -0.01892
## poly(aridity, 2, raw = T)1:Micrositeshrub
## 1.07604
## poly(aridity, 2, raw = T)2:Micrositeshrub
## -0.07123
## poly(aridity, 2, raw = T)1:NutrientAdded
## 0.64101
## poly(aridity, 2, raw = T)2:NutrientAdded
## -0.06021
## Micrositeshrub:NutrientAdded
## -12.86898
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## 3.16481
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## -0.20027car::Anova(m2.occ, test="Chisq", type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: pla.occ
## Chisq Df Pr(>Chisq)
## (Intercept) 32.0971 1 1.467e-08 ***
## poly(aridity, 2) 1.5843 2 0.452876
## Microsite 0.9069 1 0.340940
## Nutrient 1.4245 1 0.232665
## Year 0.0083 1 0.927585
## poly(aridity, 2):Microsite 7.6005 2 0.022365 *
## poly(aridity, 2):Nutrient 4.4085 2 0.110333
## Microsite:Nutrient 6.9106 1 0.008569 **
## poly(aridity, 2):Microsite:Nutrient 3.9588 2 0.138154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Plantago Biomass
m2.bio <- lmer(log(Plantago.biomass) ~ poly(aridity,2) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, pla.occ==1))
lmer(log(Plantago.biomass) ~ poly(aridity,2, raw=T) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, pla.occ==1)) ## produce co-efficient## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: log(Plantago.biomass) ~ poly(aridity, 2, raw = T) * Microsite *
## Nutrient + as.factor(Year) + (1 | ID)
## Data: subset(phyto, pla.occ == 1)
## REML criterion at convergence: 412.2272
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 0.4899
## Residual 0.9537
## Number of obs: 131, groups: ID, 92
## Fixed Effects:
## (Intercept)
## -1.39024
## poly(aridity, 2, raw = T)1
## -0.56821
## poly(aridity, 2, raw = T)2
## 0.02616
## Micrositeshrub
## -0.07088
## NutrientAdded
## 0.12524
## as.factor(Year)2017
## 2.55458
## poly(aridity, 2, raw = T)1:Micrositeshrub
## -0.16182
## poly(aridity, 2, raw = T)2:Micrositeshrub
## 0.02239
## poly(aridity, 2, raw = T)1:NutrientAdded
## -0.15517
## poly(aridity, 2, raw = T)2:NutrientAdded
## 0.02833
## Micrositeshrub:NutrientAdded
## -7.14305
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## 2.06838
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## -0.15018anova(m2.bio, test="Chisq")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## poly(aridity, 2) 0.140 0.070 2 113.68 0.0770
## Microsite 0.903 0.903 1 112.17 0.9928
## Nutrient 0.012 0.012 1 117.58 0.0134
## as.factor(Year) 109.035 109.035 1 117.58 119.8769
## poly(aridity, 2):Microsite 0.634 0.317 2 103.30 0.3486
## poly(aridity, 2):Nutrient 1.163 0.582 2 116.69 0.6394
## Microsite:Nutrient 0.952 0.952 1 112.58 1.0471
## poly(aridity, 2):Microsite:Nutrient 1.755 0.878 2 103.15 0.9649
## Pr(>F)
## poly(aridity, 2) 0.9259
## Microsite 0.3212
## Nutrient 0.9081
## as.factor(Year) <2e-16 ***
## poly(aridity, 2):Microsite 0.7065
## poly(aridity, 2):Nutrient 0.5294
## Microsite:Nutrient 0.3084
## poly(aridity, 2):Microsite:Nutrient 0.3844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## calcualte partial coefficents and confidence interval
ee <- Effect(c("aridity","Microsite"), m2.occ, se=T, xlevels=data.frame(aridity=seq(0,11, by=0.2), Microsite = c(rep("shrub",56),rep("open",56)))) %>% data.frame()
## Salvia occurrence
m3.occ <- glmer(sal.occ ~ poly(aridity,2) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial", nAGQ=0)
glmer(sal.occ ~ poly(aridity,2, raw=T) * Microsite * Nutrient + Year + (1|ID), data=phyto, family="binomial", nAGQ=0)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: binomial ( logit )
## Formula: sal.occ ~ poly(aridity, 2, raw = T) * Microsite * Nutrient +
## Year + (1 | ID)
## Data: phyto
## AIC BIC logLik deviance df.resid
## 1047.1696 1112.3817 -509.5848 1019.1696 765
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 0.3558
## Number of obs: 779, groups: ID, 210
## Fixed Effects:
## (Intercept)
## -0.259346
## poly(aridity, 2, raw = T)1
## 0.449431
## poly(aridity, 2, raw = T)2
## -0.054111
## Micrositeshrub
## -0.127842
## NutrientAdded
## -0.361166
## Year2017
## -0.448877
## poly(aridity, 2, raw = T)1:Micrositeshrub
## 0.146263
## poly(aridity, 2, raw = T)2:Micrositeshrub
## -0.010086
## poly(aridity, 2, raw = T)1:NutrientAdded
## 0.177820
## poly(aridity, 2, raw = T)2:NutrientAdded
## -0.013361
## Micrositeshrub:NutrientAdded
## -0.003606
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## -0.266290
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## 0.022996car::Anova(m3.occ, test="Chisq", type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: sal.occ
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0193 1 0.889521
## poly(aridity, 2) 13.3878 2 0.001238 **
## Microsite 2.1402 1 0.143482
## Nutrient 0.2897 1 0.590416
## Year 7.6415 1 0.005704 **
## poly(aridity, 2):Microsite 0.4538 2 0.796993
## poly(aridity, 2):Nutrient 0.5158 2 0.772662
## Microsite:Nutrient 3.7155 1 0.053908 .
## poly(aridity, 2):Microsite:Nutrient 0.4687 2 0.791092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Salvia Biomass
m3.bio <- lmer(log(Salvia.biomass) ~ poly(aridity,2) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, sal.occ==1))
lmer(log(Salvia.biomass) ~ poly(aridity,2, raw=T) * Microsite * Nutrient + as.factor(Year) + (1|ID), data=subset(phyto, sal.occ==1))## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: log(Salvia.biomass) ~ poly(aridity, 2, raw = T) * Microsite *
## Nutrient + as.factor(Year) + (1 | ID)
## Data: subset(phyto, sal.occ == 1)
## REML criterion at convergence: 1106.583
## Random effects:
## Groups Name Std.Dev.
## ID (Intercept) 0.8122
## Residual 1.0765
## Number of obs: 318, groups: ID, 181
## Fixed Effects:
## (Intercept)
## -2.362177
## poly(aridity, 2, raw = T)1
## 0.417876
## poly(aridity, 2, raw = T)2
## -0.059955
## Micrositeshrub
## -0.299705
## NutrientAdded
## 0.306645
## as.factor(Year)2017
## 1.039144
## poly(aridity, 2, raw = T)1:Micrositeshrub
## 0.015854
## poly(aridity, 2, raw = T)2:Micrositeshrub
## 0.007856
## poly(aridity, 2, raw = T)1:NutrientAdded
## -0.451572
## poly(aridity, 2, raw = T)2:NutrientAdded
## 0.069380
## Micrositeshrub:NutrientAdded
## -0.187874
## poly(aridity, 2, raw = T)1:Micrositeshrub:NutrientAdded
## 0.448441
## poly(aridity, 2, raw = T)2:Micrositeshrub:NutrientAdded
## -0.062226anova(m3.bio, test="Chisq")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## poly(aridity, 2) 12.954 6.477 2 253.00 5.5891
## Microsite 0.165 0.165 1 178.35 0.1420
## Nutrient 8.690 8.690 1 152.79 7.4991
## as.factor(Year) 52.781 52.781 1 295.70 45.5470
## poly(aridity, 2):Microsite 1.492 0.746 2 198.88 0.6438
## poly(aridity, 2):Nutrient 8.905 4.453 2 239.21 3.8424
## Microsite:Nutrient 0.412 0.412 1 178.29 0.3559
## poly(aridity, 2):Microsite:Nutrient 5.127 2.563 2 198.97 2.2121
## Pr(>F)
## poly(aridity, 2) 0.004215 **
## Microsite 0.706758
## Nutrient 0.006907 **
## as.factor(Year) 7.863e-11 ***
## poly(aridity, 2):Microsite 0.526388
## poly(aridity, 2):Nutrient 0.022777 *
## Microsite:Nutrient 0.551565
## poly(aridity, 2):Microsite:Nutrient 0.112156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## plot the phytometers
pha.ee <- Effect(c("aridity","Microsite"), m1.occ, se=T, xlevels=data.frame(aridity=seq(0,11, by=0.2), Microsite = c(rep("shrub",56),rep("open",56)))) %>% data.frame()
pla.ee <- Effect(c("aridity","Microsite"), m2.occ, se=T, xlevels=data.frame(aridity=seq(0,11, by=0.2), Microsite = c(rep("shrub",56),rep("open",56)))) %>% data.frame()
sal.ee <- Effect(c("aridity","Microsite"), m3.occ, se=T, xlevels=data.frame(aridity=seq(0,11, by=0.2), Microsite = c(rep("shrub",56),rep("open",56)))) %>% data.frame()
## phytometer occurrence
plot1 <- ggplot(data=as.data.frame(pha.ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() +
ylab(expression(italic("P. tanacetifolia")*"occurrence"))+ xlab("") +
geom_line(lwd=2) + geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c("#E69F00", "#56B4E9")) + scale_fill_manual(values=c("#E69F00", "#56B4E9"))+ xlim(0,10.5)
plot2 <- ggplot(data=as.data.frame(pla.ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() +
ylab(expression(italic("P. insularis")*"occurrence"))+ xlab("") +
geom_line(lwd=2) + geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c("#E69F00", "#56B4E9")) + scale_fill_manual(values=c("#E69F00", "#56B4E9"))+ xlim(0,10.5)
plot3 <- ggplot(data=as.data.frame(sal.ee), aes(x=aridity, y=fit, color=Microsite)) + theme_Publication() +
ylab(expression(italic("S. columbariae")*"occurrence"))+
geom_line(lwd=2) + geom_ribbon(aes(ymin = lower, ymax=upper, fill=Microsite), alpha=0.3, color=NA) +
scale_color_manual(values=c("#E69F00", "#56B4E9")) + scale_fill_manual(values=c("#E69F00", "#56B4E9"))+ xlim(0,10.5)
gridExtra::grid.arrange(plot1, plot2, plot3, nrow=3)
Analyze nutrients along aridity gradient in shrub canopies
### abiotic characteristics
nutrient <- read.csv("Data//ERG.soilnutrients.csv")
nutArid <- merge(nutrient, subset(arid, Year==2017), by=c("Site"))
nutMean <- nutArid %>% group_by(aridity, microsite) %>% summarize(n=mean(N),p=mean(P),k=mean(K), n.se=se(N), p.se=se(P), k.se=se(K))
m1 <- lm(log(N) ~ microsite * aridity, data=nutArid)
anova(m1)## Analysis of Variance Table
##
## Response: log(N)
## Df Sum Sq Mean Sq F value Pr(>F)
## microsite 1 40.160 40.160 28.3397 1.312e-06 ***
## aridity 1 18.617 18.617 13.1376 0.0005637 ***
## microsite:aridity 1 0.290 0.290 0.2045 0.6526263
## Residuals 66 93.528 1.417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(m1$residuals)##
## Shapiro-Wilk normality test
##
## data: m1$residuals
## W = 0.98748, p-value = 0.7143r.squared(m1)## Class Family Link Marginal Conditional AIC
## 1 lm gaussian identity 0.3870823 NA 228.9347m2 <- lm(log(P) ~ microsite * poly(aridity,2), data=nutArid)
anova(m2)## Analysis of Variance Table
##
## Response: log(P)
## Df Sum Sq Mean Sq F value Pr(>F)
## microsite 1 10.5455 10.5455 26.802 2.434e-06 ***
## poly(aridity, 2) 2 13.6050 6.8025 17.289 9.929e-07 ***
## microsite:poly(aridity, 2) 2 2.8558 1.4279 3.629 0.03214 *
## Residuals 64 25.1819 0.3935
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(m2$residuals)##
## Shapiro-Wilk normality test
##
## data: m2$residuals
## W = 0.93387, p-value = 0.001134r.squared(m2)## Class Family Link Marginal Conditional AIC
## 1 lm gaussian identity 0.5174789 NA 141.0856m3 <- lm(log(K) ~ microsite * aridity, data=nutArid)
anova(m3)## Analysis of Variance Table
##
## Response: log(K)
## Df Sum Sq Mean Sq F value Pr(>F)
## microsite 1 4.4454 4.4454 22.7483 1.059e-05 ***
## aridity 1 11.3446 11.3446 58.0533 1.271e-10 ***
## microsite:aridity 1 0.0001 0.0001 0.0004 0.9833
## Residuals 66 12.8976 0.1954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(m3$residuals)##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.96769, p-value = 0.06742r.squared(m3)## Class Family Link Marginal Conditional AIC
## 1 lm gaussian identity 0.5504149 NA 90.24938plot1 <- ggplot(data=nutArid, aes(x=aridity, y= log(N))) + theme_Publication() +
ylab("Nitrogen")+
geom_jitter(size=2, width = 0.2, alpha=1, aes(color=microsite)) +
geom_smooth(method="lm", lwd=2, lty=1, color="black") +
scale_color_manual(values=c(scol, ocol)) +
guides(color=guide_legend(override.aes=list(fill=NA)))
plot2 <- ggplot(data=nutArid, aes(x=aridity, y= log(P), color=microsite, fill=microsite)) + theme_Publication() +
ylab("Phosphorus")+
geom_jitter(size=2, width = 0.2, alpha=1) +
geom_smooth(method="lm", lwd=2, formula = y ~ poly(x,2)) +
scale_color_manual(values=c(scol, ocol)) +
scale_fill_manual(values=c(scol, ocol)) +
guides(color=guide_legend(override.aes=list(fill=NA)))
plot3 <- ggplot(data=nutArid, aes(x=aridity, y= log(K))) + theme_Publication() +
ylab("Potassium")+
geom_jitter(size=2, width = 0.2, alpha=1, aes(color=microsite)) +
geom_smooth(method="lm", lwd=2, lty=1, color="black") +
scale_color_manual(values=c(scol, ocol)) +
guides(color=guide_legend(override.aes=list(fill=NA)))
require(gridExtra)
grid.arrange(plot1, plot2, plot3, nrow=3)