#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# manuscript: Interactive persistent effects of past land-cover and its trajectory on tropical freshwater biodiversity
# authors: Edineusa Pereira Santos, Sílvio F. B. Ferraz, Tadeu Siqueira
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#loading packages

library(easypackages)
libraries("dplyr","vegan","lme4","car","MuMIn","lmerTest","effects","ggplot2", "patchwork")

######################################landscape data
corumb_rb<-read.table("corumbatai_rb.txt", header=T)

#id class
corumb_rb$catchment<-as.factor(corumb_rb$catchment)
corumb_rb$stream<-as.factor(corumb_rb$stream)

#new variables

# ***trajectory = forest cover final minus initial
corumb_rb$trajectory<-corumb_rb$forest2011-corumb_rb$forest1962

# forest trend across years
new_corumb<-corumb_rb[corumb_rb$trajectory!=0,] 
new_corumb$traj_ctg<-rep("Gain")
new_corumb$traj_ctg[new_corumb$trajectory<0]<-'Loss'

trend<-new_corumb%>%
  group_by(traj_ctg)%>%
  select(-catchment, -stream, -trajectory, -traj_ctg)%>%
  summarise_all(funs(mean,sd))

All=c(apply(new_corumb[,3:7], MARGIN = 2, FUN="mean"),apply(new_corumb[,3:7], MARGIN = 2, FUN="sd"))
trend=rbind(trend, c("All",All))
mean1= as.vector(unlist(c(trend[1,2:6],trend[2,2:6],trend[3,2:6])))
sd1=as.numeric(as.vector(unlist(c(trend[1,7:11],trend[2,7:11],trend[3,7:11]))))
trend=as.data.frame(cbind(Trajectory=rep(c("Gain", "Loss","All"), each=5), Year=rep(c("1962", "1972","1978","2003","2011"),times=3),mean1, sd1)) 
trend$mean1<-as.numeric(as.character(trend$mean1))
trend$sd1<-as.numeric(as.character(trend$sd1))

# repeated items across plots
same<-list (
  theme_bw(),
  theme(legend.position = "none"),
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
  ))


#figure_3
  ggplot(trend, aes(x=Year, y=mean1, group=Trajectory, colour=Trajectory))+
  geom_errorbar(aes(ymin=mean1-sd1, ymax=mean1+sd1), width=.1)+
  geom_line(aes(color=Trajectory))+
  geom_point()+
  scale_colour_manual(labels = c("All","Forest gain","Forest loss"),
                      values=c("#000000","#009E73","#D55E00"))+
  same+
  theme(legend.key.size = unit(0.2, 'lines'),
        legend.position=c(0.8,0.85),
        legend.background = element_rect(fill = "transparent"))+
  labs(y="Forest cover (%)", x="Year", color=NULL)



## ***forest cover mean and standard deviation
corumb_rb$mean.fc<-apply(select(corumb_rb, forest1962,forest2011,forest1972,forest1978,forest2003),MARGIN = 1, FUN=mean)
corumb_rb$sd.fc<-apply(select(corumb_rb, forest1962,forest2011,forest1972,forest1978,forest2003),MARGIN = 1, FUN=sd)

################################## community data
comm<-read.table("aquatic_insects.txt", header=T)

# abundance
comm$abund<-rowSums(comm[,-c(1,2)])

# rarefied richness of samples with abundance>100
raref.rich<-rarefy(comm%>%
                     filter(abund>100)%>% 
                     select(-stream, -catchment,-abund), sample=100)

# diversity indices based on TSallis entropy
indices<-eventstar(comm%>%
                     filter(abund>100)%>% 
                     select(-stream, -catchment,-abund))

# joining landscape and community data
data.all<-cbind(corumb_rb, abund=comm$abund)
data.all<-filter(data.all, abund>100)
data.all$raref.rich<-raref.rich
data.all$tsal.div<-indices$Hstar


#----------------------------PART I------------------------
#----------------------------------------------------------
#--------- models including forest cover of each year
# values below 10% of the total variation (sd.fc) were removed because stability over 
# the all years blurs the association of community response with a specific year.
#----------------------------------------------------------
# --------response variable==RAREFIED RICHNESS

mfyear<-lmer(raref.rich~forest2011+forest2003+forest1978+forest1972+forest1962 +(1|catchment), data=data.all[data.all$sd>=4.6,], REML=T)
summary(mfyear)

## residuals
plot(residuals(mfyear))
qqnorm(resid(mfyear)) 
qqline(resid(mfyear))

# % varation explained by model
r.squaredGLMM(mfyear)

#betas
stand_beta<-sjstats::std_beta(mfyear)
betasf<-as.data.frame(stand_beta[-1,])

# Anova (sum of square type II)--wt interaction
Anova(mfyear, type= "II")

#---------------response variable==ABUNDANCE
mfyeab<-glmer(abund~forest2011+forest2003+forest1978+forest1972+forest1962 +(1|catchment), data=data.all[data.all$sd>=4.6,], family=poisson(link="log"))
summary(mfyeab) 

## residuals
plot(residuals(mfyeab))
qqnorm(resid(mfyeab)) #variancia aumenta com a media
qqline(resid(mfyeab))

# % varation explained by model
#Warning message:
#'r.squaredGLMM' now calculates a revised statistic. We used previous version.
r.squaredGLMM(mfyeab)

#betas
std_beta_abundf<-sjstats::std_beta(mfyeab)
betas_abundf<-as.data.frame(std_beta_abundf[-1,])

# --------response variable==TSALLIS DIVERSITY
mfyediv<-lmer(tsal.div~forest2011+forest2003+forest1978+forest1972+forest1962 +(1|catchment), data=data.all[data.all$sd>=4.6,], REML=T)
summary(mfyediv)

#betas
sjstats::std_beta(mfyediv)

#------- plotting results
# richness
betasf$term<-as.factor(betasf$term)

rich<-ggplot(betasf, aes(betasf$term, betasf$std.estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=betasf$conf.low, ymax=betasf$conf.high), width=0.1) + 
  geom_hline(yintercept=0, linetype="dashed")+
  same+
  labs(y="Standardized beta coefficient", x=NULL) +
  annotate(geom="text", x=2.5, y=-0.85, label="y= rarefied richness")+
  theme(axis.title.x = element_blank(),axis.text.x = element_blank())+
  scale_y_continuous(limits = c(-1,1))

# abundance
abund<-ggplot(betas_abundf, aes(betas_abundf$term ,betas_abundf$std.estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=betas_abundf$conf.low, ymax=betas_abundf$conf.high), width=0.1) + 
  geom_hline(yintercept=0, linetype="dashed")+
  same+
  annotate(geom="text", x=2, y=-0.00085, label="y= abundance")+
   ylab("Standardized beta coefficient") +xlab("Year") +scale_x_discrete(labels=c("1962", "1972", "1978","2003","2011"))+
  scale_y_continuous(limits = c(-9.537982e-04,0.001),breaks=pretty(betas_abundf$std.estimate, n=4), labels = function(x) format(x, scientific = FALSE))

#Figure_4
rich/abund

#------------------------------ PART II--------------------------

# removing null trajectory
data.all<-
  filter(data.all,trajectory!=0) 

#--------- models including central tendency measures and 
# spread of forest cover values across years
cor(select(data.all,mean.fc,sd.fc, trajectory))


# Scaling is needed because variables have different ranges 
data.allsc<-data.all
data.allsc[,8:10]<-scale(data.allsc[,8:10])

# models
# All models presented in table [supp. material] were REML=T.
# Here, identical models show REML=F to compare models with anova 
# models with RICHNESS as response variable
m01<-lmer(raref.rich~mean.fc*trajectory*sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m01)
sjstats::std_beta(m01) ##first, replace REML= FALSE by TRUE

m02<-lmer(raref.rich~mean.fc*trajectory+sd.fc*trajectory+(1|catchment), data=data.allsc, REML=F)
summary(m02)
sjstats::std_beta(m02)

m03<-lmer(raref.rich~mean.fc+trajectory+sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m03)
sjstats::std_beta(m03)

m04<-lmer(raref.rich~mean.fc*trajectory+(1|catchment), data=data.allsc, REML=F)
summary(m04) # the lowest AIC=402,3
sjstats::std_beta(m04)

m05<-lmer(raref.rich~mean.fc*sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m05)
sjstats::std_beta(m05)

m06<-lmer(raref.rich~sd.fc*trajectory +(1|catchment), data=data.allsc, REML=F)
summary(m06)
sjstats::std_beta(m06)

m07<-lmer(raref.rich~mean.fc +(1|catchment), data=data.allsc, REML=F)
summary(m07)
sjstats::std_beta(m07)

m08<-lmer(raref.rich~trajectory +(1|catchment), data=data.allsc, REML=F)
summary(m08)
sjstats::std_beta(m08)

m09<-lmer(raref.rich~sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m09)
sjstats::std_beta(m09)

# model selection
#anova
anova(m01,m02, m03,m04, m05,m06, m07,m08, m09)

library(bbmle)
AICctab(m02, m04,  m07, m08, base = T, weights = T) 

# final model
m04f<-lmer(raref.rich~mean.fc*trajectory+(1|catchment), data=data.allsc, REML=T)
summary(m04f)

#betas 
sjstats::std_beta(m04f)

## residuals
plot(residuals(m04f))
qqnorm(resid(m04f)) #no problem
qqline(resid(m04f))

# Anova (sum of square type III)
Anova(m04f, type="III")

# % varation explained by model
r.squaredGLMM(m04f)

# models with ABUNDANCE as response variable
#
m1.1<-glmer(abund~mean.fc*trajectory*sd.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.1) #Correlation of Fixed Effects>0.70
vif(m1.1)# vif of sd.fc 
sjstats::std_beta(m1.1)

m1.2<-glmer(abund~mean.fc*trajectory+mean.fc*sd.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.2) #Correlation of Fixed Effects>0.68
vif(m1.2)
sjstats::std_beta(m1.2)

m1.3<-glmer(abund~mean.fc+trajectory+sd.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.3) #AIC= 9212.2
vif(m1.2)
sjstats::std_beta(m1.3)

m1.4<-glmer(abund~mean.fc*trajectory +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.4) #AIC= 9359.1 
vif(m1.4)
sjstats::std_beta(m1.4)

m1.5<-glmer(abund~mean.fc*sd.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.5)#AIC=10444.6
sjstats::std_beta(m1.5)

m1.6<-glmer(abund~sd.fc*trajectory+(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.6)#AIC=10443.3
sjstats::std_beta(m1.6)

m1.7<-glmer(abund~mean.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.7) #AIC=9934.4
sjstats::std_beta(m1.7)

m1.8<-glmer(abund~trajectory +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.8)#10445.3
sjstats::std_beta(m1.7)

m1.9<-glmer(abund~sd.fc +(1|catchment), data=data.allsc,  family=poisson(link="log"))
summary(m1.9)#10445.3
sjstats::std_beta(m1.9)

# choosing the model
# model 2.6 seems the best because 2.2 has stronger correlation between mean and sd.fc 
anova (m1.3, m1.4, m1.5, m1.6, m1.7, m1.8, m1.9)
anova(m1.3, m1.4, m1.9)  
AICtab(m1.3, m1.4, base = T, weights = T) 

#betas
sjstats::std_beta(m1.4)

#residuals distribution
plot(residuals(m1.4))
qqnorm(resid(m1.4))
qqline(resid(m1.4))

# % explained by model
r.squaredGLMM(m1.4)

#sum of square type III
Anova(m1.4, type="III")

# models with TSALLIS DIVERSITY as response variable
# Here, models show REML=F to compare models with anova 
# All models presented in table [supp. material] were REML=T, just replace it.
# and calculate the standized beta
m2.1<-lmer(tsal.div~mean.fc*trajectory*sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.1) 
sjstats::std_beta(m2.1) #first, replace REML= FALSE by TRUE

m2.2<-lmer(tsal.div~mean.fc*trajectory + mean.fc*sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.2)
sjstats::std_beta(m2.2)

m2.3<-lmer(tsal.div~mean.fc+ trajectory + sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.3)
sjstats::std_beta(m2.3)

m2.4<-lmer(tsal.div~mean.fc*trajectory +(1|catchment), data=data.allsc, REML=F)
summary(m2.4) 
sjstats::std_beta(m2.4)

m2.5<-lmer(tsal.div~mean.fc*sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.5) # AIC=
sjstats::std_beta(m2.5)

m2.6<-lmer(tsal.div~sd.fc*trajectory +(1|catchment), data=data.allsc, REML=F)
summary(m2.6) # AIC=
sjstats::std_beta(m2.6)

m2.7<-lmer(tsal.div~mean.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.7)
sjstats::std_beta(m2.7)

m2.8<-lmer(tsal.div~trajectory +(1|catchment), data=data.allsc, REML=F)
summary(m2.8)
sjstats::std_beta(m2.8)

m2.9<-lmer(tsal.div~sd.fc +(1|catchment), data=data.allsc, REML=F)
summary(m2.9)
sjstats::std_beta(m2.9)

# model selection
anova (m2.1, m2.2, m2.3, m2.4, m2.5, m2.6, m2.7, m2.8, m2.9)
anova(m2.2,m2.4,m2.7,m2.8)  
anova(m2.4, m2.8)
AICctab(m2.4, m2.8, base = T, weights = T) 

#final model
m2.4f<-lmer(tsal.div~mean.fc*trajectory +(1|catchment), data=data.allsc, REML=T)
summary(m2.4f)# REML=T
sjstats::std_beta(m2.4f)

## residuals
plot(residuals(m2.4f))
qqnorm(resid(m2.4f))
qqline(resid(m2.4f))

# % varation explained by model
r.squaredGLMM(m2.4f)

#betas
stand_beta_m2.4f<-sjstats::std_beta(m2.4f)

###############################################################
############ Graphics of interaction terms
# Predict plot with firt 10% of trajectory

#----richness
m04f<-lmer(raref.rich~scale(mean.fc)*scale(trajectory)+(1|catchment), data=data.all, REML=T)
pred_rich_10<-as.data.frame(effect("mean.fc:trajectory",m04f,
                                    xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-10,0,10))))
pred_rich_10$trajectory<-as.factor(pred_rich_10$trajectory)
#----abudance
m1.4<-glmer(abund~scale(mean.fc)*scale(trajectory) +(1|catchment), data=data.all,  family=poisson(link="log"))
pred_abund_10<-as.data.frame(effect("mean.fc:trajectory",m1.4,
                                  xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-10,0,10))))
pred_abund_10$trajectory<-as.factor(pred_abund_10$trajectory)

#----Tsallis
m2.4f<-lmer(tsal.div~scale(mean.fc)*scale(trajectory) +(1|catchment), data=data.all, REML=T)
pred_tsal_10<-as.data.frame(effect("mean.fc:trajectory",m2.4f,
                                  xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-10,0,10))))
pred_tsal_10$trajectory<-as.factor(pred_tsal_10$trajectory)

####plots 10
# repeated color
same2<-list (geom_line(size = 1),geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .3, colour = NA),
               scale_color_manual(values=c("#D55E00","#999999","#009E73")),
               scale_fill_manual(values=c("#D55E00","#999999","#009E73")))
r10<-
  ggplot(pred_rich_10, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  same2+
  labs(y=NULL,x=NULL)+
  scale_y_continuous(limits = c(6.5,20))+
  same

ab10<-
  ggplot(pred_abund_10, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  same2+
  labs(y=NULL, x=NULL)+ 
  scale_y_continuous(limits = c(80,1850))+  
  same

ts10<-
  ggplot(pred_tsal_10, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  same2+
  labs(y=NULL , x=NULL , colour="Trajecory", fill="Trajecory")+#,x="% of forest cover (mean)" "Tsallis diversity"
  scale_y_continuous(limits = c(0.5,5))+
  same

# Predict plot with firt 20% of trajectory
#----species richness
pred_rich_20<-as.data.frame(effect("mean.fc:trajectory",m04f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-20,0,20))))
pred_rich_20$trajectory<-as.factor(pred_rich_20$trajectory)

#----abudance
pred_abund_20<-as.data.frame(effect("mean.fc:trajectory",m1.4,
                                    xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-20,0,20))))
pred_abund_20$trajectory<-as.factor(pred_abund_20$trajectory)

#----Tsallis
pred_tsal_20<-as.data.frame(effect("mean.fc:trajectory",m2.4f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-20,0,20))))
pred_tsal_20$trajectory<-as.factor(pred_tsal_20$trajectory)

#### plots 20
r20<-
  ggplot(pred_rich_20, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_rich_20$mean.fc, y=0, xend=pred_rich_20$mean.fc, yend=6.5, colour="black")+
  same2+
  labs(y=NULL,x=NULL)+
  scale_y_continuous(limits = c(6.5,20))+
  same

ab20<-
  ggplot(pred_abund_20, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_abund_20$mean.fc, y=-10, xend=pred_abund_20$mean.fc, yend=80, colour="black")+
  same2+
  labs(y=NULL, x=NULL)+ #, x="% of forest cover (mean)" "% of forest cover (mean)"
  scale_y_continuous(limits = c(80,1850))+  #4,1850
  same

ts20<-
  ggplot(pred_tsal_20, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_tsal_20$mean.fc, y=0, xend=pred_tsal_20$mean.fc, yend=0.5, colour="black")+
  same2+
  labs(y=NULL , x=NULL , colour=NULL, fill=NULL)+#,x="% of forest cover (mean)" "Tsallis diversity"
  scale_y_continuous(limits = c(0.5,5))+
  same


# Predict plot with firt 30% of trajectory
#----species richness
pred_rich_30<-as.data.frame(effect("mean.fc:trajectory",m04f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-30,0,30))))
pred_rich_30<-pred_rich_30[pred_rich_30$mean.fc<=70,]
pred_rich_30$trajectory<-as.factor(pred_rich_30$trajectory)

#-----abundance
pred_abund_30<-as.data.frame(effect("mean.fc:trajectory",m1.4,
                                  xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-30,0,30))))
pred_abund_30<-pred_abund_30[pred_abund_30$mean.fc<=70,]
pred_abund_30$trajectory<-as.factor(pred_abund_30$trajectory)

#----- Tsallis diversity
pred_tsal_30<-as.data.frame(effect("mean.fc:trajectory",m2.4f,
                                  xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-30,0,30))))
pred_tsal_30<-pred_tsal_30[pred_tsal_30$mean.fc<=70,]
pred_tsal_30$trajectory<-as.factor(pred_tsal_30$trajectory)

#### plots 30
r30<-
  ggplot(pred_rich_30, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_rich_30$mean.fc, y=0, xend=pred_rich_30$mean.fc, yend=6.5, colour="black")+
  same2+
  labs(y="Rarefied richness" ,x=NULL, colour=NULL, fill=NULL)+ #"Rarefied richness" "% of forest cover (mean)"
  scale_y_continuous(limits = c(6.5,20))+
  same
  

ab30<-
  ggplot(pred_abund_30, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_abund_30$mean.fc, y=0, xend=pred_abund_30$mean.fc, yend=80, colour="black")+
  same2+
  labs(y="Abundance",x=NULL)+#
  scale_y_continuous(limits = c(80,1850))+
  same

ts30<-
  ggplot(pred_tsal_30, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_tsal_30$mean.fc, y=0, xend=pred_tsal_30$mean.fc, yend=0.5, colour="black")+
  same2+
  labs(y="Tsallis diversity",x=NULL, colour=NULL, fill=NULL)+ #
  scale_y_continuous(limits = c(0.5,5))+
  same

# Predict plot with firt 40% of trajectory
#----species richness
pred_rich_40<-as.data.frame(effect("mean.fc:trajectory",m04f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-40,0,40))))
pred_rich_40<-pred_rich_40[pred_rich_40$mean.fc<=60,]
pred_rich_40$trajectory<-as.factor(pred_rich_40$trajectory)

#-----abundance
pred_abund_40<-as.data.frame(effect("mean.fc:trajectory",m1.4,
                                    xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-40,0,40))))
pred_abund_40<-pred_abund_40[pred_abund_40$mean.fc<=60,]
pred_abund_40$trajectory<-as.factor(pred_abund_40$trajectory)

#----- Tsallis diversity
pred_tsal_40<-as.data.frame(effect("mean.fc:trajectory",m2.4f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-40,0,40))))
pred_tsal_40<-pred_tsal_40[pred_tsal_40$mean.fc<=60,]
pred_tsal_40$trajectory<-as.factor(pred_tsal_40$trajectory)

##### Plots 40
r40<-
  ggplot(pred_rich_40, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_rich_40$mean.fc, y=0, xend=pred_rich_40$mean.fc, yend=6.5, colour="black")+
  same2+
  labs(y=NULL,x=NULL)+
  scale_y_continuous(limits = c(6.5,20))+
  same

ab40<-
  ggplot(pred_abund_40, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_abund_40$mean.fc, y=-10, xend=pred_abund_40$mean.fc, yend=80, colour="black")+
  same2+
  labs(y=NULL, x=NULL)+ #, x="% of forest cover (mean)" "% of forest cover (mean)"
  scale_y_continuous(limits = c(80,1850))+  #4,1850
  same

ts40<-
  ggplot(pred_tsal_40, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_tsal_40$mean.fc, y=0, xend=pred_tsal_40$mean.fc, yend=0.5, colour="black")+
  same2+
  labs(y=NULL , x=NULL , colour=NULL, fill=NULL)+#,x="% of forest cover (mean)" "Tsallis diversity"
  scale_y_continuous(limits = c(0.5,5))+
  same


# Predict plot with firt 50% of trajectory
#----species richness
pred_rich_50<-as.data.frame(effect("mean.fc:trajectory",m04f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-50,0,50))))
pred_rich_50<-pred_rich_50[pred_rich_50$mean.fc<=50,]
pred_rich_50$trajectory<-as.factor(pred_rich_50$trajectory)

#-----abundance
pred_abund_50<-as.data.frame(effect("mean.fc:trajectory",m1.4,
                                    xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-50,0,50))))
pred_abund_50<-pred_abund_50[pred_abund_50$mean.fc<=50,]
pred_abund_50$trajectory<-as.factor(pred_abund_50$trajectory)

#----- Tsallis diversity
pred_tsal_50<-as.data.frame(effect("mean.fc:trajectory",m2.4f,
                                   xlevels=list(mean.fc=data.all$mean.fc, trajectory=c(-50,0,50))))
pred_tsal_50<-pred_tsal_50[pred_tsal_50$mean.fc<=50,]
pred_tsal_50$trajectory<-as.factor(pred_tsal_50$trajectory)


# Plots 50
##### Plots 50
r50<-
  ggplot(pred_rich_50, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_rich_50$mean.fc, y=0, xend=pred_rich_50$mean.fc, yend=6.5, colour="black")+
  same2+
  labs(y=NULL,x="% of forest cover (mean)")+
  scale_y_continuous(limits = c(6.5,20))+
  same

ab50<-
  ggplot(pred_abund_50, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_abund_50$mean.fc, y=-10, xend=pred_abund_50$mean.fc, yend=80, colour="black")+
  same2+
  labs(y=NULL, x="% of forest cover (mean)")+ #, x="% of forest cover (mean)" "% of forest cover (mean)"
  scale_y_continuous(limits = c(80,1850))+  #4,1850
  same


ts50<-
  ggplot(pred_tsal_50, aes(mean.fc, fit, color = trajectory, fill = trajectory)) +
  geom_segment(x =pred_tsal_50$mean.fc, y=0, xend=pred_tsal_50$mean.fc, yend=0.5, colour="black")+
  same2+
  labs(y=NULL, x="% of forest cover (mean)", colour=NULL, fill=NULL)+ #, x="% of forest cover (mean)" "% of forest cover (mean)"
  scale_y_continuous(limits = c(0.5,5))+
  same

# Figure_5
(r10|ab10|ts10)/(r20| ab20|ts20)/
  (r30|ab30|ts30)/(r40|ab40|ts40)/
  (r50|ab50|ts50)