---
# Script for: "Cryptogam plant community stability: warming weakens influences of species richness but enhances effects of evenness"
 # Qian Gu et al. Ecology
---

Load packages  
```{r}
rm(list=ls())
# install.packages("") # Install if needed
# Load libraries
library(tidyverse)
library(lavaan)
library(patchwork)
library(nlme)
library(MuMIn)
library(ggeffects)

```


Load data
```{r}
# .csv files must be in current working directory
lichen_moss <- read.csv("organized cryptogam cover data.csv",header=TRUE) 
all.metrics <- read.csv("calculated metrics.csv",header=TRUE)
vascular <-  read.csv('vascular abundance.csv',header = TRUE)
mydata <- read.csv("biomass-hits data.csv",header=TRUE)
vas_nonvas <- read.csv('vascular abundance and cryptogam richness and evenness.csv',header = TRUE) 

```


Table 1
```{r}
spp <- lichen_moss %>% 
  group_by(year, trt, category, plot, species) %>% 
  summarise(abundance=sum(abundance)) %>% 
  ungroup() %>% 
   spread(species, abundance)
spp[is.na(spp)] <-  0

spp <- spp %>% 
  gather(key='species',value='abundance',c(5:ncol(spp)))

# Table 1 results
con_warm <- spp %>% 
  mutate(year=factor(year),trt=factor(trt),plot=factor(plot),species=factor(species)) %>% 
  group_by(trt,plot,species) %>%
  summarise(variance=sd(abundance),mean=mean(abundance)) %>%
  ungroup() %>% 
  arrange(species) %>% 
  filter(variance!=0) %>% 
  group_by(trt,species) %>% 
  summarise(abundance=mean(mean)) %>% 
  pivot_wider(names_from = trt,values_from=abundance) %>% 
  replace(is.na(.),0) %>% 
  rename(abundance_C=Control,abundance_W=Warming) %>% 
  mutate(dominance_C=100*abundance_C/sum(abundance_C),dominance_W=100*abundance_W/sum(abundance_W)) %>% 
  select(species,dominance_C,dominance_W) %>% 
    arrange(desc(dominance_C)) %>% 
  mutate(across(where(is.numeric), round, 2))

```


Fig 1.
```{r}
### Fig 1. Statistics
all.metrics <- all.metrics %>% 
  mutate(trt=case_when(
    trt=='WC'~'Control',
    trt=='W'~'Warming'
    )) 
  
###
# Fig.1(a): Warming effect on community stability
myfit <- lm(log(stability)~trt,data=all.metrics)
summary(myfit) 
# Warming decreased community stability.
all.metrics %>% 
  group_by(trt) %>% 
  summarise(mean_stability=mean(stability))
(3.11-1.56)/3.11

# Warming effect on community mean
myfit <- lm(log(community_mean)~trt,data=all.metrics)
summary(myfit) 
# Warming decreased community mean.
all.metrics %>% 
  group_by(trt) %>% 
  summarise(mean=mean(community_mean))
(127-87.1)/127

# Warming effect on community sd
myfit <- lm(log(community_sd)~trt,data=all.metrics)
summary(myfit) 
# Warming had no effect on community sd.


###
# Fig.1(b): Warming effect on species richness
myfit <- lm(log(richness_mean)~trt, data=all.metrics)
summary(myfit) 
# Warming had no effect on species richness.


###
# Fig.1(c): Warming effect on species evenness
myfit <- lm(log(Evar)~trt, data=all.metrics)
summary(myfit) 
# Warming increased species evenness.
all.metrics %>% 
  group_by(trt) %>% 
  summarise(mean_even=mean(Evar))
(0.39-0.487)/0.39
  

###
# Fig.1(d): Warming effect on synchrony
myfit <- lm(log(synchrony)~trt,data=all.metrics)
summary(myfit) 
# Warming increased synchrony.
all.metrics %>% 
  group_by(trt) %>% 
  summarise(mean_synchrony=mean(synchrony))
(0.181-0.356)/0.181



### Fig 1. Graphs
global_theme <- theme_bw()+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text=element_text(size=50,colour='black'))+
  theme(axis.title=element_text(size=50))+
  theme(axis.ticks.length = unit(10,'pt'))+
  theme(legend.text = element_text(size = 50))+
  theme(legend.title = element_text(size = 50))+
  theme(legend.position = 'none')+
  theme(plot.margin=unit(c(1,1.5,1,1),"cm"))+
  theme(text=element_text(family="Times New Roman")) 

fig1.a <-  ggplot(all.metrics, aes(x=trt, y=stability,fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
   labs(x="", y="Community stability\n") +
  global_theme+
  scale_fill_manual(values=c('gray','red'))+
  scale_y_continuous(limits = c(0,10.5),breaks = c(2,4,6,8,10),expand=c(0,0))

fig1.b <- ggplot(all.metrics, aes(x=trt, y=richness_mean,fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
   labs(x="", y="Species richness\n") +
  global_theme+
  scale_fill_manual(values=c('gray','red'))+
  scale_y_continuous(limits = c(3.5,12.5),breaks = c(4,6,8,10,12),expand=c(0,0))

fig1.c <- ggplot(all.metrics, aes(x=trt, y=Evar,fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
   labs(x="", y="Species evenness\n") +
  global_theme+
  scale_fill_manual(values=c('gray','red'))+
  scale_y_continuous(limits = c(0.35,0.85),breaks = c(0.4,0.5,0.6,0.7,0.8),expand=c(0,0))

fig1.d <- ggplot(all.metrics, aes(x=trt, y=synchrony,fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
   labs(x="", y="Species synchrony\n") +
  global_theme+
  scale_fill_manual(values=c('gray','red'))+
  scale_y_continuous(limits = c(0,0.75),breaks = c(0.1,0.3,0.5,0.7),expand=c(0,0))

# Combine all panels together:
fig1 <- 
(fig1.a|fig1.c)/
(fig1.b|fig1.d)

```


Fig 2
```{r}
### Fig 2 statistics
### Fig.2(a)
# Relationship between community stability and species richness
myfit <- lm(log(stability)~log(richness_mean)*trt,data=all.metrics)
plot(myfit)
summary(myfit) 
# Significant interactive effect between richness and treatment on community stability.
# (1) Control 
myfit <- lm(log(stability)~log(richness_mean),data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# Positive correlation.
# Perform the Shapiro-Wilk test of normality:
sresid <-  MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant. 
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit) 
# Not significant.
#(2) Warming
myfit <- lm(log(stability)~log(richness_mean),data=subset(all.metrics,trt=='Warming'))
summary(myfit) 
# No correlation.
# Perform the Shapiro-Wilk test of normality:
sresid <-  MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit) 
# Not significant


### Fig.2(b)
# Relationship between community mean and richness
myfit <- lm(log(community_mean)~log(richness_mean)*trt, data=all.metrics)
summary(myfit) 
# No significant interactive effect between species richness and treatment on community mean.
# (1) Control 
myfit <- lm(log(community_mean)~log(richness_mean), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# No correlation.
# (2) Warming
myfit <- lm(log(community_mean)~log(richness_mean), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
#No correlation.


### Fig.2(c)
# Relationship between community SD and richness
myfit <- lm(log(community_sd)~log(richness_mean)*trt, data=all.metrics)
summary(myfit) 
# No significant interactive effect between species richness and treatment on community sd.
# (1) Control 
myfit <- lm(log(community_sd)~log(richness_mean), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# No correlation.
# (2) Warming
myfit <- lm(log(community_sd)~log(richness_mean), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
# No correlation.


### Fig.2(d)
# Relationship between community stability and species evenness
myfit <- lm(log(stability)~log(Evar)*trt, data=na.omit(all.metrics))
summary(myfit) 
# No significant interactive effect between species evenness and treatment on community stability.
# (1) Control
myfit <- lm(log(stability)~log(Evar), data=na.omit(subset(all.metrics,trt=='Control')))
summary(myfit) 
# No correlation.
# Perform the Shapiro-Wilk test of normality:
sresid <-  MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant 
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit)
# Not significant
# (2) Warming
myfit <- lm(log(stability)~log(Evar), data=na.omit(subset(all.metrics,trt=='Warming')))
summary(myfit) 
# Negative correlation.
# Perform the Shapiro-Wilk test of normality:
sresid <-  MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit) 
# Not significant


### Fig.2(e)
# Relationship between community mean and evenness
myfit <- lm(log(community_mean)~log(Evar)*trt, data=all.metrics)
summary(myfit) 
# No significant interactive effect.
# (1) Control 
myfit <- lm(log(community_mean)~log(Evar), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# Negative correlation.
# (2) Warming
myfit <- lm(log(community_mean)~log(Evar), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
# Negative correlation.


### Fig.2(f)
# Relationship between community SD and evenness
myfit <- lm(log(community_sd)~log(Evar)*trt, data=all.metrics)
summary(myfit) 
# No significant interactive effect.
# (1) Control 
myfit <- lm(log(community_sd)~log(Evar), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# No correlation.
# (2) Warming
myfit <- lm(log(community_sd)~log(Evar), data=subset(all.metrics,trt=='Warming'))
summary(myfit) 
# No correlation.


### Fig.2(g)
# Relationship between community stability and synchrony 
myfit <- lm(log(stability)~log(synchrony)*trt, data=all.metrics)
summary(myfit) 
# No significant interactive effect.
# (1) Control 
myfit <- lm(log(stability)~log(synchrony), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# Negative relationship. 
# Perform the Shapiro-Wilk test of normality:
sresid <-  MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit) 
# Not significant
# (2) Warming
myfit <- lm(log(stability)~log(synchrony), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
# Negative relationship.
# Perform the Shapiro-Wilk test of normality:
sresid <- MASS::studres(myfit)
shapiro.test(sresid) 
# Not significant
# Perform Breusch-Pagan test of homoscedasticity:
lmtest::bptest(myfit) 
# Not significant


### Fig.2(h)
# Relationship between community mean and synchrony
myfit <- lm(log(community_mean)~log(synchrony)*trt, data=all.metrics)
summary(myfit) 
# Negative interactive effect.
# (1) Control 
myfit <- lm(log(community_mean)~log(synchrony), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# No correlation.
# (2) Warming
myfit <- lm(log(community_mean)~log(synchrony), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
# Negative correlation.


### Fig.2(i)
# Relationship between community SD and synchrony
myfit <- lm(log(community_sd)~log(synchrony)*trt, data=all.metrics)
summary(myfit) 
# Negative interactive effect.
# (1) Control 
myfit <- lm(log(community_sd)~log(synchrony), data=subset(all.metrics,trt=='Control'))
summary(myfit) 
# Positive correlation.
# (2) Warming
myfit <- lm(log(community_sd)~log(synchrony), data=subset(all.metrics,trt=='Warming'))
summary(myfit)
# No correlation.




### Fig 2. Graphs
global_theme <- theme_bw()+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text=element_text(size=60,colour='black'))+
  theme(axis.title=element_text(size=60))+
  theme(axis.ticks.length = unit(10,'pt'))+
  theme(legend.text = element_text(size = 60))+
  theme(legend.title = element_text(size = 60))+
  theme(legend.position = 'none')+
  theme(plot.margin=unit(c(1,1.5,1,1),"cm"))+
  theme(text=element_text(family="Times New Roman")) 

# Set x,y scale used in Fig.2
stability_y <- scale_y_continuous(limits = c(-0.4,2.2),breaks = c(0,0.5,1,1.5,2), expand=c(0,0))
mean_y <- scale_y_continuous(limits = c(3.3,5.7),breaks = c(3.5,4,4.5,5,5.5), expand=c(0,0))
sd_y <- scale_y_continuous(limits = c(2.6,4.7),breaks = c(3,3.5,4,4.5), expand=c(0,0))
richness_x <- scale_x_continuous(limits = c(1.4, 2.6),breaks = c(1.5,1.7,1.9,2.1,2.3,2.5), expand = c(0,0))
evenness_x <- scale_x_continuous(limits = c(-1.3,-0.3),breaks = c(-1.2,-0.8,-0.4), expand = c(0,0))
synchrony_x <- scale_x_continuous(limits = c(-4.1,0.1),breaks = c(-4,-3,-2,-1,0), expand = c(0,0))


### Fig 2(a)
# Relationships between community stability and species richness
stability_richness <- ggplot(all.metrics, aes(x=log(richness_mean), y=log(stability))) + 
  geom_point(aes(color=trt),size=10) +
  stat_smooth(method='lm',color='black',se=F,size=2,data=subset(all.metrics,trt=='Control')) +
   labs(x="", y="\nln(Community stability)\n") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  richness_x + 
  stability_y


### Fig 2(b)
# Relationships between community mean and richness
mean_richness <- ggplot(all.metrics, aes(x=log(richness_mean), y=log(community_mean))) + 
  geom_point(aes(color=trt),size=10) +
   labs(x="", y="\nln(Community mean)\n") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  richness_x + 
  mean_y


### Fig 2(c)
# Relationships between community SD and richness
sd_richness <- ggplot(all.metrics, aes(x=log(richness_mean), y=log(community_sd))) + 
  geom_point(aes(color=trt),size=10) +
   labs(x="\nln(Species richness)\n", y="\nln(Community SD)\n") +
  scale_color_manual(values=c('black','red'))+
  global_theme +
  richness_x + 
  sd_y


### Fig 2(d)
# Relationships between community stability and species evenness
stability_evenness <- ggplot(all.metrics, aes(x=log(Evar), y=log(stability))) + 
  geom_point(aes(color=trt),size=10) +
  stat_smooth(method='lm',color='red',se=F,size=2,data=subset(all.metrics,trt=='Warming'))+
   labs(x="", y="") +
  scale_color_manual(values=c('black','red'))+
  global_theme +
  evenness_x + 
  stability_y


### Fig 2(e)
# Relationships between community mean and evenness
mean_evenness <- ggplot(all.metrics, aes(x=log(Evar), y=log(community_mean))) + 
  geom_point(aes(color=trt),size=10) +
 stat_smooth(method='lm',color='black',se=F,size=2,data=subset(all.metrics,trt=='Control')) +
  stat_smooth(method='lm',color='red',se=F,size=2,data=subset(all.metrics,trt=='Warming')) +
   labs(x="", y="") +
  scale_color_manual(values=c('black','red'))+
  global_theme +
  evenness_x + 
  mean_y


### Fig 2(f)
# Relationships between community SD and evenness
sd_evenness <- ggplot(all.metrics, aes(x=log(Evar), y=log(community_sd))) + 
  geom_point(aes(color=trt),size=10) +
   labs(x="\nln(Species evenness)\n", y="") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  evenness_x + 
  sd_y


### Fig 2(g)
# Relationships between community stability and synchrony 
stability_synchrony <- ggplot(all.metrics, aes(x=log(synchrony), y=log(stability))) +
  geom_point(aes(color=trt),size=10) +
  stat_smooth(method='lm',color='black',se=F,size=2,data=subset(all.metrics,trt=='Control')) +
  stat_smooth(method='lm',color='red',se=F,size=2,data=subset(all.metrics,trt=='Warming')) +
   labs(x="", y="") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  synchrony_x +
  stability_y


### Fig 2(h)
# Relationships between community mean and synchrony
mean_synchrony <-  ggplot(all.metrics, aes(x=log(synchrony), y=log(community_mean))) + 
  geom_point(aes(color=trt),size=10) +
  stat_smooth(method='lm',color='red',se=F,size=2,data=subset(all.metrics,trt=='Warming'))+
   labs(x="", y="") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  synchrony_x +
  mean_y


### Fig 2(i) 
# Relationships between community SD and synchrony
sd_synchrony <-  ggplot(all.metrics, aes(x=log(synchrony), y=log(community_sd))) + 
  geom_point(aes(color=trt),size=10) +
  stat_smooth(method='lm',color='black',se=F,size=2,data=subset(all.metrics,trt=='Control'))+
   labs(x="\nln(Species synchrony)\n", y="") +
  scale_color_manual(values=c('black','red')) +
  global_theme +
  synchrony_x +
  sd_y


# Combine all panels together
fig2 <- 
(stability_richness|stability_evenness|stability_synchrony)/
(mean_richness|mean_evenness|mean_synchrony)/
(sd_richness|sd_evenness|sd_synchrony)

```


Fig 3, Fig S5 (Appendix), and Fig S6 (Appendix).
```{r}
semdata <- inner_join(all.metrics,vascular, by=c('trt','plot'))
semdata <- semdata %>% 
  mutate(ln.vas=log(vas_abundance), ln.syn=log(synchrony),ln.richness=log(richness_mean), ln.Evar=log(Evar), ln.stab=log(stability)) %>% 
  select(trt, plot, ln.stab, ln.syn, ln.richness, ln.Evar, ln.vas)


###
# Control structural equation models (SEM)
control <- semdata %>% 
  filter(trt=='Control')
# (1) model 1C: the 'full' model:
model_1C <- '
ln.stab ~ ln.syn + ln.richness + ln.Evar + ln.vas 
ln.syn ~ ln.richness + ln.Evar + ln.vas
ln.Evar ~ ln.vas 
ln.richness ~ ln.vas
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar
ln.richness ~~ ln.richness
ln.vas ~~ ln.vas
ln.richness ~~ ln.Evar'
set.seed(123)
fit <- lavaan(model_1C, data=control, auto.var=TRUE, se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(2) model 2C: Delete the vascular abundance (ln.vas) from model 1C
model_2C <- '
ln.stab ~ ln.syn + ln.richness + ln.Evar 
ln.syn ~ ln.richness + ln.Evar
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar
ln.richness ~~ ln.richness
ln.richness ~~ ln.Evar'
set.seed(123)
fit <-  lavaan(model_2C, data=control, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

# (3) model 3C: Delete species evenness (ln.Evar) from model 2C
model_3C <- '
ln.stab ~ ln.syn + ln.richness
ln.syn ~ ln.richness
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.richness ~~ ln.richness'
set.seed(123)
fit <- lavaan(model_3C, data=control, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(4) model 4C: Delete species richness (ln.richness) from model 2C
model_4C <- '
ln.stab ~ ln.syn + ln.Evar 
ln.syn ~ ln.Evar
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar
ln.richness ~~ ln.Evar'
set.seed(123)
fit = lavaan(model_4C, data=control, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(5) model 5C: Delete the direct path between community stability and species richness from model 3C
model_5C <- '
ln.stab ~ ln.syn
ln.syn ~ ln.richness
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.richness ~~ ln.richness'
set.seed(123)
fit <-  lavaan(model_5C, data=control, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)



###
# Warming structural equation models (SEM) 
warming <- semdata %>% 
  filter(trt=='Warming')
#(1) model 1W: the 'full' model
model_1W <- '
ln.stab ~ ln.syn + ln.richness + ln.Evar + ln.vas 
ln.syn ~ ln.richness + ln.Evar + ln.vas
ln.Evar ~ ln.vas 
ln.richness ~ ln.vas
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar
ln.richness ~~ ln.richness
ln.vas ~~ ln.vas
ln.richness ~~ ln.Evar'
set.seed(123)
fit <-  lavaan(model_1W, data=warming, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(2) model 2W: Delete the vascular abundance (ln.vas) from model 1C:
model_2W <- '
ln.stab ~ ln.syn + ln.richness + ln.Evar
ln.syn ~ ln.richness + ln.Evar
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar
ln.richness ~~ ln.richness
ln.richness ~~ ln.Evar'
set.seed(123)
fit <-  lavaan(model_2W, data=warming, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(3) model 3W: Delete species richness (ln.richness) from model 2W
model_3W <- '
ln.stab ~ ln.syn + ln.Evar
ln.syn ~ ln.Evar
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar'
set.seed(123)
fit <-  lavaan(model_3W, data=warming, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(4) model 4W: Delete species evenness (ln.Evar) from model 2W
model_4W <- '
ln.stab ~ ln.syn + ln.richness
ln.syn ~ ln.richness
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.richness ~~ ln.richness'
set.seed(123)
fit <-  lavaan(model_4W, data=warming, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

#(5) model 5W: Delete the direct path between community stability and species evenness from model 3W
model_5W <- '
ln.stab ~ ln.syn
ln.syn ~ ln.Evar
ln.stab ~~ ln.stab
ln.syn ~~ ln.syn
ln.Evar ~~ ln.Evar'
set.seed(123)
fit <-  lavaan(model_5W, data=warming, auto.var=TRUE,se='bootstrap')
summary(fit, fit.measures=TRUE, standardized = TRUE, rsquare=TRUE)

```


Fig S2 (Appendix)
```{r}
###
# Fig S2. Statistics (calculate the calibration factor for bryophytes versus lichen cover)
# model for lichen
lichen <- subset(mydata,group=='lichen')
myfit_lichen <- lme(log10(biomass) ~log10(hits), random=~1|spp,data = lichen)
r.squaredGLMM(myfit_lichen)
lichen.pred <- ggpredict(myfit_lichen)
lichen.pred <-  lichen.pred$hits
lichen.pred$category <-  'lichen'
lichen.pred <- mutate(lichen.pred, "y_t"=10^predicted,
                           "conf.high_t"=10^conf.high, "conf.low_t"=10^conf.low)
pred_lichen <- lm((y_t) ~(x), data = lichen.pred)
summary(pred_lichen)

# model for bryophytes
moss <-  subset(mydata,group=='moss')
myfit_moss <- lm(log10(biomass) ~log10(hits),data = moss)
summary(myfit_moss)
moss.pred <- ggpredict(myfit_moss)
moss.pred <-  moss.pred$hits
moss.pred$category <-  'moss'
moss.pred <- mutate(moss.pred, "y_t"=10^predicted,
                           "conf.high_t"=10^conf.high, "conf.low_t"=10^conf.low)
pred_moss <- lm((y_t) ~(x), data = moss.pred)
summary(pred_moss)

#calibration factor:
(1.1040/0.243090) 


### 
# Fig S2. Graphs
mydata <- rbind(lichen.pred,moss.pred)
global_theme <- theme_bw()+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text=element_text(size=20,colour='black'))+
  theme(axis.title=element_text(size=20))+
  theme(axis.ticks.length = unit(5,'pt'))+
  theme(legend.text = element_text(size = 20))+
  theme(legend.title = element_text(size = 20))+
  theme(legend.position = 'none')+
  theme(plot.margin=unit(c(1,1.5,1,1),"cm"))+
  theme(text=element_text(family="Times New Roman")) 

ggplot(mydata, aes(x=x, y=y_t)) + 
  geom_line(data = subset(mydata,category=='lichen'), size = 1, color = "blue") +
  geom_ribbon(data = subset(mydata,category=='lichen'), aes(x = x, ymin = conf.low_t, ymax = conf.high_t), fill = "lightblue", alpha = 0.5) +
  geom_line(data = subset(mydata,category=='moss'), size = 1, color = "green") +
  geom_ribbon(data = subset(mydata,category=='moss'), aes(x = x, ymin = conf.low_t, ymax = conf.high_t), fill = "lightblue", alpha = 0.5) +
   labs(x=expression(paste('Hits (','m'^"2",')')), y=expression(paste('Biomass (','g/m'^"2",')')), title= "") +
  scale_color_manual(values=c('green','blue')) +
  global_theme


```

  
Fig S3 (Appendix)
```{r}
###
vas_nonvas <- vas_nonvas %>% 
  mutate(plot=factor(plot),year=factor(year))

########### Fig S3. Statistics
# Fig S3(a): Relationships between cryptogam species richness and vascular abundance
#(1) Control
myfit <- lme(log(richness) ~ log(vas_abundance), random=~1|plot, data=subset(vas_nonvas,trt=='Control'))
summary(myfit)
#(2) Warming
myfit <- lme(log(richness) ~ log(vas_abundance), random=~1|plot, data=subset(vas_nonvas,trt=='Warming'))
summary(myfit)

# Fig S3(b): Relationships between cryptogam species evenness and vascular abundance
#(1) Control
myfit <- lme(log(Evar) ~ log(vas_abundance), random=~1|plot, data=subset(vas_nonvas,trt=='Control'))
summary(myfit)
#(2) Warming
myfit <- lme(log(Evar) ~ log(vas_abundance), random=~1|plot, data=subset(vas_nonvas,trt=='Warming'))
summary(myfit)


### Fig S3. Graphs
global_theme <- theme_bw()+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text=element_text(size=30,colour='black'))+
  theme(axis.title=element_text(size=30))+
  theme(axis.ticks.length = unit(10,'pt'))+
  theme(legend.text = element_text(size = 30))+
  theme(legend.title = element_text(size = 30))+
  theme(legend.position = 'right')+
  theme(legend.title = element_blank())+
  theme(plot.margin=unit(c(1,1.5,1,1),"cm"))+
  theme(text=element_text(family="Times New Roman")) 

#(1) Relationships between vascular abundance and cryptogam species richness
vas_richness <-  ggplot(vas_nonvas, aes(x=log(vas_biomass), y=log(richness))) + 
  geom_point(aes(color=trt,fill=trt,shape=year),size=6) +
    labs(x="\nln(Vascular plant abundance)\n", y="\nln(Species richness)\n") +
  scale_fill_manual(values=c('black','red'))+
  scale_color_manual(values=c('black','red'))+
  scale_shape_manual(values=c(22,23,24,25))+
  global_theme +
  scale_x_continuous(limits = c(3.4, 6.2),breaks = c(3.5, 4,4.5,5,5.5,6), expand = c(0,0))+
  scale_y_continuous(limits = c(0.8, 3.2),breaks = c(1,1.5,2,2.5,3), expand = c(0,0))


#(2) Relationships between vascular abundance and  cryptogam species evenness
vas_evenness <-  ggplot(vas_nonvas, aes(x=log(vas_biomass), y=log(Evar))) + 
  geom_point(aes(color=trt,fill=trt,shape=year),size=6) +
  stat_smooth(method='lm',color='red',se=F,size=2,data=subset(vas_nonvas,trt=='W'))+
  labs(x="\nln(Vascular plant abundance)\n", y="\nln(Species evenness)\n") +
  scale_fill_manual(values=c('black','red'))+
  scale_color_manual(values=c('black','red'))+
  scale_shape_manual(values=c(22,23,24,25))+
  global_theme +
  scale_x_continuous(limits = c(3.6, 6.2),breaks = c( 4,4.5,5,5.5,6), expand = c(0,0)) +
  scale_y_continuous(limits = c(-1.8, 0.2),breaks = c(-1.5,-1,-0.5,0), expand = c(0,0))

fig.S3 <- (vas_richness|vas_evenness)

```



fig. S4 (Appendix).
```{r}
##
#Fig S4(a): Species richness changes over time
# 2005
myfit <- lm(log(richness)~trt, data=subset(vas_nonvas,year==2005))
summary(myfit)
# No difference between control and warming plots.
# 2008
myfit <- lm(log(richness)~trt, data=subset(vas_nonvas,year==2008))
summary(myfit)
# No difference between control and warming plots.
# 2011
myfit <- lm(log(richness)~trt, data=subset(vas_nonvas,year==2011))
summary(myfit)
# No difference between control and warming plots.
# 2017
myfit <- lm(log(richness)~trt, data=subset(vas_nonvas,year==2017))
summary(myfit)
# Warming lower than control.


#Fig S4(b): Species richness changes over time
# 2005
myfit <- lm(log(Evar)~trt, data=subset(vas_nonvas,year==2005))
summary(myfit)
# No difference between control and warming plots.
# 2008
myfit <- lm(log(Evar)~trt, data=subset(vas_nonvas,year==2008))
summary(myfit)
# No difference between control and warming plots.
# 2011
myfit <- lm(log(Evar)~trt, data=subset(vas_nonvas,year==2011))
summary(myfit)
# Warming higher than control.
# 2017
myfit <- lm(log(Evar)~trt, data=subset(vas_nonvas,year==2017))
summary(myfit)
# Warming higher than control.


###
# Fig S4. Graphs
fig.S4a <-  ggplot(vas_nonvas, aes(x=year, y=log(richness),fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
  labs(x="", y="\nln(Species richness)\n") +
  scale_fill_manual(values=c('gray','red'))+
  global_theme +
  scale_y_continuous(limits = c(0.8, 3.2),breaks = c(1,1.5,2,2.5,3), expand = c(0,0))

fig.S4b <-  ggplot(vas_nonvas, aes(x=year, y=log(Evar),fill=trt)) + 
  geom_boxplot(lwd=2,outlier.size=3)+
  labs(x="", y="\nln(Species evenness)\n") +
  scale_fill_manual(values=c('gray','red'))+
  global_theme +
  scale_y_continuous(limits = c(-1.8, 0.2),breaks = c(-1.5,-1,-0.5,0), expand = c(0,0))
# Combine panels together:
fig.S4 <- (a|b)


```


