In AnneLew/Lewerentz-etal_2021_Macrophytes-DDG: Analysis of the Depth Diversity Gradient of macrophytes
knitr::opts_chunk$set(echo = T, message=FALSE, warning=FALSE, dpi=300)
devtools::load_all(path = ".")
#Load packages
library(tidyverse)
library(knitr)
library(ggpubr)
library(data.table)
library(kableExtra)
library(knitr)
library(vegan)
library(gamm4)
library(voxel)
library(gridExtra)
library(ggpmisc)
library(grid)
## ggplot theme
source(paste0(here::here(),"/R/theme.R"))
Material & methods
Data overview
Abiotic data
Chem_overview_SmallDatasetMorph<-Chem_uniform_LOIx %>%
filter(dataset=="SAK") %>%
filter(datasetWLF=="WLF") %>%
group_by(Lake)%>%
select(Area_ha, maxDepth_m, WLF) %>% unique() %>% ungroup() %>%
rename_(.dots=setNames(names(.), tolower(gsub("\\_", ".", names(.)))))%>%
summarise_each(funs(mean, median, min, max, sd),-lake) %>%
tidyr::gather(variable, value) %>%
tidyr::separate(variable, c("var", "stat"), sep = "\\_") %>%
tidyr::spread(var, value)
Chem_overview_SmallDataset<-Chem_uniform_LOIx %>% rowwise() %>%
filter(dataset=="SAK") %>%
filter(datasetWLF=="WLF") %>%
rename_(.dots=setNames(names(.), tolower(gsub("\\_", ".", names(.)))))%>%ungroup()%>%
select(chloride,conduct,nh4n,no3n,ntot,o2diss,ph,ptot,sac,sio2,temp,tempsd,transp,wlf)%>%
summarise_each(funs(mean, median, min, max, sd))%>%
tidyr::gather(variable, value) %>%
tidyr::separate(variable, c("var", "stat"), sep = "\\_") %>%
tidyr::spread(var, value)#%>%t()
Chem_overview_SmallD<-cbind(Chem_overview_SmallDatasetMorph,Chem_overview_SmallDataset[,2:15]) %>%
tibble::rownames_to_column() %>%
tidyr::gather(var, value, -rowname) %>%
tidyr::spread(rowname, value) %>%
rename(max="1",mean="2", median="3", min="4", sd="5")
df<-cbind(Chem_overview_SmallDatasetMorph,Chem_overview_SmallDataset[,2:15]) %>%
column_to_rownames(var = "stat")
dt<-type_convert(as_tibble(cbind(Parameters = names(df), t(df)))) %>% mutate_if(is.numeric, round, 2)
dt %>%
select(Parameters,min, max,mean,sd,median)
Biotic data, dataset levels overview
complete<-MakrophS_ALL %>%
dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N =n_distinct(Probestelle),Years_N = n_distinct(YEAR))
Level1<-MakrophS_ALL %>% filter(datasettot=="LEVEL3") %>% #group_by(datasettot) %>%
dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N = n_distinct(Probestelle),Years_N = n_distinct(YEAR))
Level_Timeseries<-MakrophS_ALL %>% group_by(Lake)%>%filter(n_distinct(YEAR)>3)%>% ungroup() %>%
dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N = n_distinct(Probestelle),Years_N = n_distinct(YEAR))
Level<- rbind(complete,Level1,Level_Timeseries)
Level$dataset<-c("Biotic","Biotic+Abiotic", "Timeseries")
Level[,c(6,1,2,3,4,5)]
Results
Species richness per dataset group
M<-MakrophS_ALL[5:96]
M[colSums(M)!=0] %>% ncol() #75
M3<-MakrophS_ALL %>% filter(datasettot=="LEVEL3")
M3<-M3[5:96]
M3[colSums(M3)!=0] %>% ncol() #57
MT<-MakrophS_ALL %>% group_by(Lake)%>%filter(n_distinct(YEAR)>3)%>% ungroup()
MT<-MT[5:96]
MT[colSums(MT)!=0] %>% ncol() #66
Species numbers in taxonomic groups
M<-MakrophS_ALL[5:96]
specpres<-M[colSums(M)!=0]
species <- as.data.frame(colnames(specpres))
species
names(species)<-"Taxon"
unique(left_join(species,Makroph[c(8,11)],by=c("Taxon"))) %>%
group_by(Erscheinungsform) %>% dplyr::summarise(N=n_distinct(Taxon))
Gamma richness overview: min, max, mean, sd
min(Makroph_Lake_ALL$GAMMA)
max(Makroph_Lake_ALL$GAMMA)
mean(Makroph_Lake_ALL$GAMMA)
sd(Makroph_Lake_ALL$GAMMA)
Depth diversity gradients
A1<-ggplot(data=Makroph_Depth)+
geom_ribbon(aes(x = Tiefe, ymax =mAlpha+sdAlpha, ymin =mAlpha-sdAlpha), alpha = 0.6, fill = "grey")+
geom_line(aes(x=(Tiefe), y=mAlpha), size=1)+
scale_x_reverse()+
ylab("Alpha richness")+xlab("Depth (m)")
formula1 = y ~ poly(x, 3, raw=TRUE)
A2<-ggplot(data = PEAK, aes(x=AlphaPeakDepth, y=AlphaPeakRichness))+
geom_point(aes(shape=datasettotsimpl),alpha=0.2)+
scale_shape_discrete(name = "s",
labels = c("DDG peak, Biodiversity dataset",
"DDG peak, Environmental & biodiversity dataset"))+
geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+
xlab(expression(D[alpha][max](m)))+
ylab(expression(R[alpha][max](N)))+
theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+
scale_x_reverse()+
xlim(-0.5,-5)
A3<-ggplot(data=PEAK, aes(x=AlphaPeakDepth))+
geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1
scale_x_reverse()+xlab("Depth (m)")+ylab(expression(D[alpha][max](counts)))
B1<-ggplot(data=Makroph_Depth)+
geom_ribbon(aes(x = Tiefe, ymax =mBeta+sdBeta, ymin =mBeta-sdBeta), alpha = 0.6, fill = "grey")+
geom_line(aes(x=(Tiefe), y=mBeta), size=1)+ scale_x_reverse()+
ylab("Beta richness")+xlab(("Depth (m)"))
B2<-ggplot(data = PEAK, aes(x=BetaPeakDepth, y=BetaPeakRichness))+
geom_point(aes(shape=datasettotsimpl),alpha=0.2)+
scale_shape_discrete(name = "s",
labels = c("DDG peak, Biodiversity dataset",
"DDG peak, Environmental & biodiversity dataset"))+
geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+
xlab(expression(D[beta][max](m)))+
ylab(expression(R[beta][max](N)))+
theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+
scale_x_reverse()+
xlim(-0.5,-5)
B3<-ggplot(data=PEAK, aes(x=BetaPeakDepth))+
geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1
scale_x_reverse()+xlab("Depth (m)")+
ylab(expression(D[beta][max](counts)))
C1<-ggplot(data=Makroph_Depth)+
geom_ribbon(aes(x = Tiefe, ymax =mGamma+sdGamma, ymin =mGamma-sdGamma), alpha = 0.6, fill = "grey")+
geom_line(aes(x=(Tiefe), y=mGamma), size=1, colour="black")+ scale_x_reverse()+
ylab("Gamma richness")+xlab(("Depth (m)"))
C2<-ggplot(data = PEAK, aes(x=GammaPeakDepth, y=GammaPeakRichness))+
geom_point(aes(shape=datasettotsimpl),alpha=0.2)+
scale_shape_discrete(name = "s",
labels = c("DDG peak, Biodiversity dataset",
"DDG peak, Environmental & biodiversity dataset"))+
geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+
xlab(expression(D[gamma][max](m)))+
ylab(expression(R[gamma][max](N)))+
theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+
scale_x_reverse()+
xlim(-0.5,-5)
C3<-ggplot(data=PEAK, aes(x=GammaPeakDepth))+
geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1
scale_x_reverse()+xlab("Depth (m)")+ylab(expression(D[gamma][max](counts)))
Fig2<-ggarrange(A1,B1,C1,A2,B2,C2,A3,B3,C3,
ncol=3,nrow=3,labels=c("(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)"),
heights = c(2,1.5,1.8),common.legend = T, legend="none", align = "hv")
# png("Fig2.png", width=16.6, height=18,units = "cm",res=300)
# Fig2
# dev.off
#
# ggsave("Fig2.pdf",plot=Fig1,width=16.6, height = 18, device=cairo_pdf, units = "cm")
windowsFonts("Arial" = windowsFont("Arial"))
Fig2
Drivers
Cor.test
cor.test(PEAK_Chem_norm$AlphaPeakRichness, PEAK_Chem_norm$Ntot)
GAMM incl PCA
################ PCA ###############################################
label_lake<- PEAK_Chem_norm[complete.cases(PEAK_Chem_norm[,c(3:16,19)]),] #LEVEL 3 data
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"))
Rotation <-lak.pca$rotation %>% as.data.frame()
Rotation$variable <- row.names(Rotation)
R1<-ggplot(data=Rotation)+
geom_bar(aes(y=variable,x=PC1,fill = PC1 > 0.4 | PC1< -0.4), stat='identity')+
xlim(-0.65,0.65)+
#theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
ylab("")+xlab("loading")+
ggtitle(" PC1 - 30.1% \nSiO2 & Cond axis")
R2<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC2,fill = PC2 > 0.4 | PC2< -0.4), stat='identity')+
xlim(-0.65,0.65)+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ylab("")+
ggtitle(" PC2 - 26.1%\nTemp & Ptot axis")+xlab("loading")+ylab("")
R3<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC3,fill = PC3 > 0.4 | PC3< -0.4), stat='identity')+
xlim(-0.65,0.65)+ylab("")+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
ggtitle(" PC3 - 13.3%\nTempsd – Chl axis")+xlab("loading")
R4<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC4,fill = PC4 > 0.4 | PC4< -0.4), stat='identity')+
xlim(-0.65,0.65)+ylab("")+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
ggtitle(" PC4 - 10.5%\nO2diss – SAC axis")+xlab("loading")
## GAMM
gam_AlphaPeakDepth <- gamm4(AlphaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4),
random= ~(1|Lake), # package gamm4
data=PEAK_PCA)
summary(gam_AlphaPeakDepth$gam)
gam_AlphaPeakRichness <- gamm4(AlphaPeakRichness ~ s(PC1),
random= ~(1|Lake), # package gamm4
data=PEAK_PCA)
summary(gam_AlphaPeakRichness$gam)
vars <- c("PC2", "PC4","PC3","PC1")
grob2 <- grobTree(textGrob("*** 36.6% dc", x=0.1, y=0.9, hjust=0,
gp=gpar(col="grey12", fontsize=13, fontface="italic")))
plotPC2<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC2") +
geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC2"), alpha = 0.2) +
geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC2"), alpha = 0.2) +
#scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) +
theme(legend.position="none")+
ylab("")+ylim(-3.5,-0.5)+
ggtitle("")+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
theme(legend.title = element_blank())+
annotation_custom(grob2)
grob4 <- grobTree(textGrob("** 30.3% dc", x=0.1, y=0.9, hjust=0,
gp=gpar(col="grey12", fontsize=13, fontface="italic")))
vars <- c("PC4")
plotPC4<- plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC4") +
geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC4"), alpha = 0.2) +
geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC4"), alpha = 0.2) +
#scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) +
theme(legend.position="none")+theme(legend.title = element_blank())+
ylab("")+ylim(-3.5,-0.5)+
#ylab(expression(predicted_D[alpha][max]))+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
ggtitle("")+
annotation_custom(grob4)
grob3 <- grobTree(textGrob("** 28.6% dc", x=0.1, y=0.9, hjust=0,
gp=gpar(col="grey12", fontsize=13, fontface="italic")))
vars <- c("PC3")
plotPC3<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC3") +
geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC3"), alpha = 0.2) +
geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC3"), alpha = 0.2) +
#scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) +
theme(legend.position="none")+theme(legend.title = element_blank())+
ylab("")+ylim(-3.5,-0.5)+
#ylab(expression(predicted_D[alpha][max]))+
theme(axis.title.y=element_blank(),axis.text.y=element_blank())+
ggtitle("")+
annotation_custom(grob3)
grob1 <- grobTree(textGrob("* 25.6% dc", x=0.1, y=0.9, hjust=0,
gp=gpar(col="grey12", fontsize=13, fontface="italic")))
vars <- c("PC1")
plotPC1<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC1") +
geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC1"), alpha = 0.2) +
geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC1"), alpha = 0.2) +
#scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) +
theme(legend.position="none")+
ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+
ggtitle("")+ theme(legend.title = element_blank())+
annotation_custom(grob1)
grob5 <- grobTree(textGrob("**", x=0.1, y=0.9, hjust=0,
gp=gpar(col="grey12", fontsize=13, fontface="italic")))
vars2 <- c("PC1")
plot2<-plotGAMM(gam_AlphaPeakRichness, smooth.cov = "PC1") +
geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakRichness", x = "PC1"), alpha = 0.2) +
geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakRichness", x = "PC1"), alpha = 0.2) +
#scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) +
theme(legend.position="none")+
ylab(expression(predicted_R[alpha][max]))+
#labs(tag = "E")+
ggtitle("")+ theme(legend.title = element_blank())+
annotation_custom(grob5)
empty<-plot(0, xaxt = 'n', yaxt = 'n', bty = 'n', pch = '', ylab = '', xlab = '')
#ggarrange(R1,R2,R3,R4,
# ncol=4, common.legend = T, legend="top", labels = c("(a)","(b)","(c)","(d)"), widths = c(1.3,1,1,1))
# ggarrange(plotPC1, plotPC2, plotPC3, plotPC4, plot2,
# ncol=4, nrow=3,
# heights = c(4,3,3),
# common.legend = T, legend="top",
# labels = c("(e)","(f)","(g)","(h)","(i)"),
# #widths = c(1.4,1,1,1),
# align = "hv"
# )
Fig3<-ggarrange(R1,R2,R3,R4, plotPC1, plotPC2, plotPC3, plotPC4, plot2,empty, empty, empty,
ncol=4, nrow=3,
heights = c(4,3,3),
common.legend = T, legend="top",
labels = c("(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)"),
#widths = c(1.4,1,1,1),
align = "hv"
)
Fig3
# ggsave("Fig3.pdf",plot=Fig3,width = 9, height = 8, device=cairo_pdf)
#
# png("Fig3.png", width=16.6, height=18,units = "cm",res=300)
# Fig3
# dev.off
Temporal change
Invariability analysis
LAKECHANGEInv<-PEAK %>%
dplyr::group_by(Lake) %>% #summarize informtion for lakes (over timeseries)
dplyr::summarise(NYEAR=n_distinct(YEAR),
AlphaPeakDepthInv=mean(-AlphaPeakDepth)/sd(-AlphaPeakDepth),
AlphaPeakRichnessInv=mean(AlphaPeakRichness)/sd(AlphaPeakRichness),
BetaPeakDepthInv=mean(-BetaPeakDepth)/sd(-BetaPeakDepth),
BetaPeakRichnessInv=mean(BetaPeakRichness)/sd(BetaPeakRichness),
GammaPeakDepthInv=mean(-GammaPeakDepth)/sd(-GammaPeakDepth),
GammaPeakRichnessInv=mean(GammaPeakRichness)/sd(GammaPeakRichness),
GammaRichnessInv=mean(GAMMA)/sd(GAMMA)
)%>%
filter(NYEAR>3) %>% #For timeseries dataset
arrange(AlphaPeakDepthInv) %>%
filter_all(all_vars(!is.infinite(.)))
#LAKECHANGEInv
# ggplot(LAKECHANGEInv %>% gather(Type, "Inv", 3:9)) +
# geom_histogram(aes(x=Inv), col="white", fill="black",binwidth = 1)+
# facet_wrap(~Type, ncol = 2)+
# xlab("Inverse Variation coefficient = Invariability coefficient")
LAKECHANGEInv%>%
dplyr::ungroup() %>%
select(-Lake, -NYEAR)%>%
dplyr::summarise_all(list(mean=mean,sd=sd), na.rm=TRUE)
Linear models
PEAK<-PEAK %>%mutate(YEAR=as.numeric(YEAR))
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
attributes(p) <- NULL
return(p)
}
GammaModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(GAMMA ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="GammaSlope", "pValue"="GammaPValue"))
AlphaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(AlphaPeakRichness ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="AlphaPeakRichnessSlope", "pValue"="AlphaPeakRichnessPValue"))
AlphaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(AlphaPeakDepth ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="AlphaPeakDepthSlope", "pValue"="AlphaPeakDepthPValue"))
AlphaPeakModel <- merge (AlphaPeakDepthModel, AlphaPeakRichnessModel, by="Lake")
BetaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(BetaPeakRichness ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="BetaPeakRichnessSlope", "pValue"="BetaPeakRichnessPValue"))
BetaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(BetaPeakDepth ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="BetaPeakDepthSlope", "pValue"="BetaPeakDepthPValue"))
BetaPeakModel <- merge (BetaPeakDepthModel, BetaPeakRichnessModel, by="Lake")
GammaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(GammaPeakRichness ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="GammaPeakRichnessSlope", "pValue"="GammaPeakRichnessPValue"))
GammaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>%
do({
mod = lm(GammaPeakDepth ~ YEAR, data = .)
data.frame(Slope = coef(mod)[2],pValue=lmp(mod))
}) %>% plyr::rename(c("Slope"="GammaPeakDepthSlope", "pValue"="GammaPeakDepthPValue"))
GammaPeakModel <- merge (GammaPeakDepthModel, GammaPeakRichnessModel, by="Lake")
PeakModel <- merge (AlphaPeakModel,BetaPeakModel,by="Lake")
PeakModel <- merge (PeakModel,GammaPeakModel,by="Lake")
Model <- merge (GammaModel,PeakModel,by="Lake")
ModelT<-Model %>% mutate(GammaTrend=ifelse(Model$GammaSlope>0,"+",ifelse(Model$GammaSlope<0,"-","0"))) %>%
mutate(GammapV=ifelse(GammaPValue<0.001,"***",ifelse(GammaPValue<0.01,"**",ifelse(GammaPValue<0.05,"*",ifelse(GammaPValue<0.1,".","NA")))))%>%
mutate(AlphaPeakDepthTrend=ifelse(Model$AlphaPeakDepthSlope>0.5*0,"+",ifelse(Model$AlphaPeakDepthSlope<0,"-","0"))) %>%
mutate(AlphaPeakDepthpV=ifelse(AlphaPeakDepthPValue<0.001,"***",ifelse(AlphaPeakDepthPValue<0.01,"**",ifelse(AlphaPeakDepthPValue<0.05,"*",ifelse(AlphaPeakDepthPValue<0.1,".","NA")))))%>%
mutate(AlphaPeakRichnessTrend=ifelse(Model$AlphaPeakRichnessSlope>0,"+",ifelse(Model$AlphaPeakRichnessSlope<0,"-","0")))%>%
mutate(AlphaPeakRichnesspV=ifelse(AlphaPeakRichnessPValue<0.001,"***",ifelse(AlphaPeakRichnessPValue<0.01,"**",ifelse(AlphaPeakRichnessPValue<0.05,"*",ifelse(AlphaPeakRichnessPValue<0.1,".","NA")))))%>%
mutate(BetaPeakDepthTrend=ifelse(Model$BetaPeakDepthSlope>0,"+",ifelse(Model$BetaPeakDepthSlope<0,"-","0"))) %>%
mutate(BetaPeakDepthpV=ifelse(BetaPeakDepthPValue<0.001,"***",ifelse(BetaPeakDepthPValue<0.01,"**",ifelse(BetaPeakDepthPValue<0.05,"*",ifelse(BetaPeakDepthPValue<0.1,".","NA")))))%>%
mutate(BetaPeakRichnessTrend=ifelse(Model$BetaPeakRichnessSlope>0,"+",ifelse(Model$BetaPeakRichnessSlope<0,"-","0")))%>%
mutate(BetaPeakRichnesspV=ifelse(BetaPeakRichnessPValue<0.001,"***",ifelse(BetaPeakRichnessPValue<0.01,"**",ifelse(BetaPeakRichnessPValue<0.05,"*",ifelse(BetaPeakRichnessPValue<0.1,".","NA")))))%>%
mutate(GammaPeakDepthTrend=ifelse(Model$GammaPeakDepthSlope>0,"+",ifelse(Model$GammaPeakDepthSlope<0,"-","0"))) %>%
mutate(GammaPeakDepthpV=ifelse(GammaPeakDepthPValue<0.001,"***",ifelse(GammaPeakDepthPValue<0.01,"**",ifelse(GammaPeakDepthPValue<0.05,"*",ifelse(GammaPeakDepthPValue<0.1,".","NA")))))%>%
mutate(GammaPeakRichnessTrend=ifelse(Model$GammaPeakRichnessSlope>0,"+",ifelse(Model$GammaPeakRichnessSlope<0,"-","0")))%>%
mutate(GammaPeakRichnesspV=ifelse(GammaPeakRichnessPValue<0.001,"***",ifelse(GammaPeakRichnessPValue<0.01,"**",ifelse(GammaPeakRichnessPValue<0.05,"*",ifelse(GammaPeakRichnessPValue<0.1,".","NA")))))
ModelTshort<-ModelT[c(1,16:29)] %>%
tidyr::unite(GammaLM, c("GammaTrend", "GammapV"))%>%
tidyr::unite(AlphaPeakDepthLM, c("AlphaPeakDepthTrend", "AlphaPeakDepthpV"))%>%
tidyr::unite(AlphaPeakRichnessLM, c("AlphaPeakRichnessTrend", "AlphaPeakRichnesspV"))%>%
tidyr::unite(BetaPeakDepthLM, c("BetaPeakDepthTrend", "BetaPeakDepthpV"))%>%
tidyr::unite(BetaPeakRichnessLM, c("BetaPeakRichnessTrend", "BetaPeakRichnesspV"))%>%
tidyr::unite(GammaPeakDepthLM, c("GammaPeakDepthTrend", "GammaPeakDepthpV"))%>%
tidyr::unite(GammaPeakRichnessLM, c("GammaPeakRichnessTrend", "GammaPeakRichnesspV"))
ModelTshort <- data.frame(lapply(ModelTshort, function(x) {gsub("_NA", "", x)}))
ModelTshort <- data.frame(lapply(ModelTshort, function(x) {gsub("_", "", x)}))
ModelTshort
Discussion
Summary figure
ggplot(data=Makroph_Depth)+
geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mAlpha), col="darkgreen", size=1)+
geom_line(data=Makroph_Depth, aes(x=(Tiefe), y=mBeta), col="green2", size=1)+
geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mGamma), col="lightgreen", size=1, colour="black")+
ylab("Richness")+xlab(("Depth (m)"))+ scale_x_reverse()
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 = T, message=FALSE, warning=FALSE, dpi=300)
devtools::load_all(path = ".") #Load packages library(tidyverse) library(knitr) library(ggpubr) library(data.table) library(kableExtra) library(knitr) library(vegan) library(gamm4) library(voxel) library(gridExtra) library(ggpmisc) library(grid) ## ggplot theme source(paste0(here::here(),"/R/theme.R"))
Material & methods
Data overview
Abiotic data
Chem_overview_SmallDatasetMorph<-Chem_uniform_LOIx %>% filter(dataset=="SAK") %>% filter(datasetWLF=="WLF") %>% group_by(Lake)%>% select(Area_ha, maxDepth_m, WLF) %>% unique() %>% ungroup() %>% rename_(.dots=setNames(names(.), tolower(gsub("\\_", ".", names(.)))))%>% summarise_each(funs(mean, median, min, max, sd),-lake) %>% tidyr::gather(variable, value) %>% tidyr::separate(variable, c("var", "stat"), sep = "\\_") %>% tidyr::spread(var, value) Chem_overview_SmallDataset<-Chem_uniform_LOIx %>% rowwise() %>% filter(dataset=="SAK") %>% filter(datasetWLF=="WLF") %>% rename_(.dots=setNames(names(.), tolower(gsub("\\_", ".", names(.)))))%>%ungroup()%>% select(chloride,conduct,nh4n,no3n,ntot,o2diss,ph,ptot,sac,sio2,temp,tempsd,transp,wlf)%>% summarise_each(funs(mean, median, min, max, sd))%>% tidyr::gather(variable, value) %>% tidyr::separate(variable, c("var", "stat"), sep = "\\_") %>% tidyr::spread(var, value)#%>%t() Chem_overview_SmallD<-cbind(Chem_overview_SmallDatasetMorph,Chem_overview_SmallDataset[,2:15]) %>% tibble::rownames_to_column() %>% tidyr::gather(var, value, -rowname) %>% tidyr::spread(rowname, value) %>% rename(max="1",mean="2", median="3", min="4", sd="5") df<-cbind(Chem_overview_SmallDatasetMorph,Chem_overview_SmallDataset[,2:15]) %>% column_to_rownames(var = "stat") dt<-type_convert(as_tibble(cbind(Parameters = names(df), t(df)))) %>% mutate_if(is.numeric, round, 2) dt %>% select(Parameters,min, max,mean,sd,median)
Biotic data, dataset levels overview
complete<-MakrophS_ALL %>% dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N =n_distinct(Probestelle),Years_N = n_distinct(YEAR)) Level1<-MakrophS_ALL %>% filter(datasettot=="LEVEL3") %>% #group_by(datasettot) %>% dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N = n_distinct(Probestelle),Years_N = n_distinct(YEAR)) Level_Timeseries<-MakrophS_ALL %>% group_by(Lake)%>%filter(n_distinct(YEAR)>3)%>% ungroup() %>% dplyr::summarise(N = n_distinct(Lake,YEAR),Lakes_N = n_distinct(Lake),Transects_N = n_distinct(MST_NR),Depths_N = n_distinct(Probestelle),Years_N = n_distinct(YEAR)) Level<- rbind(complete,Level1,Level_Timeseries) Level$dataset<-c("Biotic","Biotic+Abiotic", "Timeseries") Level[,c(6,1,2,3,4,5)]
Results
Species richness per dataset group
M<-MakrophS_ALL[5:96] M[colSums(M)!=0] %>% ncol() #75 M3<-MakrophS_ALL %>% filter(datasettot=="LEVEL3") M3<-M3[5:96] M3[colSums(M3)!=0] %>% ncol() #57 MT<-MakrophS_ALL %>% group_by(Lake)%>%filter(n_distinct(YEAR)>3)%>% ungroup() MT<-MT[5:96] MT[colSums(MT)!=0] %>% ncol() #66
Species numbers in taxonomic groups
M<-MakrophS_ALL[5:96] specpres<-M[colSums(M)!=0] species <- as.data.frame(colnames(specpres)) species names(species)<-"Taxon" unique(left_join(species,Makroph[c(8,11)],by=c("Taxon"))) %>% group_by(Erscheinungsform) %>% dplyr::summarise(N=n_distinct(Taxon))
Gamma richness overview: min, max, mean, sd
min(Makroph_Lake_ALL$GAMMA) max(Makroph_Lake_ALL$GAMMA) mean(Makroph_Lake_ALL$GAMMA) sd(Makroph_Lake_ALL$GAMMA)
Depth diversity gradients
A1<-ggplot(data=Makroph_Depth)+ geom_ribbon(aes(x = Tiefe, ymax =mAlpha+sdAlpha, ymin =mAlpha-sdAlpha), alpha = 0.6, fill = "grey")+ geom_line(aes(x=(Tiefe), y=mAlpha), size=1)+ scale_x_reverse()+ ylab("Alpha richness")+xlab("Depth (m)") formula1 = y ~ poly(x, 3, raw=TRUE) A2<-ggplot(data = PEAK, aes(x=AlphaPeakDepth, y=AlphaPeakRichness))+ geom_point(aes(shape=datasettotsimpl),alpha=0.2)+ scale_shape_discrete(name = "s", labels = c("DDG peak, Biodiversity dataset", "DDG peak, Environmental & biodiversity dataset"))+ geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+ xlab(expression(D[alpha][max](m)))+ ylab(expression(R[alpha][max](N)))+ theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+ scale_x_reverse()+ xlim(-0.5,-5) A3<-ggplot(data=PEAK, aes(x=AlphaPeakDepth))+ geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1 scale_x_reverse()+xlab("Depth (m)")+ylab(expression(D[alpha][max](counts))) B1<-ggplot(data=Makroph_Depth)+ geom_ribbon(aes(x = Tiefe, ymax =mBeta+sdBeta, ymin =mBeta-sdBeta), alpha = 0.6, fill = "grey")+ geom_line(aes(x=(Tiefe), y=mBeta), size=1)+ scale_x_reverse()+ ylab("Beta richness")+xlab(("Depth (m)")) B2<-ggplot(data = PEAK, aes(x=BetaPeakDepth, y=BetaPeakRichness))+ geom_point(aes(shape=datasettotsimpl),alpha=0.2)+ scale_shape_discrete(name = "s", labels = c("DDG peak, Biodiversity dataset", "DDG peak, Environmental & biodiversity dataset"))+ geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+ xlab(expression(D[beta][max](m)))+ ylab(expression(R[beta][max](N)))+ theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+ scale_x_reverse()+ xlim(-0.5,-5) B3<-ggplot(data=PEAK, aes(x=BetaPeakDepth))+ geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1 scale_x_reverse()+xlab("Depth (m)")+ ylab(expression(D[beta][max](counts))) C1<-ggplot(data=Makroph_Depth)+ geom_ribbon(aes(x = Tiefe, ymax =mGamma+sdGamma, ymin =mGamma-sdGamma), alpha = 0.6, fill = "grey")+ geom_line(aes(x=(Tiefe), y=mGamma), size=1, colour="black")+ scale_x_reverse()+ ylab("Gamma richness")+xlab(("Depth (m)")) C2<-ggplot(data = PEAK, aes(x=GammaPeakDepth, y=GammaPeakRichness))+ geom_point(aes(shape=datasettotsimpl),alpha=0.2)+ scale_shape_discrete(name = "s", labels = c("DDG peak, Biodiversity dataset", "DDG peak, Environmental & biodiversity dataset"))+ geom_smooth(method=lm,model=lm, formula = formula1,span=1, col="black")+ xlab(expression(D[gamma][max](m)))+ ylab(expression(R[gamma][max](N)))+ theme(legend.position = c(0.9, 0.2),legend.title=element_blank())+ scale_x_reverse()+ xlim(-0.5,-5) C3<-ggplot(data=PEAK, aes(x=GammaPeakDepth))+ geom_histogram(breaks=c(0,-1,-2,-4,-5), color="black",fill="white")+#binwidth=1 scale_x_reverse()+xlab("Depth (m)")+ylab(expression(D[gamma][max](counts))) Fig2<-ggarrange(A1,B1,C1,A2,B2,C2,A3,B3,C3, ncol=3,nrow=3,labels=c("(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)"), heights = c(2,1.5,1.8),common.legend = T, legend="none", align = "hv") # png("Fig2.png", width=16.6, height=18,units = "cm",res=300) # Fig2 # dev.off # # ggsave("Fig2.pdf",plot=Fig1,width=16.6, height = 18, device=cairo_pdf, units = "cm") windowsFonts("Arial" = windowsFont("Arial")) Fig2
Drivers
Cor.test
cor.test(PEAK_Chem_norm$AlphaPeakRichness, PEAK_Chem_norm$Ntot)
GAMM incl PCA
################ PCA ############################################### label_lake<- PEAK_Chem_norm[complete.cases(PEAK_Chem_norm[,c(3:16,19)]),] #LEVEL 3 data 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")) Rotation <-lak.pca$rotation %>% as.data.frame() Rotation$variable <- row.names(Rotation) R1<-ggplot(data=Rotation)+ geom_bar(aes(y=variable,x=PC1,fill = PC1 > 0.4 | PC1< -0.4), stat='identity')+ xlim(-0.65,0.65)+ #theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ ylab("")+xlab("loading")+ ggtitle(" PC1 - 30.1% \nSiO2 & Cond axis") R2<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC2,fill = PC2 > 0.4 | PC2< -0.4), stat='identity')+ xlim(-0.65,0.65)+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ylab("")+ ggtitle(" PC2 - 26.1%\nTemp & Ptot axis")+xlab("loading")+ylab("") R3<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC3,fill = PC3 > 0.4 | PC3< -0.4), stat='identity')+ xlim(-0.65,0.65)+ylab("")+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ ggtitle(" PC3 - 13.3%\nTempsd – Chl axis")+xlab("loading") R4<-ggplot(data=Rotation)+geom_bar(aes(y=variable,x=PC4,fill = PC4 > 0.4 | PC4< -0.4), stat='identity')+ xlim(-0.65,0.65)+ylab("")+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ ggtitle(" PC4 - 10.5%\nO2diss – SAC axis")+xlab("loading") ## GAMM gam_AlphaPeakDepth <- gamm4(AlphaPeakDepth ~ s(PC1)+s(PC2)+s(PC3)+s(PC4), random= ~(1|Lake), # package gamm4 data=PEAK_PCA) summary(gam_AlphaPeakDepth$gam) gam_AlphaPeakRichness <- gamm4(AlphaPeakRichness ~ s(PC1), random= ~(1|Lake), # package gamm4 data=PEAK_PCA) summary(gam_AlphaPeakRichness$gam) vars <- c("PC2", "PC4","PC3","PC1") grob2 <- grobTree(textGrob("*** 36.6% dc", x=0.1, y=0.9, hjust=0, gp=gpar(col="grey12", fontsize=13, fontface="italic"))) plotPC2<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC2") + geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC2"), alpha = 0.2) + geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC2"), alpha = 0.2) + #scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) + theme(legend.position="none")+ ylab("")+ylim(-3.5,-0.5)+ ggtitle("")+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ theme(legend.title = element_blank())+ annotation_custom(grob2) grob4 <- grobTree(textGrob("** 30.3% dc", x=0.1, y=0.9, hjust=0, gp=gpar(col="grey12", fontsize=13, fontface="italic"))) vars <- c("PC4") plotPC4<- plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC4") + geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC4"), alpha = 0.2) + geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC4"), alpha = 0.2) + #scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) + theme(legend.position="none")+theme(legend.title = element_blank())+ ylab("")+ylim(-3.5,-0.5)+ #ylab(expression(predicted_D[alpha][max]))+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ ggtitle("")+ annotation_custom(grob4) grob3 <- grobTree(textGrob("** 28.6% dc", x=0.1, y=0.9, hjust=0, gp=gpar(col="grey12", fontsize=13, fontface="italic"))) vars <- c("PC3") plotPC3<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC3") + geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC3"), alpha = 0.2) + geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC3"), alpha = 0.2) + #scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) + theme(legend.position="none")+theme(legend.title = element_blank())+ ylab("")+ylim(-3.5,-0.5)+ #ylab(expression(predicted_D[alpha][max]))+ theme(axis.title.y=element_blank(),axis.text.y=element_blank())+ ggtitle("")+ annotation_custom(grob3) grob1 <- grobTree(textGrob("* 25.6% dc", x=0.1, y=0.9, hjust=0, gp=gpar(col="grey12", fontsize=13, fontface="italic"))) vars <- c("PC1") plotPC1<-plotGAMM(gam_AlphaPeakDepth, smooth.cov = "PC1") + geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC1"), alpha = 0.2) + geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakDepth", x = "PC1"), alpha = 0.2) + #scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) + theme(legend.position="none")+ ylab(expression(predicted_D[alpha][max]))+ylim(-3.5,-0.5)+ ggtitle("")+ theme(legend.title = element_blank())+ annotation_custom(grob1) grob5 <- grobTree(textGrob("**", x=0.1, y=0.9, hjust=0, gp=gpar(col="grey12", fontsize=13, fontface="italic"))) vars2 <- c("PC1") plot2<-plotGAMM(gam_AlphaPeakRichness, smooth.cov = "PC1") + geom_point(data = PEAK_PCA, aes_string(y = "AlphaPeakRichness", x = "PC1"), alpha = 0.2) + geom_rug(data = PEAK_PCA, aes_string(y = "AlphaPeakRichness", x = "PC1"), alpha = 0.2) + #scale_color_manual("Private", values = c("#868686FF", "#0073C2FF")) + theme(legend.position="none")+ ylab(expression(predicted_R[alpha][max]))+ #labs(tag = "E")+ ggtitle("")+ theme(legend.title = element_blank())+ annotation_custom(grob5) empty<-plot(0, xaxt = 'n', yaxt = 'n', bty = 'n', pch = '', ylab = '', xlab = '') #ggarrange(R1,R2,R3,R4, # ncol=4, common.legend = T, legend="top", labels = c("(a)","(b)","(c)","(d)"), widths = c(1.3,1,1,1)) # ggarrange(plotPC1, plotPC2, plotPC3, plotPC4, plot2, # ncol=4, nrow=3, # heights = c(4,3,3), # common.legend = T, legend="top", # labels = c("(e)","(f)","(g)","(h)","(i)"), # #widths = c(1.4,1,1,1), # align = "hv" # ) Fig3<-ggarrange(R1,R2,R3,R4, plotPC1, plotPC2, plotPC3, plotPC4, plot2,empty, empty, empty, ncol=4, nrow=3, heights = c(4,3,3), common.legend = T, legend="top", labels = c("(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)","(i)"), #widths = c(1.4,1,1,1), align = "hv" ) Fig3 # ggsave("Fig3.pdf",plot=Fig3,width = 9, height = 8, device=cairo_pdf) # # png("Fig3.png", width=16.6, height=18,units = "cm",res=300) # Fig3 # dev.off
Temporal change
Invariability analysis
LAKECHANGEInv<-PEAK %>% dplyr::group_by(Lake) %>% #summarize informtion for lakes (over timeseries) dplyr::summarise(NYEAR=n_distinct(YEAR), AlphaPeakDepthInv=mean(-AlphaPeakDepth)/sd(-AlphaPeakDepth), AlphaPeakRichnessInv=mean(AlphaPeakRichness)/sd(AlphaPeakRichness), BetaPeakDepthInv=mean(-BetaPeakDepth)/sd(-BetaPeakDepth), BetaPeakRichnessInv=mean(BetaPeakRichness)/sd(BetaPeakRichness), GammaPeakDepthInv=mean(-GammaPeakDepth)/sd(-GammaPeakDepth), GammaPeakRichnessInv=mean(GammaPeakRichness)/sd(GammaPeakRichness), GammaRichnessInv=mean(GAMMA)/sd(GAMMA) )%>% filter(NYEAR>3) %>% #For timeseries dataset arrange(AlphaPeakDepthInv) %>% filter_all(all_vars(!is.infinite(.))) #LAKECHANGEInv # ggplot(LAKECHANGEInv %>% gather(Type, "Inv", 3:9)) + # geom_histogram(aes(x=Inv), col="white", fill="black",binwidth = 1)+ # facet_wrap(~Type, ncol = 2)+ # xlab("Inverse Variation coefficient = Invariability coefficient") LAKECHANGEInv%>% dplyr::ungroup() %>% select(-Lake, -NYEAR)%>% dplyr::summarise_all(list(mean=mean,sd=sd), na.rm=TRUE)
Linear models
PEAK<-PEAK %>%mutate(YEAR=as.numeric(YEAR)) lmp <- function (modelobject) { if (class(modelobject) != "lm") stop("Not an object of class 'lm' ") f <- summary(modelobject)$fstatistic p <- pf(f[1],f[2],f[3],lower.tail=FALSE) attributes(p) <- NULL return(p) } GammaModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(GAMMA ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="GammaSlope", "pValue"="GammaPValue")) AlphaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(AlphaPeakRichness ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="AlphaPeakRichnessSlope", "pValue"="AlphaPeakRichnessPValue")) AlphaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(AlphaPeakDepth ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="AlphaPeakDepthSlope", "pValue"="AlphaPeakDepthPValue")) AlphaPeakModel <- merge (AlphaPeakDepthModel, AlphaPeakRichnessModel, by="Lake") BetaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(BetaPeakRichness ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="BetaPeakRichnessSlope", "pValue"="BetaPeakRichnessPValue")) BetaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(BetaPeakDepth ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="BetaPeakDepthSlope", "pValue"="BetaPeakDepthPValue")) BetaPeakModel <- merge (BetaPeakDepthModel, BetaPeakRichnessModel, by="Lake") GammaPeakRichnessModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(GammaPeakRichness ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="GammaPeakRichnessSlope", "pValue"="GammaPeakRichnessPValue")) GammaPeakDepthModel<-PEAK%>%group_by(Lake)%>%filter(n_distinct(YEAR)>3) %>% do({ mod = lm(GammaPeakDepth ~ YEAR, data = .) data.frame(Slope = coef(mod)[2],pValue=lmp(mod)) }) %>% plyr::rename(c("Slope"="GammaPeakDepthSlope", "pValue"="GammaPeakDepthPValue")) GammaPeakModel <- merge (GammaPeakDepthModel, GammaPeakRichnessModel, by="Lake") PeakModel <- merge (AlphaPeakModel,BetaPeakModel,by="Lake") PeakModel <- merge (PeakModel,GammaPeakModel,by="Lake") Model <- merge (GammaModel,PeakModel,by="Lake") ModelT<-Model %>% mutate(GammaTrend=ifelse(Model$GammaSlope>0,"+",ifelse(Model$GammaSlope<0,"-","0"))) %>% mutate(GammapV=ifelse(GammaPValue<0.001,"***",ifelse(GammaPValue<0.01,"**",ifelse(GammaPValue<0.05,"*",ifelse(GammaPValue<0.1,".","NA")))))%>% mutate(AlphaPeakDepthTrend=ifelse(Model$AlphaPeakDepthSlope>0.5*0,"+",ifelse(Model$AlphaPeakDepthSlope<0,"-","0"))) %>% mutate(AlphaPeakDepthpV=ifelse(AlphaPeakDepthPValue<0.001,"***",ifelse(AlphaPeakDepthPValue<0.01,"**",ifelse(AlphaPeakDepthPValue<0.05,"*",ifelse(AlphaPeakDepthPValue<0.1,".","NA")))))%>% mutate(AlphaPeakRichnessTrend=ifelse(Model$AlphaPeakRichnessSlope>0,"+",ifelse(Model$AlphaPeakRichnessSlope<0,"-","0")))%>% mutate(AlphaPeakRichnesspV=ifelse(AlphaPeakRichnessPValue<0.001,"***",ifelse(AlphaPeakRichnessPValue<0.01,"**",ifelse(AlphaPeakRichnessPValue<0.05,"*",ifelse(AlphaPeakRichnessPValue<0.1,".","NA")))))%>% mutate(BetaPeakDepthTrend=ifelse(Model$BetaPeakDepthSlope>0,"+",ifelse(Model$BetaPeakDepthSlope<0,"-","0"))) %>% mutate(BetaPeakDepthpV=ifelse(BetaPeakDepthPValue<0.001,"***",ifelse(BetaPeakDepthPValue<0.01,"**",ifelse(BetaPeakDepthPValue<0.05,"*",ifelse(BetaPeakDepthPValue<0.1,".","NA")))))%>% mutate(BetaPeakRichnessTrend=ifelse(Model$BetaPeakRichnessSlope>0,"+",ifelse(Model$BetaPeakRichnessSlope<0,"-","0")))%>% mutate(BetaPeakRichnesspV=ifelse(BetaPeakRichnessPValue<0.001,"***",ifelse(BetaPeakRichnessPValue<0.01,"**",ifelse(BetaPeakRichnessPValue<0.05,"*",ifelse(BetaPeakRichnessPValue<0.1,".","NA")))))%>% mutate(GammaPeakDepthTrend=ifelse(Model$GammaPeakDepthSlope>0,"+",ifelse(Model$GammaPeakDepthSlope<0,"-","0"))) %>% mutate(GammaPeakDepthpV=ifelse(GammaPeakDepthPValue<0.001,"***",ifelse(GammaPeakDepthPValue<0.01,"**",ifelse(GammaPeakDepthPValue<0.05,"*",ifelse(GammaPeakDepthPValue<0.1,".","NA")))))%>% mutate(GammaPeakRichnessTrend=ifelse(Model$GammaPeakRichnessSlope>0,"+",ifelse(Model$GammaPeakRichnessSlope<0,"-","0")))%>% mutate(GammaPeakRichnesspV=ifelse(GammaPeakRichnessPValue<0.001,"***",ifelse(GammaPeakRichnessPValue<0.01,"**",ifelse(GammaPeakRichnessPValue<0.05,"*",ifelse(GammaPeakRichnessPValue<0.1,".","NA"))))) ModelTshort<-ModelT[c(1,16:29)] %>% tidyr::unite(GammaLM, c("GammaTrend", "GammapV"))%>% tidyr::unite(AlphaPeakDepthLM, c("AlphaPeakDepthTrend", "AlphaPeakDepthpV"))%>% tidyr::unite(AlphaPeakRichnessLM, c("AlphaPeakRichnessTrend", "AlphaPeakRichnesspV"))%>% tidyr::unite(BetaPeakDepthLM, c("BetaPeakDepthTrend", "BetaPeakDepthpV"))%>% tidyr::unite(BetaPeakRichnessLM, c("BetaPeakRichnessTrend", "BetaPeakRichnesspV"))%>% tidyr::unite(GammaPeakDepthLM, c("GammaPeakDepthTrend", "GammaPeakDepthpV"))%>% tidyr::unite(GammaPeakRichnessLM, c("GammaPeakRichnessTrend", "GammaPeakRichnesspV")) ModelTshort <- data.frame(lapply(ModelTshort, function(x) {gsub("_NA", "", x)})) ModelTshort <- data.frame(lapply(ModelTshort, function(x) {gsub("_", "", x)})) ModelTshort
Discussion
Summary figure
ggplot(data=Makroph_Depth)+ geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mAlpha), col="darkgreen", size=1)+ geom_line(data=Makroph_Depth, aes(x=(Tiefe), y=mBeta), col="green2", size=1)+ geom_line(data=Makroph_Depth,aes(x=(Tiefe), y=mGamma), col="lightgreen", size=1, colour="black")+ ylab("Richness")+xlab(("Depth (m)"))+ scale_x_reverse()
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.