In AnneLew/Lewerentz-etal_2021_Macrophytes-DDG: Analysis of the Depth Diversity Gradient of macrophytes
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
#Load packages
devtools::load_all(path = ".")
library(tidyverse)
library(knitr)
library(ggpubr)
library(PerformanceAnalytics)
library(data.table)
library(kableExtra)
library(knitr)
library(ggpmisc)
library(vegan)
library(gamm4)
#library(corrplot)
library(voxel) #GAMM plotting using ggplot2
library(multcomp) #Herberich Test
library(multcompView) #Herberich Test
library(sandwich) #Herberich Test
library(broom) #Herberich Test
#devtools::install_github("vqv/ggbiplot")
library(ggbiplot)
library(lubridate)
## ggplot theme updates
source(paste0(here::here(),"/R/theme.R"))
General information
Lake morphology
Table 1: Morphology of selected lakes: Lake surface area [ha] (Area_ha) and depths of deepest point of the lake [m] (maxDepth_m)
lakes<-unique(Makroph_Lake_DepthS$Lake)
Morpho<-Morphology[Morphology$Name_Makro_short %in% lakes, ]
Morpho %>% dplyr::select(Lake, Area_ha, maxDepth_m)
Depth diversity gradients (DDG) of macrophytes: shape
Correlations between richness components
richness<-Makroph_Lake_DepthS[c(97,101,99)]
colnames(richness) <- c("alpha richness", "beta richness", "gamma richness")
chart.Correlation(richness, histogram=F, pch=9, method = "pearson")
Figure 1: Correlations between diversity metrices (method = pearson). Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
DDG of all field campaigns
A1<-ggplot(data=Makroph_Lake_DepthS)+
geom_line(aes(x=Tiefe, y=ALPHA, group=interaction(Lake, YEAR)),col="grey")+
scale_x_reverse()+
geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mAlpha))+
geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mAlpha-sdAlpha, ymax=mAlpha+sdAlpha), width=.1)+
ylab("Alpha richness")+xlab(("Depth (m)"))
B1<-ggplot(data=Makroph_Lake_DepthS)+
geom_line(aes(x=Tiefe, y=BETA, group=interaction(Lake, YEAR)),col="grey")+ scale_x_reverse()+
geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mBeta))+
geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mBeta-sdBeta, ymax=mBeta+sdBeta), width=.1)+
ylab("Beta richness")+xlab(("Depth (m)"))
C1<-ggplot(data=Makroph_Lake_DepthS)+
geom_line(aes(x=Tiefe, y=GAMMA, group=interaction(Lake, YEAR)),col="grey")+ scale_x_reverse()+
geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mGamma))+
geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mGamma-sdGamma, ymax=mGamma+sdGamma), width=.1)+
ylab("Gamma richness")+xlab(("Depth (m)"))
ggarrange(A1,B1,C1,ncol=3,nrow=1,labels=c("(a)","(b)","(c)"))
Figure 2: Depth diversity gradients of macrophytes for alpha (a), beta (b) and gamma richness (c) with mean as black line and sd as black bars. Each single grey line is one field campaign (lake*year).
DDG per lake
### ALPHA & Peak Plot
ggplot(data=Makroph_Lake_DepthS, aes(x=(Tiefe), y=ALPHA, col=factor(YEAR),
group=interaction(Lake,YEAR)))+
geom_line(aes(linetype=datasettotsimpl))+
scale_linetype_manual(values=c("dotted", "solid"))+
facet_wrap(~ Lake, ncol=5)+
ylab("Alpha richness")+ xlab("")+labs(fill = "Pattern type")+ylim(0,15)+
geom_errorbar(data=Makroph_Lake_DepthS,aes(ymin=ALPHA-ALPHAsd, ymax=ALPHA+ALPHAsd), width=.1)+
geom_point(data=PEAK, aes(x=(AlphaPeakDepth), y=AlphaPeakRichness, shape=datasettotsimpl))+
scale_x_reverse()+
labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+
ggtitle("(a) Alpha richness")
### BETA & Peak Plot
ggplot(Makroph_Lake_DepthS, aes(x=Tiefe, y=BETA, group=interaction(Lake,YEAR), col=factor(YEAR),
linetype=datasettotsimpl))+
geom_line()+facet_wrap(~ Lake, ncol=5)+
scale_linetype_manual(values=c("dotted", "solid"))+
ylab("Beta richness")+ xlab("Depth (m)")+labs(fill = "Pattern type")+
geom_point(data=PEAK, aes(x=BetaPeakDepth, y=BetaPeakRichness, shape=datasettotsimpl))+
scale_x_reverse()+
labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+
ggtitle("(b) Beta richness")
### GAMMA & Peak Plot
ggplot(Makroph_Lake_DepthS, aes(x=Tiefe, y=GAMMA, group=interaction(Lake,YEAR), col=factor(YEAR),
linetype=datasettotsimpl))+
geom_line()+facet_wrap(~ Lake, ncol=5)+
scale_linetype_manual(values=c("dotted", "solid"))+
ylab("Gamma richness")+ xlab("Depth (m)")+labs(fill = "Pattern type")+
geom_point(data=PEAK, aes(x=GammaPeakDepth, y=GammaPeakRichness, shape=datasettotsimpl))+
scale_x_reverse()+
labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+
ggtitle("(c) Gamma richness")
Figure 3: DDG of submerged macrophytes for alpha (a), beta (b) and gamma richness (c). For alpha richness, lines show the mean alpha richness per lake and year with their corresponding standard deviation; the single richness peaks (=DGG measures) are depicted as points. The different dataset levels can be distinguished by line type and point shape. Points and dashed line: Biotic dataset of all available macrophyte mapping (biodiversity dataset); triangles and solid line: subset of biotic dataset, where also abiotic data is available (environmental & biodiversity dataset).
Herberich test
Makroph_Lake_DepthS$TiefeFact <- as.factor(Makroph_Lake_DepthS$Tiefe)
# Alpha richness
aov1 = aov(ALPHA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model
Heteroaov1 <- glht(aov1,mcp(TiefeFact="Tukey") , vcov=vcovHC)
#summary(Heteroaov1) #Studies sites do not significantly differ in ....
Herb1<-confint(Heteroaov1) %>%
broom::tidy() %>%
filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups
ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) +
geom_hline(yintercept=0, linetype="11", colour="grey60") +
geom_errorbar(width=0.1) +
geom_point() +
coord_flip()+
ggtitle("95% family-wise confidence level - alpha richness")
#beta richness
aov2 = aov(BETA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model
Heteroaov2 <- glht(aov2,mcp(TiefeFact="Tukey") , vcov=vcovHC)
#summary(Heteroaov2) #Studies sites do not significantly differ in ....
Herb2<-confint(Heteroaov2) %>%
broom::tidy() %>%
filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups
ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) +
geom_hline(yintercept=0, linetype="11", colour="grey60") +
geom_errorbar(width=0.1) +
geom_point() +
coord_flip()+
ggtitle("95% family-wise confidence level - beta richness")
#Gamma richness
aov3 = aov(GAMMA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model
Heteroaov3 <- glht(aov3,mcp(TiefeFact="Tukey") , vcov=vcovHC)
#summary(Heteroaov3) #Studies sites do not significantly differ in ....
Herb3<-confint(Heteroaov3) %>%
broom::tidy() %>%
filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups
ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) +
geom_hline(yintercept=0, linetype="11", colour="grey60") +
geom_errorbar(width=0.1) +
geom_point() +
coord_flip() +
ggtitle("95% family-wise confidence level - gamma richness")
ggarrange(Herb1,Herb2,Herb3, nrow = 3,labels=c("(a)","(b)","(c)"))
Figure 4: Simultaneous tests for linear models with multiple comparisons of means using Tukey contrasts that are robust under non-normality, heteroscedasticity and variable sample size (Herberich et al. 2010) to check significant differences between depths within richness components. Results are plotted for alpha richness (a), beta richness (b) and gamma richness(c).
DDG pattern types
Table 2: Overview about DDG patterns of decreasing curves, hump-shaped curves (different peak depths) and increasing curves. Shows number of field campaigns (lake*year) showing a distinct depth pattern for each richness component.
## Organize DDG into Depth patterns
#For Alpha Richness peak
PEAK$DepthGroup<-if_else(PEAK$AlphaPeakDepth>-1,
"Decreasing curve (Peak: >-1m)",
if_else(PEAK$AlphaPeakDepth>-2,
"Hump-shaped (Peak: -1- -2m)",
if_else(PEAK$AlphaPeakDepth>-4,
"Hump-shaped (Peak: -2- -4m)",
"Increasing curve (Peak: <-4m)")))
PEAKCLASSAlpha<-PEAK %>% group_by(DepthGroup) %>%
dplyr::summarise(Alpha_N=n_distinct(interaction(Lake,YEAR)))
#For Gamma Richness peak
PEAKCLASS2<-PEAK %>% group_by(GammaPeakDepth) %>%
dplyr::summarise(datasets=n_distinct(interaction(Lake,YEAR)))
PEAKCLASS2$DepthGroup<-if_else(PEAKCLASS2$GammaPeakDepth>-1,
"Decreasing curve (Peak: >-1m)",
if_else(PEAKCLASS2$GammaPeakDepth>-2,
"Hump-shaped (Peak: -1- -2m)",
if_else(PEAKCLASS2$GammaPeakDepth>-4,
"Hump-shaped (Peak: -2- -4m)",
"Increasing curve (Peak: <-4m)")))
PEAKCLASSGamma<-PEAKCLASS2 %>% group_by(DepthGroup) %>% dplyr::summarise(Gamma_N=sum(datasets))
#For Beta Richness peak
PEAKCLASS3<-PEAK %>% group_by(BetaPeakDepth) %>%
dplyr::summarise(datasets=n_distinct(interaction(Lake,YEAR)))
PEAKCLASS3$DepthGroup<-if_else(PEAKCLASS3$BetaPeakDepth>-1,
"Decreasing curve (Peak: >-1m)",
if_else(PEAKCLASS3$BetaPeakDepth>-2,
"Hump-shaped (Peak: -1- -2m)",
if_else(PEAKCLASS3$BetaPeakDepth>-4,
"Hump-shaped (Peak: -2- -4m)",
"Increasing curve (Peak: <-4m)")))
PEAKCLASSBeta<-PEAKCLASS3 %>% group_by(DepthGroup) %>% dplyr::summarise(Beta_N=sum(datasets))
#Combine results
PEAKCLASS<-transpose(left_join(left_join(PEAKCLASSAlpha,PEAKCLASSBeta,
by="DepthGroup"),PEAKCLASSGamma,by="DepthGroup"),
keep.names = "col", make.names = "DepthGroup")
rownames(PEAKCLASS) <- PEAKCLASS[,1]
PEAKCLASS[,1] <- NULL
PEAKCLASS
Correlations between DDG metrices
peakvalues<-PEAK[c(18,3,5,14,15,16,17)]
colnames(peakvalues) <- c("gamma richness", "D(alpha,max)", "R(alpha,max)", "D(gamma,max)", "R(gamma,max)", "D(beta,max)", "R(beta,max)")
chart.Correlation(peakvalues, histogram=F, pch=9, method = "pearson")
Figure 5: Pearson correlations between DDG measures of different species richness components. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
Chi-square test of DDG pattern
Chi-square is done to see siginificant differences in frequency among the DDG pattern types (Table 2) to see if they are significantly different for each richness component.
#Chi-square test
chi<-chisq.test(PEAKCLASS, simulate.p.value=TRUE)
chi
Depth diversity gradients of macrophytes: drivers
Data representativeness of nested subsets
To show that the diversity metrics of the environmental & biodiversity dataset are representative for the diversity metrics of biodiversity dataset we applied the PERMANOVA test adonis2, using the R package ‘vegan’ which compares centroids and the variance (Oksanen et al. 2019). A non-significant result (p>0.05) confirms that centroids and variance of two groups are not different. The results show that the Environmental & biodiversity dataset (N=27) is representative for the Biodiversity dataset (N=100).
adonis2(PEAK[,c(3,5)]~datasettotsimpl, data=PEAK, by = NULL) #ALPHAPEAK
adonis2((PEAK[,c(14,15)])~datasettotsimpl, data=PEAK, by = NULL) #GammaPEAK
adonis2(scale(PEAK[,c(16,17)])~datasettotsimpl, data=PEAK, by = NULL) #BetaPEAK
Correlations between drivers
PEAK_Chem_norm_BA <- PEAK_Chem_norm %>% ungroup()%>% dplyr::filter(WLF!="NA") %>% filter(SAC!="NA")
values<-PEAK_Chem_norm_BA[c(25,21,23,31,32,29,30,16,19,3:15)]
colnames(values) <- c("gamma", "D(alp,max)", "R(alp,max)", "D(gam,max)", "R(gam,max)", "D(bet,max)", "R(bet,max)", "Area","WLF","Chloride","Cond","Ntot","NH4N", "NO3N","Osdiss","Ptot","pH","SiO2","Temp","Transp","SAC","Tempsd")
chart.Correlation(values,
histogram=F, pch=8, method = "pearson")
Figure 6: Pearson correlation between normalized chemical-physical values & DDG measures for all richness components. Environmental & biodiversity dataset is used. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
values<-PEAK_Chem_norm[c(25,21,23,31,32,29,30,16,19,3:15)]
colnames(values) <- c("gamma", "D(alp,max)", "R(alp,max)", "D(gam,max)", "R(gam,max)", "D(bet,max)", "R(bet,max)", "Area","WLF","Chloride","Cond","Ntot","NH4N", "NO3N","Osdiss","Ptot","pH","SiO2","Temp","Transp","SAC","Tempsd")
chart.Correlation(values,
histogram=F, pch=9, method = "pearson")
Figure 7: Pearson correlation between normalized chemical-physical values & DDG measures for all richness components. Biodiversity dataset is used, zero values are ignored. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
Principle component analysis
#environmental & biodiversity dataset
label_lake<- PEAK_Chem_norm[complete.cases(PEAK_Chem_norm[,c(3:16,19)]),]
lak.pca <- prcomp(na.omit(label_lake[,c(3:16,19)]),center = TRUE, scale. = TRUE)
#print(lak.pca)
#summary(lak.pca)
PCA <- (data.frame(label_lake$Lake))
PCA$YEAR <- label_lake$YEAR
PCA$PC1 <- lak.pca$x[,1]
PCA$PC2 <- lak.pca$x[,2]
PCA$PC3 <- lak.pca$x[,3]
PCA$PC4 <- lak.pca$x[,4]
names(PCA)[1]<-"Lake"
PEAK_PCA<-merge(PCA, PEAK_Chem_norm, by=c("Lake", "YEAR"))
#Plot
p1<- ggbiplot(lak.pca, choices = 1:2, obs.scale = 1, var.scale = 1,
labels=label_lake$Lake, arrow.color = "#FF0000",
ellipse = TRUE, cicle = TRUE)
p2<- ggbiplot(lak.pca, choices = 3:4, obs.scale = 1, var.scale = 1,
labels=label_lake$Lake, arrow.color = "#FF0000",
ellipse = TRUE, cicle = TRUE)
figure <- ggarrange(p1,p2,
labels = c("(a)","(b)"),
ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom", align = "hv")
figure
Figure 8: PCA with all environmental parameters of the environemnatal & biodiversity dataset. PCA axes 1-4 cover 0.8086% of variation. Axes 1 & 2 are plotted in panel (a), axes 3 & 4 in panel (b).
GAMM for gamma richness & beta and gamma DDG measures
#Beta
# gam_BetaPeakDepth <- gamm4(BetaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4),
# random= ~(1|Lake), # package gamm4
# data=PEAK_PCA)
# summary(gam_BetaPeakDepth$gam)
gam_BetaPeakRichness <- gamm4(BetaPeakRichness ~ s(PC3),
random= ~(1|Lake), # package gamm4
data=PEAK_PCA)
#summary(gam_BetaPeakRichness$gam)
#Gamma
# gam_GammaPeakDepth <- gamm4(GammaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4),
# random= ~(1|Lake), # package gamm4
# data=PEAK_PCA)
# summary(gam_GammaPeakDepth$gam)
gam_GammaPeakRichness <- gamm4(GammaPeakRichness ~ s(PC1),
random= ~(1|Lake), # package gamm4
data=PEAK_PCA)
#summary(gam_GammaPeakRichness$gam)
# Gamma richness
gam_Gamma <- gamm4(GAMMA ~ s(PC3), #s(PC1)+s(PC2)+s(PC3)+s(PC4)
random= ~(1|Lake), # package gamm4
data=PEAK_PCA)
#summary(gam_Gamma$gam)
#Plot
par(mfrow=c(3,1))
plot(gam_Gamma$gam, residuals=TRUE, main = "Gamma richness",shade = T,seWithMean=T)
text(-2.5, 4.0, "(a)", cex=2)
plot(gam_BetaPeakRichness$gam,residuals=TRUE, main = expression(R[beta][max]),
shade = T,seWithMean=T)
text(-2.5, 4.0, "(b)", cex=2)
plot(gam_GammaPeakRichness$gam,residuals=TRUE, main = expression(R[gamma][max]),
shade = T,seWithMean=T)
text(-3.5, 5.0, "(c)", cex=2)
Figure 9: Minimal GAMMs for response variables Gamma richness (a), R(beta,max) (b) and R(gamma, max) (c). Low r_square (adj) were found: 0.213 (a); 0.0124 (b); 0.265 (c).
Alternative GAMM with parameters selected on expert knowledge
gam_AlphaPeakDepth_alternative <- gamm4(AlphaPeakDepth ~
s(Conduct)+
#s(Ntot) +
#s(Ptot)+
#s(Temp)+
s(Transp)+
#s(Area_ha)+
s(WLF),
random= ~(1|Lake), # package gamm4
data=PEAK_Chem_norm_BA)
#summary(gam_AlphaPeakDepth_alternative$gam)
gam_AlphaPeakRichness_alternative <- gamm4(AlphaPeakRichness ~
#s(Conduct)+
#s(Ntot) +
#s(Ptot)+
#s(Temp)+
s(Transp)+
s(Area_ha),#+
#s(WLF),
random= ~(1|Lake), # package gamm4
data=PEAK_Chem_norm_BA)
#summary(gam_AlphaPeakRichness_alternative$gam)
# plot(gam_AlphaPeakRichness_alternative$gam, residuals=TRUE, main = "AlphaPeakRichness",shade = T,seWithMean=T)
p1<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "Transp") +
geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Transp"), alpha = 0.2) +
geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Transp"), alpha = 0.2) +
theme(legend.position="none")+
ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+
ggtitle("")+ theme(legend.title = element_blank())
p2<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "WLF") +
geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "WLF"), alpha = 0.2) +
geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "WLF"), alpha = 0.2) +
theme(legend.position="none")+
ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+
ggtitle("")+ theme(legend.title = element_blank())
p3<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "Conduct") +
geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Conduct"), alpha = 0.2) +
geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Conduct"), alpha = 0.2) +
theme(legend.position="none")+
ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+
ggtitle("")+ theme(legend.title = element_blank())
p4<-plotGAMM(gam_AlphaPeakRichness_alternative, smooth.cov = "Transp") +
geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Transp"), alpha = 0.2) +
geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Transp"), alpha = 0.2) +
theme(legend.position="none")+
ylab(expression(predicted_R[alpha][max]))+
ggtitle("")+ theme(legend.title = element_blank())
p5<-plotGAMM(gam_AlphaPeakRichness_alternative, smooth.cov = "Area_ha") +
geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Area_ha"), alpha = 0.2) +
geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Area_ha"), alpha = 0.2) +
theme(legend.position="none")+
ylab(expression(predicted_R[alpha][max]))+
ggtitle("")+ theme(legend.title = element_blank())
ggarrange(p1,p2,p3,p4,p5,
ncol=3, nrow=2,
#heights = c(4,3,3),
common.legend = T, legend="top",
labels = c("(a)","(b)","(c)", "(d)","(e)"),
#widths = c(1.4,1,1,1),
align = "hv"
)
Figure 10: Alternative GAMMs without PCA for selected parameters. The selection was based on expert knowledge on highly likly influencing factors and excluding high correlations (see Figure 6 and 7). We selected Conductivity (representative for Chloride), Ntot (representative for NH4 and NO3), Ptot, Temperature (representative for Tempsd), Transp (representative for SAC), Area and WLF. pH, O2diss & SiO2 were excluded. The resulting GAMM for D(alpha,max) 0.697 (a-c) has a R² lower (R²=0.697) than in the analysis using the PCA axes as variables (R²=0.73). GAMM for R(alpha,max) (d-f) shows also a slightly lower R² (R²=0.432) than in the analysis using the PCA axes as variables (R²=0.44).
Depth diversity gradients of macrophytes: recent shifts
Invariability analysis
LAKECHANGEInv<-PEAK %>%
dplyr::group_by(Lake) %>% #summarize informtion for lakes (over timeseries)
dplyr::summarise(NYEAR=n_distinct(YEAR),
Dalphamax=mean(-AlphaPeakDepth)/sd(-AlphaPeakDepth),
Ralphamax=mean(AlphaPeakRichness)/sd(AlphaPeakRichness),
Dbetamax=mean(-BetaPeakDepth)/sd(-BetaPeakDepth),
Rbetamax=mean(BetaPeakRichness)/sd(BetaPeakRichness),
Dgammamax=mean(-GammaPeakDepth)/sd(-GammaPeakDepth),
Rgammamax=mean(GammaPeakRichness)/sd(GammaPeakRichness),
xGammaRichness=mean(GAMMA)/sd(GAMMA)
)%>%
filter(NYEAR>3) %>% #For timeseries dataset
arrange(Dalphamax) %>%
filter_all(all_vars(!is.infinite(.)))
#Plot
ggplot(LAKECHANGEInv %>% gather(Type, "Inv", 3:9)) +
geom_histogram(aes(x=Inv), col="white", fill="black",binwidth = 1)+
facet_wrap(~Type, ncol = 3)+
xlab("Invariability coefficient")
Figure 11: Histograms for Invariability coefficients for single DDG measures (Dmax and Rmax of alpha, beta and gamma richness) and gamma richness for the biodiversity timeseries dataset. High invariability means a high stability within the timeseries.
Temporal trend of gamma richness
General
BTD <- PEAK %>% group_by(Lake) %>% filter(n_distinct(YEAR)>3) %>%
mutate(YEAR = as.numeric(as.character(YEAR)))
formula <- y ~ x
ggplot(data=BTD,
aes(x=YEAR, y=GAMMA, label=sprintf("%0.2f", round(YEAR, digits = 0)))) +
xlab("year")+ylab("GAMMA richness")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
#label.x.npc = "left", label.y.npc = 0.15,
formula = formula, parse = TRUE, size = 3)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 4), sep = "")),
label.x = 'left', label.y = 0.35, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016))
Figure 12: Temporal change of gamma richness from biodiversity time series dataset for all lakes together. Points show individual values and the blue line is a linear model.
Individual lakes
formula <- y ~ x
ggplot(BTD, aes(x=factor(YEAR), y=GAMMA, group=factor(Lake))) +
xlab("")+ylab("GAMMA richness")+ylim(0,40)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +
facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F, formula = y~x) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "left", label.y = 30,
formula = formula, parse = TRUE, size = 3)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x ="right", label.y = "bottom", size = 3,
vstep = 0.0)
Figure 13: Temporal change of gamma richness from biodiversity time series dataset for all individual lakes. Points show single values and the blue line is a linear model per lake.
Temporal trend of DDG measures
General
#formula <- y ~ x
A1<-ggplot(BTD,
aes(x=YEAR, y=AlphaPeakDepth)) +
xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
#label.x = "left", label.y = 0.15,
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom", size = 3,vstep = 0.0)
A2<-ggplot(BTD,
aes(x=YEAR, y=AlphaPeakRichness)) +
xlab("year")+ylab(expression(R[gamma][max]))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) + ylim(0,30)+
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom",
size = 3,vstep = 0.0)
B1<-ggplot(BTD,
aes(x=YEAR, y=BetaPeakDepth)) +
xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
#label.x = "left", label.y = 0.15,
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom", size = 3,vstep = 0.0)
B2<-ggplot(BTD,
aes(x=YEAR, y=BetaPeakRichness)) +
xlab("year")+ylab(expression(R[gamma][max]))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) + ylim(0,30)+
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom",
size = 3,vstep = 0.0)
C1<-ggplot(BTD,
aes(x=YEAR, y=GammaPeakDepth)) +
xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
#label.x = "left", label.y = 0.15,
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom", size = 3,vstep = 0.0)
C2<-ggplot(BTD,
aes(x=YEAR, y=GammaPeakRichness)) +
xlab("year")+ylab(expression(R[gamma][max]))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) + ylim(0,30)+
geom_smooth(method = "lm", se = T) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
formula = formula, parse = TRUE, size = 3)+
scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")),
label.x = "right", label.y = "bottom",
size = 3,vstep = 0.0)
ggarrange(A1,A2,B1,B2,C1,C2, ncol=2, nrow = 3,labels=c("(a)","(b)","(c)","(d)","(e)","(f)"))
Figure 14: Temporal change of DDG metrices from biodiversity time series dataset for all lakes at once. Points show single values and the blue line is a linear model. In panel (a,c,e) temporal change of D(max) is shown and in panel (b,d,f) the corresponding course of R(max) is depicted.
Individual lakes
formula <- y ~ x
T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=AlphaPeakDepth, group=factor(Lake)))+
xlab("")+ylab("AlphaPeakDepth")+ylim(-4,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=AlphaPeakRichness, group=factor(Lake)))+
xlab("")+ylab("AlphaPeakRichness")+ylim(0,15)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 15: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(alpha,max) is shown for all lakes and in panel (b) the corresponding course of R(alpha,max) is depicted.
formula <- y ~ x
T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=BetaPeakDepth, group=factor(Lake)))+
xlab("")+ylab("BetaPeakDepth")+ylim(-4,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=BetaPeakRichness, group=factor(Lake)))+
xlab("")+ylab("BetaPeakRichness")+ylim(0,25)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 16: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(beta,max) is shown for all lakes and in panel (b) the corresponding course of R(beta,max) is depicted.
formula <- y ~ x
T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=GammaPeakDepth, group=factor(Lake)))+
xlab("")+ylab("GammaPeakDepth")+ylim(-4,1)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=GammaPeakRichness, group=factor(Lake)))+
xlab("")+ylab("GammaPeakRichness")+ylim(0,30)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) +
geom_smooth(method = "lm", se = F) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x = "left", label.y = "top",
formula = formula, parse = TRUE, size = 3,vstep = 0.0)+
stat_fit_glance(method = 'lm',
method.args = list(formula = formula),
#geom = 'text',
aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")),
label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0)
ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 17: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(gamma,max) is shown for all lakes and in panel (b) the corresponding course of R(gamma,max) is depicted.
Literature
Herberich, E. et al. 2010. A Robust Procedure for Comparing Multiple Means under Heteroscedasticity in Unbalanced Designs (F Rapallo, Ed.). - PLoS ONE 5: e9788.
AnneLew/Lewerentz-etal_2021_Macrophytes-DDG documentation built on Dec. 17, 2021, 8:52 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
#Load packages devtools::load_all(path = ".") library(tidyverse) library(knitr) library(ggpubr) library(PerformanceAnalytics) library(data.table) library(kableExtra) library(knitr) library(ggpmisc) library(vegan) library(gamm4) #library(corrplot) library(voxel) #GAMM plotting using ggplot2 library(multcomp) #Herberich Test library(multcompView) #Herberich Test library(sandwich) #Herberich Test library(broom) #Herberich Test #devtools::install_github("vqv/ggbiplot") library(ggbiplot) library(lubridate) ## ggplot theme updates source(paste0(here::here(),"/R/theme.R"))
General information
Lake morphology
Table 1: Morphology of selected lakes: Lake surface area [ha] (Area_ha) and depths of deepest point of the lake [m] (maxDepth_m)
lakes<-unique(Makroph_Lake_DepthS$Lake) Morpho<-Morphology[Morphology$Name_Makro_short %in% lakes, ] Morpho %>% dplyr::select(Lake, Area_ha, maxDepth_m)
Depth diversity gradients (DDG) of macrophytes: shape
Correlations between richness components
richness<-Makroph_Lake_DepthS[c(97,101,99)] colnames(richness) <- c("alpha richness", "beta richness", "gamma richness") chart.Correlation(richness, histogram=F, pch=9, method = "pearson")
Figure 1: Correlations between diversity metrices (method = pearson). Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
DDG of all field campaigns
A1<-ggplot(data=Makroph_Lake_DepthS)+ geom_line(aes(x=Tiefe, y=ALPHA, group=interaction(Lake, YEAR)),col="grey")+ scale_x_reverse()+ geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mAlpha))+ geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mAlpha-sdAlpha, ymax=mAlpha+sdAlpha), width=.1)+ ylab("Alpha richness")+xlab(("Depth (m)")) B1<-ggplot(data=Makroph_Lake_DepthS)+ geom_line(aes(x=Tiefe, y=BETA, group=interaction(Lake, YEAR)),col="grey")+ scale_x_reverse()+ geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mBeta))+ geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mBeta-sdBeta, ymax=mBeta+sdBeta), width=.1)+ ylab("Beta richness")+xlab(("Depth (m)")) C1<-ggplot(data=Makroph_Lake_DepthS)+ geom_line(aes(x=Tiefe, y=GAMMA, group=interaction(Lake, YEAR)),col="grey")+ scale_x_reverse()+ geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mGamma))+ geom_errorbar(data=Makroph_Depth,aes(x=(Tiefe),ymin=mGamma-sdGamma, ymax=mGamma+sdGamma), width=.1)+ ylab("Gamma richness")+xlab(("Depth (m)")) ggarrange(A1,B1,C1,ncol=3,nrow=1,labels=c("(a)","(b)","(c)"))
Figure 2: Depth diversity gradients of macrophytes for alpha (a), beta (b) and gamma richness (c) with mean as black line and sd as black bars. Each single grey line is one field campaign (lake*year).
DDG per lake
### ALPHA & Peak Plot ggplot(data=Makroph_Lake_DepthS, aes(x=(Tiefe), y=ALPHA, col=factor(YEAR), group=interaction(Lake,YEAR)))+ geom_line(aes(linetype=datasettotsimpl))+ scale_linetype_manual(values=c("dotted", "solid"))+ facet_wrap(~ Lake, ncol=5)+ ylab("Alpha richness")+ xlab("")+labs(fill = "Pattern type")+ylim(0,15)+ geom_errorbar(data=Makroph_Lake_DepthS,aes(ymin=ALPHA-ALPHAsd, ymax=ALPHA+ALPHAsd), width=.1)+ geom_point(data=PEAK, aes(x=(AlphaPeakDepth), y=AlphaPeakRichness, shape=datasettotsimpl))+ scale_x_reverse()+ labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+ ggtitle("(a) Alpha richness") ### BETA & Peak Plot ggplot(Makroph_Lake_DepthS, aes(x=Tiefe, y=BETA, group=interaction(Lake,YEAR), col=factor(YEAR), linetype=datasettotsimpl))+ geom_line()+facet_wrap(~ Lake, ncol=5)+ scale_linetype_manual(values=c("dotted", "solid"))+ ylab("Beta richness")+ xlab("Depth (m)")+labs(fill = "Pattern type")+ geom_point(data=PEAK, aes(x=BetaPeakDepth, y=BetaPeakRichness, shape=datasettotsimpl))+ scale_x_reverse()+ labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+ ggtitle("(b) Beta richness") ### GAMMA & Peak Plot ggplot(Makroph_Lake_DepthS, aes(x=Tiefe, y=GAMMA, group=interaction(Lake,YEAR), col=factor(YEAR), linetype=datasettotsimpl))+ geom_line()+facet_wrap(~ Lake, ncol=5)+ scale_linetype_manual(values=c("dotted", "solid"))+ ylab("Gamma richness")+ xlab("Depth (m)")+labs(fill = "Pattern type")+ geom_point(data=PEAK, aes(x=GammaPeakDepth, y=GammaPeakRichness, shape=datasettotsimpl))+ scale_x_reverse()+ labs(linetype="Peak dataset type",shape="DDG dataset type", colour="Year")+ ggtitle("(c) Gamma richness")
Figure 3: DDG of submerged macrophytes for alpha (a), beta (b) and gamma richness (c). For alpha richness, lines show the mean alpha richness per lake and year with their corresponding standard deviation; the single richness peaks (=DGG measures) are depicted as points. The different dataset levels can be distinguished by line type and point shape. Points and dashed line: Biotic dataset of all available macrophyte mapping (biodiversity dataset); triangles and solid line: subset of biotic dataset, where also abiotic data is available (environmental & biodiversity dataset).
Herberich test
Makroph_Lake_DepthS$TiefeFact <- as.factor(Makroph_Lake_DepthS$Tiefe) # Alpha richness aov1 = aov(ALPHA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model Heteroaov1 <- glht(aov1,mcp(TiefeFact="Tukey") , vcov=vcovHC) #summary(Heteroaov1) #Studies sites do not significantly differ in .... Herb1<-confint(Heteroaov1) %>% broom::tidy() %>% filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) + geom_hline(yintercept=0, linetype="11", colour="grey60") + geom_errorbar(width=0.1) + geom_point() + coord_flip()+ ggtitle("95% family-wise confidence level - alpha richness") #beta richness aov2 = aov(BETA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model Heteroaov2 <- glht(aov2,mcp(TiefeFact="Tukey") , vcov=vcovHC) #summary(Heteroaov2) #Studies sites do not significantly differ in .... Herb2<-confint(Heteroaov2) %>% broom::tidy() %>% filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) + geom_hline(yintercept=0, linetype="11", colour="grey60") + geom_errorbar(width=0.1) + geom_point() + coord_flip()+ ggtitle("95% family-wise confidence level - beta richness") #Gamma richness aov3 = aov(GAMMA ~ (TiefeFact), data=Makroph_Lake_DepthS) #Fit an Analysis of Variance Model Heteroaov3 <- glht(aov3,mcp(TiefeFact="Tukey") , vcov=vcovHC) #summary(Heteroaov3) #Studies sites do not significantly differ in .... Herb3<-confint(Heteroaov3) %>% broom::tidy() %>% filter(contrast %in% c("-0.5 - -1.5","-1.5 - -3","-3 - -5"))%>% #to select just neighboring groups ggplot(aes(x=contrast, y=estimate, ymin=conf.low, ymax=conf.high)) + geom_hline(yintercept=0, linetype="11", colour="grey60") + geom_errorbar(width=0.1) + geom_point() + coord_flip() + ggtitle("95% family-wise confidence level - gamma richness") ggarrange(Herb1,Herb2,Herb3, nrow = 3,labels=c("(a)","(b)","(c)"))
Figure 4: Simultaneous tests for linear models with multiple comparisons of means using Tukey contrasts that are robust under non-normality, heteroscedasticity and variable sample size (Herberich et al. 2010) to check significant differences between depths within richness components. Results are plotted for alpha richness (a), beta richness (b) and gamma richness(c).
DDG pattern types
Table 2: Overview about DDG patterns of decreasing curves, hump-shaped curves (different peak depths) and increasing curves. Shows number of field campaigns (lake*year) showing a distinct depth pattern for each richness component.
## Organize DDG into Depth patterns #For Alpha Richness peak PEAK$DepthGroup<-if_else(PEAK$AlphaPeakDepth>-1, "Decreasing curve (Peak: >-1m)", if_else(PEAK$AlphaPeakDepth>-2, "Hump-shaped (Peak: -1- -2m)", if_else(PEAK$AlphaPeakDepth>-4, "Hump-shaped (Peak: -2- -4m)", "Increasing curve (Peak: <-4m)"))) PEAKCLASSAlpha<-PEAK %>% group_by(DepthGroup) %>% dplyr::summarise(Alpha_N=n_distinct(interaction(Lake,YEAR))) #For Gamma Richness peak PEAKCLASS2<-PEAK %>% group_by(GammaPeakDepth) %>% dplyr::summarise(datasets=n_distinct(interaction(Lake,YEAR))) PEAKCLASS2$DepthGroup<-if_else(PEAKCLASS2$GammaPeakDepth>-1, "Decreasing curve (Peak: >-1m)", if_else(PEAKCLASS2$GammaPeakDepth>-2, "Hump-shaped (Peak: -1- -2m)", if_else(PEAKCLASS2$GammaPeakDepth>-4, "Hump-shaped (Peak: -2- -4m)", "Increasing curve (Peak: <-4m)"))) PEAKCLASSGamma<-PEAKCLASS2 %>% group_by(DepthGroup) %>% dplyr::summarise(Gamma_N=sum(datasets)) #For Beta Richness peak PEAKCLASS3<-PEAK %>% group_by(BetaPeakDepth) %>% dplyr::summarise(datasets=n_distinct(interaction(Lake,YEAR))) PEAKCLASS3$DepthGroup<-if_else(PEAKCLASS3$BetaPeakDepth>-1, "Decreasing curve (Peak: >-1m)", if_else(PEAKCLASS3$BetaPeakDepth>-2, "Hump-shaped (Peak: -1- -2m)", if_else(PEAKCLASS3$BetaPeakDepth>-4, "Hump-shaped (Peak: -2- -4m)", "Increasing curve (Peak: <-4m)"))) PEAKCLASSBeta<-PEAKCLASS3 %>% group_by(DepthGroup) %>% dplyr::summarise(Beta_N=sum(datasets)) #Combine results PEAKCLASS<-transpose(left_join(left_join(PEAKCLASSAlpha,PEAKCLASSBeta, by="DepthGroup"),PEAKCLASSGamma,by="DepthGroup"), keep.names = "col", make.names = "DepthGroup") rownames(PEAKCLASS) <- PEAKCLASS[,1] PEAKCLASS[,1] <- NULL PEAKCLASS
Correlations between DDG metrices
peakvalues<-PEAK[c(18,3,5,14,15,16,17)] colnames(peakvalues) <- c("gamma richness", "D(alpha,max)", "R(alpha,max)", "D(gamma,max)", "R(gamma,max)", "D(beta,max)", "R(beta,max)") chart.Correlation(peakvalues, histogram=F, pch=9, method = "pearson")
Figure 5: Pearson correlations between DDG measures of different species richness components. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
Chi-square test of DDG pattern
Chi-square is done to see siginificant differences in frequency among the DDG pattern types (Table 2) to see if they are significantly different for each richness component.
#Chi-square test chi<-chisq.test(PEAKCLASS, simulate.p.value=TRUE) chi
Depth diversity gradients of macrophytes: drivers
Data representativeness of nested subsets
To show that the diversity metrics of the environmental & biodiversity dataset are representative for the diversity metrics of biodiversity dataset we applied the PERMANOVA test adonis2, using the R package ‘vegan’ which compares centroids and the variance (Oksanen et al. 2019). A non-significant result (p>0.05) confirms that centroids and variance of two groups are not different. The results show that the Environmental & biodiversity dataset (N=27) is representative for the Biodiversity dataset (N=100).
adonis2(PEAK[,c(3,5)]~datasettotsimpl, data=PEAK, by = NULL) #ALPHAPEAK adonis2((PEAK[,c(14,15)])~datasettotsimpl, data=PEAK, by = NULL) #GammaPEAK adonis2(scale(PEAK[,c(16,17)])~datasettotsimpl, data=PEAK, by = NULL) #BetaPEAK
Correlations between drivers
PEAK_Chem_norm_BA <- PEAK_Chem_norm %>% ungroup()%>% dplyr::filter(WLF!="NA") %>% filter(SAC!="NA") values<-PEAK_Chem_norm_BA[c(25,21,23,31,32,29,30,16,19,3:15)] colnames(values) <- c("gamma", "D(alp,max)", "R(alp,max)", "D(gam,max)", "R(gam,max)", "D(bet,max)", "R(bet,max)", "Area","WLF","Chloride","Cond","Ntot","NH4N", "NO3N","Osdiss","Ptot","pH","SiO2","Temp","Transp","SAC","Tempsd") chart.Correlation(values, histogram=F, pch=8, method = "pearson")
Figure 6: Pearson correlation between normalized chemical-physical values & DDG measures for all richness components. Environmental & biodiversity dataset is used. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
values<-PEAK_Chem_norm[c(25,21,23,31,32,29,30,16,19,3:15)] colnames(values) <- c("gamma", "D(alp,max)", "R(alp,max)", "D(gam,max)", "R(gam,max)", "D(bet,max)", "R(bet,max)", "Area","WLF","Chloride","Cond","Ntot","NH4N", "NO3N","Osdiss","Ptot","pH","SiO2","Temp","Transp","SAC","Tempsd") chart.Correlation(values, histogram=F, pch=9, method = "pearson")
Figure 7: Pearson correlation between normalized chemical-physical values & DDG measures for all richness components. Biodiversity dataset is used, zero values are ignored. Significance levels of p-values: *p<0.1;**p<0.05;***p<0.01
Principle component analysis
#environmental & biodiversity dataset label_lake<- PEAK_Chem_norm[complete.cases(PEAK_Chem_norm[,c(3:16,19)]),] lak.pca <- prcomp(na.omit(label_lake[,c(3:16,19)]),center = TRUE, scale. = TRUE) #print(lak.pca) #summary(lak.pca) PCA <- (data.frame(label_lake$Lake)) PCA$YEAR <- label_lake$YEAR PCA$PC1 <- lak.pca$x[,1] PCA$PC2 <- lak.pca$x[,2] PCA$PC3 <- lak.pca$x[,3] PCA$PC4 <- lak.pca$x[,4] names(PCA)[1]<-"Lake" PEAK_PCA<-merge(PCA, PEAK_Chem_norm, by=c("Lake", "YEAR")) #Plot p1<- ggbiplot(lak.pca, choices = 1:2, obs.scale = 1, var.scale = 1, labels=label_lake$Lake, arrow.color = "#FF0000", ellipse = TRUE, cicle = TRUE) p2<- ggbiplot(lak.pca, choices = 3:4, obs.scale = 1, var.scale = 1, labels=label_lake$Lake, arrow.color = "#FF0000", ellipse = TRUE, cicle = TRUE) figure <- ggarrange(p1,p2, labels = c("(a)","(b)"), ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom", align = "hv") figure
Figure 8: PCA with all environmental parameters of the environemnatal & biodiversity dataset. PCA axes 1-4 cover 0.8086% of variation. Axes 1 & 2 are plotted in panel (a), axes 3 & 4 in panel (b).
GAMM for gamma richness & beta and gamma DDG measures
#Beta # gam_BetaPeakDepth <- gamm4(BetaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4), # random= ~(1|Lake), # package gamm4 # data=PEAK_PCA) # summary(gam_BetaPeakDepth$gam) gam_BetaPeakRichness <- gamm4(BetaPeakRichness ~ s(PC3), random= ~(1|Lake), # package gamm4 data=PEAK_PCA) #summary(gam_BetaPeakRichness$gam) #Gamma # gam_GammaPeakDepth <- gamm4(GammaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4), # random= ~(1|Lake), # package gamm4 # data=PEAK_PCA) # summary(gam_GammaPeakDepth$gam) gam_GammaPeakRichness <- gamm4(GammaPeakRichness ~ s(PC1), random= ~(1|Lake), # package gamm4 data=PEAK_PCA) #summary(gam_GammaPeakRichness$gam) # Gamma richness gam_Gamma <- gamm4(GAMMA ~ s(PC3), #s(PC1)+s(PC2)+s(PC3)+s(PC4) random= ~(1|Lake), # package gamm4 data=PEAK_PCA) #summary(gam_Gamma$gam) #Plot par(mfrow=c(3,1)) plot(gam_Gamma$gam, residuals=TRUE, main = "Gamma richness",shade = T,seWithMean=T) text(-2.5, 4.0, "(a)", cex=2) plot(gam_BetaPeakRichness$gam,residuals=TRUE, main = expression(R[beta][max]), shade = T,seWithMean=T) text(-2.5, 4.0, "(b)", cex=2) plot(gam_GammaPeakRichness$gam,residuals=TRUE, main = expression(R[gamma][max]), shade = T,seWithMean=T) text(-3.5, 5.0, "(c)", cex=2)
Figure 9: Minimal GAMMs for response variables Gamma richness (a), R(beta,max) (b) and R(gamma, max) (c). Low r_square (adj) were found: 0.213 (a); 0.0124 (b); 0.265 (c).
Alternative GAMM with parameters selected on expert knowledge
gam_AlphaPeakDepth_alternative <- gamm4(AlphaPeakDepth ~ s(Conduct)+ #s(Ntot) + #s(Ptot)+ #s(Temp)+ s(Transp)+ #s(Area_ha)+ s(WLF), random= ~(1|Lake), # package gamm4 data=PEAK_Chem_norm_BA) #summary(gam_AlphaPeakDepth_alternative$gam) gam_AlphaPeakRichness_alternative <- gamm4(AlphaPeakRichness ~ #s(Conduct)+ #s(Ntot) + #s(Ptot)+ #s(Temp)+ s(Transp)+ s(Area_ha),#+ #s(WLF), random= ~(1|Lake), # package gamm4 data=PEAK_Chem_norm_BA) #summary(gam_AlphaPeakRichness_alternative$gam) # plot(gam_AlphaPeakRichness_alternative$gam, residuals=TRUE, main = "AlphaPeakRichness",shade = T,seWithMean=T) p1<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "Transp") + geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Transp"), alpha = 0.2) + geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Transp"), alpha = 0.2) + theme(legend.position="none")+ ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+ ggtitle("")+ theme(legend.title = element_blank()) p2<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "WLF") + geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "WLF"), alpha = 0.2) + geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "WLF"), alpha = 0.2) + theme(legend.position="none")+ ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+ ggtitle("")+ theme(legend.title = element_blank()) p3<-plotGAMM(gam_AlphaPeakDepth_alternative, smooth.cov = "Conduct") + geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Conduct"), alpha = 0.2) + geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakDepth", x = "Conduct"), alpha = 0.2) + theme(legend.position="none")+ ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+ ggtitle("")+ theme(legend.title = element_blank()) p4<-plotGAMM(gam_AlphaPeakRichness_alternative, smooth.cov = "Transp") + geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Transp"), alpha = 0.2) + geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Transp"), alpha = 0.2) + theme(legend.position="none")+ ylab(expression(predicted_R[alpha][max]))+ ggtitle("")+ theme(legend.title = element_blank()) p5<-plotGAMM(gam_AlphaPeakRichness_alternative, smooth.cov = "Area_ha") + geom_point(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Area_ha"), alpha = 0.2) + geom_rug(data = PEAK_Chem_norm_BA, aes_string(y = "AlphaPeakRichness", x = "Area_ha"), alpha = 0.2) + theme(legend.position="none")+ ylab(expression(predicted_R[alpha][max]))+ ggtitle("")+ theme(legend.title = element_blank()) ggarrange(p1,p2,p3,p4,p5, ncol=3, nrow=2, #heights = c(4,3,3), common.legend = T, legend="top", labels = c("(a)","(b)","(c)", "(d)","(e)"), #widths = c(1.4,1,1,1), align = "hv" )
Figure 10: Alternative GAMMs without PCA for selected parameters. The selection was based on expert knowledge on highly likly influencing factors and excluding high correlations (see Figure 6 and 7). We selected Conductivity (representative for Chloride), Ntot (representative for NH4 and NO3), Ptot, Temperature (representative for Tempsd), Transp (representative for SAC), Area and WLF. pH, O2diss & SiO2 were excluded. The resulting GAMM for D(alpha,max) 0.697 (a-c) has a R² lower (R²=0.697) than in the analysis using the PCA axes as variables (R²=0.73). GAMM for R(alpha,max) (d-f) shows also a slightly lower R² (R²=0.432) than in the analysis using the PCA axes as variables (R²=0.44).
Depth diversity gradients of macrophytes: recent shifts
Invariability analysis
LAKECHANGEInv<-PEAK %>% dplyr::group_by(Lake) %>% #summarize informtion for lakes (over timeseries) dplyr::summarise(NYEAR=n_distinct(YEAR), Dalphamax=mean(-AlphaPeakDepth)/sd(-AlphaPeakDepth), Ralphamax=mean(AlphaPeakRichness)/sd(AlphaPeakRichness), Dbetamax=mean(-BetaPeakDepth)/sd(-BetaPeakDepth), Rbetamax=mean(BetaPeakRichness)/sd(BetaPeakRichness), Dgammamax=mean(-GammaPeakDepth)/sd(-GammaPeakDepth), Rgammamax=mean(GammaPeakRichness)/sd(GammaPeakRichness), xGammaRichness=mean(GAMMA)/sd(GAMMA) )%>% filter(NYEAR>3) %>% #For timeseries dataset arrange(Dalphamax) %>% filter_all(all_vars(!is.infinite(.))) #Plot ggplot(LAKECHANGEInv %>% gather(Type, "Inv", 3:9)) + geom_histogram(aes(x=Inv), col="white", fill="black",binwidth = 1)+ facet_wrap(~Type, ncol = 3)+ xlab("Invariability coefficient")
Figure 11: Histograms for Invariability coefficients for single DDG measures (Dmax and Rmax of alpha, beta and gamma richness) and gamma richness for the biodiversity timeseries dataset. High invariability means a high stability within the timeseries.
Temporal trend of gamma richness
General
BTD <- PEAK %>% group_by(Lake) %>% filter(n_distinct(YEAR)>3) %>% mutate(YEAR = as.numeric(as.character(YEAR))) formula <- y ~ x ggplot(data=BTD, aes(x=YEAR, y=GAMMA, label=sprintf("%0.2f", round(YEAR, digits = 0)))) + xlab("year")+ylab("GAMMA richness")+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), #label.x.npc = "left", label.y.npc = 0.15, formula = formula, parse = TRUE, size = 3)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 4), sep = "")), label.x = 'left', label.y = 0.35, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016))
Figure 12: Temporal change of gamma richness from biodiversity time series dataset for all lakes together. Points show individual values and the blue line is a linear model.
Individual lakes
formula <- y ~ x ggplot(BTD, aes(x=factor(YEAR), y=GAMMA, group=factor(Lake))) + xlab("")+ylab("GAMMA richness")+ylim(0,40)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F, formula = y~x) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x.npc = "left", label.y = 30, formula = formula, parse = TRUE, size = 3)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x ="right", label.y = "bottom", size = 3, vstep = 0.0)
Figure 13: Temporal change of gamma richness from biodiversity time series dataset for all individual lakes. Points show single values and the blue line is a linear model per lake.
Temporal trend of DDG measures
General
#formula <- y ~ x A1<-ggplot(BTD, aes(x=YEAR, y=AlphaPeakDepth)) + xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), #label.x = "left", label.y = 0.15, formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) A2<-ggplot(BTD, aes(x=YEAR, y=AlphaPeakRichness)) + xlab("year")+ylab(expression(R[gamma][max]))+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + ylim(0,30)+ geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) B1<-ggplot(BTD, aes(x=YEAR, y=BetaPeakDepth)) + xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), #label.x = "left", label.y = 0.15, formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) B2<-ggplot(BTD, aes(x=YEAR, y=BetaPeakRichness)) + xlab("year")+ylab(expression(R[gamma][max]))+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + ylim(0,30)+ geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) C1<-ggplot(BTD, aes(x=YEAR, y=GammaPeakDepth)) + xlab("year")+ylab(expression(D[gamma][max]))+ylim(-5,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), #label.x = "left", label.y = 0.15, formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) C2<-ggplot(BTD, aes(x=YEAR, y=GammaPeakRichness)) + xlab("year")+ylab(expression(R[gamma][max]))+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) + ylim(0,30)+ geom_smooth(method = "lm", se = T) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), formula = formula, parse = TRUE, size = 3)+ scale_x_continuous(breaks=c(2007,2010,2013,2016), limits = c(2006,2017))+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), label.x = "right", label.y = "bottom", size = 3,vstep = 0.0) ggarrange(A1,A2,B1,B2,C1,C2, ncol=2, nrow = 3,labels=c("(a)","(b)","(c)","(d)","(e)","(f)"))
Figure 14: Temporal change of DDG metrices from biodiversity time series dataset for all lakes at once. Points show single values and the blue line is a linear model. In panel (a,c,e) temporal change of D(max) is shown and in panel (b,d,f) the corresponding course of R(max) is depicted.
Individual lakes
formula <- y ~ x T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=AlphaPeakDepth, group=factor(Lake)))+ xlab("")+ylab("AlphaPeakDepth")+ylim(-4,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=AlphaPeakRichness, group=factor(Lake)))+ xlab("")+ylab("AlphaPeakRichness")+ylim(0,15)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 15: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(alpha,max) is shown for all lakes and in panel (b) the corresponding course of R(alpha,max) is depicted.
formula <- y ~ x T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=BetaPeakDepth, group=factor(Lake)))+ xlab("")+ylab("BetaPeakDepth")+ylim(-4,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=BetaPeakRichness, group=factor(Lake)))+ xlab("")+ylab("BetaPeakRichness")+ylim(0,25)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 16: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(beta,max) is shown for all lakes and in panel (b) the corresponding course of R(beta,max) is depicted.
formula <- y ~ x T1<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=GammaPeakDepth, group=factor(Lake)))+ xlab("")+ylab("GammaPeakDepth")+ylim(-4,1)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) T2<-ggplot(PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3), aes(x=factor(YEAR), y=GammaPeakRichness, group=factor(Lake)))+ xlab("")+ylab("GammaPeakRichness")+ylim(0,30)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ geom_point(alpha = 0.3) +facet_wrap(vars(Lake),ncol=4) + geom_smooth(method = "lm", se = F) + stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", label.y = "top", formula = formula, parse = TRUE, size = 3,vstep = 0.0)+ stat_fit_glance(method = 'lm', method.args = list(formula = formula), #geom = 'text', aes(label = paste("P-value = ", signif(..p.value.., digits = 2), sep = "")), label.x = 'left', label.y = "bottom", size = 3,vstep = 0.0) ggarrange(T1,T2, nrow=2, labels = c("(a)","(b)"))
Figure 17: Temporal change of DDG metrices from biodiversity time series dataset for all lakes. Points show annual values and the blue line is a linear model. In panel (a) temporal change of D(gamma,max) is shown for all lakes and in panel (b) the corresponding course of R(gamma,max) is depicted.
Literature
Herberich, E. et al. 2010. A Robust Procedure for Comparing Multiple Means under Heteroscedasticity in Unbalanced Designs (F Rapallo, Ed.). - PLoS ONE 5: e9788.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.