ls()
remove(list=ls())

library(ade4)
library(packfor)
library(AEM)
library(plyr)
library(nlme)

source("rsquared.glmm.R")
SpecAb<-as.data.frame(read.csv("SpeciesSubplot.csv", sep = ",",header=TRUE))
Traits<-as.data.frame(read.csv("Traitsbewerkt2.csv", sep = ",",header=TRUE))
Traits<-Traits[which(Traits$Species%in%colnames(SpecAb)),]
SpecAb<-SpecAb[,which(colnames(SpecAb)%in%Traits$Species)]
# The first is called the marginal R2 and describes the proportion of variance explained by the fixed factor(s) alone. The second is the conditional R2, which describes the proportion of variance explained by both the fixed and random factors
anova(M5.lme)

Form1<-formula(AGC ~ Occupancy+SpeciesN+Groups+Planted+Simpson, data=All)
M1.lme <- gls(Form1,   method = "REML", data = All , na.action=na.exclude)
Form4<-formula(AGC ~ Occupancy+SpeciesN+Div_Pl+Groups+Simpson, data=All)
M4.lme <- lme(Form4, random = ~1|Plot,  method = "ML", data = All , na.action=na.exclude)
Form5<-formula(AGC ~ Occupancy+SpeciesN+Div_Pl+Groups, data=All)
M5.lme <- lme(Form5, random = ~1|Plot,  method = "ML", data = All , na.action=na.exclude)
Form6<-formula(log(AGC) ~ Occupancy+SpeciesN+Groups, data=All)
M6.lme <- lme(Form6, random = ~1|Plot,  method = "ML", data = All , na.action=na.exclude)
All$Occupancy<-scale(All$Occupancy)
All$SpeciesN<-scale(All$SpeciesN)
anova(M4.lme,M5.lme,M6.lme)
summary(M6.lme)
anova.lme(M6.lme)

0.5 * (1 - pchisq(5.159065, 1))
0.5 * ((1 - pchisq(8.447942, 1)) + (1 - pchisq(8.447942, 2)))
summary(M3.lme)

Form31a<-formula(log(AGC) ~ Occupancy+Groups+SpeciesN, data=All)
M31a.lme <- lme(Form31a, random = ~1|Plot,  method = "REML",control=ctrl, data = All , na.action=na.exclude)
summary(M31a.lme)
anova(M31a.lme)
AIC(M31a.lme)
r.squaredGLMM(M31a.lme)
# lmm<-lm(AGC~Occupancy,data=All)
# anova(lmm)
# summary(lmm)
Form31b<-formula(AGC ~ Groups+Occupancy+SpeciesN, data=All)
M31b.lme <- lme(Form31b, random = ~1+Occupancy|Plot,  method = "ML",control=ctrl, data = All , na.action=na.exclude)
summary(M31b.lme)
AIC(M31b.lme)
r.squaredGLMM(M31b.lme)
anova(M31a.lme,M31b.lme)

res_lme=residuals(M31a.lme)
qqnorm(res_lme,main="M31")
qqline(res_lme)
hist(res_lme,main="M31")
r.squaredGLMM(M31a.lme)
plot(res_lme)

Form21<-formula(AGC ~ Occupancy+SpeciesN+Groups, data=All) 
M21.lme <- lme(Form21, random = ~1+Occupancy|Plot,  method = "REML", data = All , na.action=na.exclude)
summary(M21.lme)
AIC(M21.lme)
res_lme=residuals(M21.lme)
qqnorm(res_lme,main="M21")
qqline(res_lme)
hist(res_lme,main="M21")
r.squaredGLMM(M21.lme)
plot(res_lme)

Form21b<-formula(AGC ~ Occupancy+SpeciesN+Groups, data=All)
M21b.lme <- lme(Form21b, random = ~ 1|Plot,  method = "REML", data = All , na.action=na.exclude)
anova(M21.lme,M21b.lme)

ctrl <- lmeControl(opt='optim');

SpecTraits<-as.data.frame(read.csv('Species_Traits4.csv'))
par(mfrow = c(1,1 ),mar=c(2.6,3,1.2,0.4),mgp=c(1.75,0.3,0),oma=c(0.5,0.5,0.5,0.5),cex.axis=0.7)
SpecTraits$Light<-factor(SpecTraits$Light)
SpecTraits$Lebrun1<-factor(SpecTraits$Lebrun1)

dis2<-gowdis(SpecTraits[c(1:22),c(4,7)])
clus2<-hclust(dis2)
plot(clus2)
rect.hclust(clus2, 4)
plot(clus2)
rect.hclust(clus2, 5)

SpecTraits$grp2<- cutree(clus2,4)
Tab<-table(SpecTraits$grp2,SpecTraits$Light)
chisq.test(Tab) 
kruskal.test(SpecTraits$WD~SpecTraits$grp2)
aggregate(SpecTraits$WD,by=list(SpecTraits$grp2),mean)
aggregate(SpecTraits$WD,by=list(SpecTraits$grp2),sd)


All$Cl1<-NULL
All$Cl2<-NULL
All$Cl3<-NULL
All$Cl4<-NULL

Classes<-as.data.frame(read.csv('Class_10.csv'))
All<-merge.with.order(All,Classes, by=c('Plot'), all.x = T, sort=F ,keep_order = 1)
All$Cl1<-factor(All$Cl1)
All$Cl2<-factor(All$Cl2)
All$Cl3<-factor(All$Cl3)
All$Cl4<-factor(All$Cl4)
# All$Cl5<-factor(All$Cl5)
All$Nursing<-factor(All$Nursing)
# All$Class.y<-factor(All$Class.y)
# All$Cl1.y<-factor(All$Cl1.y)
# All$Cl2.y<-factor(All$Cl2.y)
# All$Cl3.y<-factor(All$Cl3.y)
# All$Class1<-factor(All$Class1)
# All$Class2<-factor(All$Class2)


Form3<-formula(log(AGC) ~ Groups+Occupancy+SpeciesN+Nursing, data=All)
M3.lme <- lme(Form3, random = ~1+Occupancy|Plot,  method = "ML",control=ctrl, data = All , na.action=na.exclude)
summary(M3.lme)
AIC(M3.lme)
qqnorm(res_lme,main="M1")
qqline(res_lme)
r.squaredGLMM(M3.lme)
Form4<-formula(log(AGC) ~ Cl1+Cl2+Cl3+Cl4+Nursing+Occupancy+SpeciesN, data=All)
M4.lme <- lme(Form4, random = ~1+Occupancy|Plot,  method = "ML",control=ctrl, data = All , na.action=na.exclude)
summary(M4.lme)
AIC(M4.lme)
r.squaredGLMM(M4.lme)
anova(M3.lme,M4.lme)

par(mfrow = c( 1,3 ),mar=c(2.6,3,1.2,0.4),mgp=c(1.75,0.3,0),oma=c(0.5,0.5,0.5,0.5),cex.axis=1)
par(xpd=TRUE)
plot(All$Occupancy,All$AGC,tck=0.005,las=1,ylab=expression(text="AGC (Mg C ha"^"-1"*")"),xlab=expression(text="BA"['pl']),cex.lab=1.1)
text(0.02,720,expression(bold(text="a")))
rp = vector('expression',2)
mod1=cor.test(All$Occupancy,All$AGC,method="spearman")
r=mod1$estimate
my.p=mod1$p.value
rp[1] = substitute(expression(italic(r) == MYVALUE), list(MYVALUE = format(r,dig=3)))[2]
rp[2] = substitute(expression(italic(p) == MYOTHERVALUE), list(MYOTHERVALUE = format(my.p, digits = 2)))[2]
legend('topright', legend = rp, bty = 'n')
plot(All$Occupancy,All$SpeciesN,tck=0.005,las=1,xlab=expression(text="BA"['pl']),ylab="Effective Species Richness",cex.lab=1.1)
text(0.02,17.92,expression(bold(text="b")))
rp = vector('expression',2)
mod1=cor.test(All$Occupancy,All$SpeciesN,method="spearman")
r=mod1$estimate
my.p=mod1$p.value
rp[1] = substitute(expression(italic(r) == MYVALUE), list(MYVALUE = format(r,dig=3)))[2]
rp[2] = substitute(expression(italic(p) == MYOTHERVALUE), list(MYOTHERVALUE = format(my.p, digits = 2)))[2]
legend('topright', legend = rp, bty = 'n')
plot(All$SpeciesN,All$AGC,tck=0.005,las=1,xlab="Effective Species Richness",ylab=expression(text="AGC (Mg C ha"^"-1"*")"),cex.lab=1.1)
text(1.21,720,expression(bold(text="c")))
rp = vector('expression',2)
mod1=cor.test(All$SpeciesN,All$AGC,method="spearman")
r=mod1$estimate
my.p=mod1$p.value
rp[1] = substitute(expression(italic(r) == MYVALUE), list(MYVALUE = format(r,dig=2)))[2]
rp[2] = substitute(expression(italic(p) == MYOTHERVALUE), list(MYOTHERVALUE = format(my.p, digits = 2)))[2]
legend('topright', legend = rp, bty = 'n')
par(xpd=FALSE)

dev.off()