library(ggplot2)
library(lme4)
library(lmerTest)
library(MuMIn)
library(ggeffects)
library(insight)
library(robustlmm)

## setwd to the correct directory

########################################
## data processing

## read in height and diameter
FABdata<-read.csv("FABData_2018.csv")

## easy way to convert factors to numeric
fac2num<-function(x) as.numeric(as.character(x))

## calculate relative growth rates from the first three years
## not bothering to account for mortality 2013-6 (from above)
## just because I'm just interested in mean growth rates of living trees
## biovolume estimates assume total woody volume can be approximated
## as a cylinder (r = basal diameter)
FABdata$V_2013_proj<-pi*(fac2num(FABdata$D_13)/20)^2*FABdata$H_13
FABdata$V_2014_proj<-pi*(FABdata$D_14/20)^2*FABdata$H_14
FABdata$V_2015_proj<-pi*(FABdata$D_15/20)^2*FABdata$H_15
FABdata$V_2016_proj<-pi*(FABdata$D_16/20)^2*FABdata$H_16
FABdata$initRGR<-(log(FABdata$V_2016_proj)-log(FABdata$V_2013_proj))/3

## change species code to get initial growth in monoculture
## for each species
mean(FABdata$initRGR[FABdata$Species=="ACNE" & FABdata$Treatment=="Monoculture"],na.rm=T)

## remove duplicates -- trees that were replanted after dying
FABdata_dup<-FABdata$Location[which(duplicated(FABdata$Location))]
elim<-list()
for(i in unique(FABdata_dup)){
  dup_entry<-FABdata[FABdata$Location==i,]
  elim[[i]]<-dup_entry$ID[-which.max(dup_entry$Year_planted)]
}

FABdata<-FABdata[-which(FABdata$ID %in% unlist(elim)),]

###########################
## build models to predict basal diameter in 2017 and 2018 when unmeasured

## this function creates a list of models -- one for each species
## each is a linear regression predicting 'current' year basal diameter
## from 'current' year height and 'previous' year basal diameter
## d.col, h.col, and d.col.prev are all names of columns
## in the data frame df
gen.dh.models<-function(df,sp.col,d.col,h.col,d.col.prev,dbh.col=NULL){
  model.list<-list()
  df.sp.list<-split(df,df[[sp.col]])
  model.list<-lapply(df.sp.list,function(x){
    diam<-x[,d.col]
    height<-x[,h.col]
    diam.prev<-x[,d.col.prev]
    return(lm(diam~height+diam.prev))
  })
  return(model.list)
}

## create the models to make predictions
models.2017<-gen.dh.models(FABdata,"Species","D_17","H_17","D_16")
models.2018<-gen.dh.models(FABdata,"Species","D_18","H_18","D_16")

## note: these linear models clearly do not work for PIBA in 2018
## because basal diameter was measured on only five individuals
## a workaround is shown later

## create estimates of projected basal diameter in 2017
## if measured, same as the measured value
## if unmeasured, predicted via the species-species
## regression models built above
FABdata$D_2017_proj<-apply(FABdata,1,function(x){
  if(!is.na(x[['D_17']])){
    return(fac2num(x[['D_17']]))
  } else {
    sp<-as.character(x[['Species']])
    model.pred<-models.2017[[match(sp,names(models.2017))]]
    pred.d<-predict(model.pred,
                    newdata=data.frame(height=fac2num(x[['H_17']]),
                                       diam.prev=fac2num(x[['D_16']])))
    return(pred.d)
  }
})

## same but in 2018
FABdata$D_2018_proj<-apply(FABdata,1,function(x){
  if(!is.na(x[['D_18']])){
    return(fac2num(x[['D_18']]))
  } else {
    sp<-as.character(x[['Species']])
    model.pred<-models.2018[[match(sp,names(models.2018))]]
    pred.d<-predict(model.pred,
                    newdata=data.frame(height=fac2num(x[['H_18']]),
                                       diam.prev=fac2num(x[['D_16']])))
    return(pred.d)
  }
})

## estimate biovolume in 2017 and calculate RGR for each year
FABdata$V_2017_proj<-pi*(FABdata$D_2017_proj/20)^2*FABdata$H_17
FABdata$RGR13.14<-log(FABdata$V_2014_proj)-log(FABdata$V_2013_proj)
FABdata$RGR14.15<-log(FABdata$V_2015_proj)-log(FABdata$V_2014_proj)
FABdata$RGR15.16<-log(FABdata$V_2016_proj)-log(FABdata$V_2015_proj)
FABdata$RGR16.17<-log(FABdata$V_2017_proj)-log(FABdata$V_2016_proj)

## estimate biovolume in 2018
FABdata$V_2018_proj<-pi*(FABdata$D_2018_proj/20)^2*FABdata$H_18

## specifically for PIBA, because the diameter prediction models
## are so bad, we assume RGR is the same from 2017 to 2018
## as from 2016 to 2017
FABdata$V_2018_proj[which(FABdata$Species=="PIBA")]<-with(FABdata[which(FABdata$Species=="PIBA"),],
                                                   V_2017_proj*exp(1)^RGR16.17)

## calculate RGR from 16-18 and 17-18
FABdata$RGR17.18<-log(FABdata$V_2018_proj)-log(FABdata$V_2017_proj)
FABdata$RGR16.18<-(log(FABdata$V_2018_proj)-log(FABdata$V_2016_proj))/2

## in the original data, dead trees have an NA
## here, I replace NAs for 0s for trees that died (or resprouted)
## between 2017 and 2018
## trees dead in both years don't really matter
FABdata$V_2017_proj[which(is.na(FABdata$H_17) & !is.na(FABdata$H_18))]<-0
FABdata$V_2018_proj[which(!is.na(FABdata$H_17) & is.na(FABdata$H_18))]<-0

## absolute growth rate from 2017 to 2018
FABdata$V_growth_proj<-FABdata$V_2018_proj-FABdata$V_2017_proj

###########################
## neighborhood level analyses

let2num<-function(let) return(match(let,LETTERS[1:26]))

## this function identifies the eight neighbors of a given individual
neighbor.id<-function(plot,row,column){
  plots<-rep(plot,9)
  rows<-c(rep(row-1,3),rep(row,3),rep(row+1,3))
  if(column!="A"){
    cols<-LETTERS[c(rep(let2num(column)+(-1:1),3))]
  } else {
    cols<-rep(c(NA,"A","B"),3)
  }
  
  plots[which(cols=="I" & plot%%7==0)]<-NA
  plots[which(is.na(cols) & (plot-1)%%7==0)]<-NA
  
  plots[which(cols=="I" & plot%%7!=0)]<-plot+1
  plots[which(is.na(cols) & (plot-1)%%7!=0)]<-plot-1
  cols[which(cols=="I" & plot%%7!=0)]<-"A"
  cols[which(is.na(cols) & (plot-1)%%7!=0)]<-"H"
  
  plots[which(rows==0 & plot %in% c(1:7,50:56,99:105))]<-NA
  plots[which(rows==9 & plot %in% c(43:49,92:98,141:147))]<-NA
  
  plots[which(rows==9 & !(plot %in% c(43:49,92:98,141:147)))]<-plots[which(rows==9 & !(plot %in% c(43:49,92:98,141:147)))]+7
  plots[which(rows==0 & !(plot %in% c(1:7,50:56,99:105)))]<-plots[which(rows==0 & !(plot %in% c(1:7,50:56,99:105)))]-7
  rows[which(rows==9 & !(plot %in% c(43:49,92:98,141:147)))]<-1
  rows[which(rows==0 & !(plot %in% c(1:7,50:56,99:105)))]<-8
  
  return(data.frame(plots=plots,rows=rows,cols=cols))
}

## this function runs through all individuals and
## applies the above function to grab the biomass of neighbors
neighbor.finder<-function(dat,outcome.var){
  NW.out<-NULL
  W.out<-NULL
  SW.out<-NULL
  N.out<-NULL
  S.out<-NULL
  NE.out<-NULL
  E.out<-NULL
  SE.out<-NULL
  
  for(i in 1:nrow(dat)){
    plot<-dat$Plot[i]
    row<-dat$Row[i]
    col<-dat$Col[i]
    neighbors<-neighbor.id(plot=plot,row=row,column=col)
    neighbors$cols<-factor(neighbors$cols,levels=LETTERS[1:8])
    
    out<-sapply(1:9,function(x) which(dat$Plot==neighbors$plot[x] & dat$Row==neighbors$rows[x] & dat$Col==neighbors$cols[x]))
    out[which(out=="integer(0)")]<-NA
    out<-unlist(out)
    
    NW.out[i]<-dat[out[1],outcome.var]
    N.out[i]<-dat[out[2],outcome.var]
    NE.out[i]<-dat[out[3],outcome.var]
    W.out[i]<-dat[out[4],outcome.var]
    E.out[i]<-dat[out[6],outcome.var]
    SW.out[i]<-dat[out[7],outcome.var]
    S.out[i]<-dat[out[8],outcome.var]
    SE.out[i]<-dat[out[9],outcome.var]
  }
  
  return(list(NW.out=NW.out,N.out=N.out,NE.out=NE.out,
              W.out=W.out,E.out=E.out,
              SW.out=SW.out,S.out=S.out,SE.out=SE.out))
}

## apply the above function and compile neighbor biomass
neighbor.bio<-neighbor.finder(FABdata,"V_2018_proj")

neighbor.df <- data.frame(matrix(unlist(neighbor.bio), nrow=9088, byrow=F),stringsAsFactors=FALSE)
colnames(neighbor.df)<-names(neighbor.bio)

## add neighbor biomass in each direction and the mean
## to each row of the main data frame
FABdata$NW.bio<-neighbor.bio$NW.out
FABdata$N.bio<-neighbor.bio$N.out
FABdata$NE.bio<-neighbor.bio$NE.out
FABdata$W.bio<-neighbor.bio$W.out
FABdata$E.bio<-neighbor.bio$E.out
FABdata$SW.bio<-neighbor.bio$SW.out
FABdata$S.bio<-neighbor.bio$S.out
FABdata$SE.bio<-neighbor.bio$SE.out

FABdata$all.neighbors<-apply(neighbor.df,1,function(x) sum(x,na.rm=T)/8)

## aggregate to the plot scale
FAB.sp.neighbor.ag<-aggregate(list(FABdata$V_growth_proj,FABdata$V_2016_proj,FABdata$V_2018_proj,
                                   FABdata$all.neighbors,FABdata$RGR16.18),
                              by=list(FABdata$Species,FABdata$Plot,FABdata$Treatment,FABdata$BEPA_Treatment),
                            FUN=function(x) mean(x,na.rm=T))
colnames(FAB.sp.neighbor.ag)<-c("Species","Plot","Treatment","BEPA_Treatment",
                                "Biomass.change","V2016","V2018",
                                "Neighbor.biomass","RGR16.18ind")
FAB.sp.neighbor.ag$RGR16.18plot<-(log(FAB.sp.neighbor.ag$V2018)-log(FAB.sp.neighbor.ag$V2016))/2

## very kludgy fix to get the color scheme right for later plotting
## since ACNE has no shaded bicultures
FAB.ag.ACNE<-FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="ACNE",]
FAB.ag.ACNE<-rbind(FAB.ag.ACNE,NA)
FAB.ag.ACNE$Treatment[24]<-"Shaded biculture"

## create a copy of the data that doesn't include
## individuals on the experiment's edge
FABdata$exp.edge<-NA
FABdata$exp.edge[FABdata$Plot %in% c(1:7,50:56,99:105) & FABdata$Row==1]<-"X"
FABdata$exp.edge[FABdata$Plot %in% c(43:49,92:98,141:147) & FABdata$Row==8]<-"X"
FABdata$exp.edge[FABdata$Plot %in% (c(0:20)*7+1) & FABdata$Col=="A"]<-"X"
FABdata$exp.edge[FABdata$Plot %in% (c(1:21)*7) & FABdata$Col=="H"]<-"X"
FABdata.noedge<-FABdata[-which(FABdata$exp.edge=="X"),]

## calculate mean absolute growth in different treatments
## not including edges
FABdata_TIAM_mono<-FABdata[with(FABdata,Species=="TIAM" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Monoculture"),]
FABdata_TIAM_shade<-FABdata[with(FABdata,Species=="TIAM" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Shaded biculture"),]
FABdata_QURU_mono<-FABdata[with(FABdata,Species=="QURU" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Monoculture"),]
FABdata_QURU_shade<-FABdata[with(FABdata,Species=="QURU" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Shaded biculture"),]
FABdata_BEPA_mono<-FABdata[with(FABdata,Species=="BEPA" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Monoculture"),]
FABdata_BEPA_shade<-FABdata[with(FABdata,Species=="BEPA" & !(Row %in% c(1,8)) & !(Col %in% c("A","H")) & Treatment=="Shaded biculture"),]
mean(FABdata_TIAM_mono$V_growth_proj,na.rm=T)

## testing for etiolation
## selecting living plants in 2016
## looking at residuals of diameter~height relationship
## and checking relationship with neighbor size
## replace BEPA with other species to test others
BEPA_2016<-FABdata[FABdata$Species=="BEPA" & !is.na(FABdata$H_16) & !is.na(FABdata$D_16) & !is.na(FABdata$all.neighbors),]
BEPA_DH_2016<-lmer(D_16~H_16+(1|Plot),data=BEPA_2016)
BEPA_DH_2016_resid<-resid(BEPA_DH_2016)
summary(lm(BEPA_DH_2016_resid~BEPA_2016$all.neighbors))

DH_ratio_2016<-lmer(log(D_16/(10*H_16))~all.neighbors+(1|Plot),data=BEPA_tmp)

###########################################
## plotting: absolute growth rate

## for each species, we get the best-fit line from a mixed-effect model
## AGR (2017 to 2018)~neighbor size with Plot random effect
ACNEbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="ACNE",]),
                            terms="all.neighbors")

ACNEbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="ACNE" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

## and then we plot, using those predictions
ACNEbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="ACNE",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=ACNEbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=ACNEbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=ACNEbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=ACNEbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.ag.ACNE,
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Acer negundo")+guides(color=F)

ACRUbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="ACRU",]),
                            terms="all.neighbors")

ACRUbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="ACRU" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

ACRUbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="ACRU",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=ACRUbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=ACRUbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=ACRUbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=ACRUbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="ACRU",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"),
        legend.position = c(0.7,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Acer rubrum")+
  guides(color=guide_legend(title="Treatment"))

BEPAbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="BEPA",]),
                            terms="all.neighbors")

BEPAbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="BEPA" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

BEPAbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="BEPA",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=BEPAbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=BEPAbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=BEPAbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=BEPAbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="BEPA",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=BEPA_Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Betula papyrifera")+guides(color=F)

QUALbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="QUAL",]),
                            terms="all.neighbors")

QUALbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="QUAL" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

QUALbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="QUAL",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=QUALbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=QUALbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=QUALbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=QUALbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUAL",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Quercus alba")+guides(color=F)

QUELbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="QUEL",]),
                            terms="all.neighbors")

QUELbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="QUEL" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

QUELbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="QUEL",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=QUELbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=QUELbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=QUELbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=QUELbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUEL",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Quercus ellipsoidalis")+guides(color=F)

QUMAbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="QUMA",]),
                            terms="all.neighbors")

QUMAbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="QUMA" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

QUMAbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="QUMA",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=QUMAbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=QUMAbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=QUMAbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=QUMAbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUMA",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Quercus macrocarpa")+guides(color=F)

QURUbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="QURU",]),
                            terms="all.neighbors")

QURUbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="QURU" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

QURUbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="QURU",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=QURUbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=QURUbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=QURUbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=QURUbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QURU",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Quercus rubra")+guides(color=F)

TIAMbiomass_rlmer<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                 data=FABdata[FABdata$Species=="TIAM",]),
                            terms="all.neighbors")

TIAMbiomass_rlmer_mort<-ggpredict(rlmer(V_growth_proj~all.neighbors+(1|Plot),
                                      data=FABdata[FABdata$Species=="TIAM" & FABdata$V_2018_proj>0,]),
                                 terms="all.neighbors")

TIAMbiomass_rlmer_plot<-ggplot(FABdata[FABdata$Species=="TIAM",],aes(x=all.neighbors,y=V_growth_proj))+
  geom_point(colour="gray")+
  geom_ribbon(data=TIAMbiomass_rlmer,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_ribbon(data=TIAMbiomass_rlmer_mort,aes(x=x,ymin=conf.low,y=predicted,ymax=conf.high),alpha=0.15)+
  geom_line(data=TIAMbiomass_rlmer,aes(x=x,y=predicted),color="blue",size=1)+
  geom_line(data=TIAMbiomass_rlmer_mort,aes(x=x,y=predicted),color="red",size=1)+
  theme_bw()+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Volume growth (",cm^3,")")))+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="TIAM",],
             aes(x=Neighbor.biomass,y=Biomass.change,color=Treatment),size=4)+
  theme(text=element_text(size=25),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  coord_cartesian(xlim=c(0,9000))+
  ggtitle("Tilia americana")+guides(color=F)

pdf("../LTER_ASM_2018/biomass.pdf",width=12,height=20,onefile = F)
egg::ggarrange(plots = list(ACNEbiomass_rlmer_plot,ACRUbiomass_rlmer_plot,
                            BEPAbiomass_rlmer_plot,QUALbiomass_rlmer_plot,
                            QUELbiomass_rlmer_plot,QUMAbiomass_rlmer_plot,
                            QURUbiomass_rlmer_plot,TIAMbiomass_rlmer_plot),
               nrow=4,ncol=2)
dev.off()

###########################################
## plotting: relative growth rate

## how many points are outside the plot bounds?
sum(FABdata$RGR16.18< -2 | FABdata$RGR16.18>2.5 | FABdata$all.neighbors>9000,na.rm=T)
## and how many points are there total?
sum(!is.na(FABdata$RGR16.18) & !is.na(FABdata$all.neighbors))

## again, we make predictions from the mixed-effects model
## RGR (2016 to 2018)~neighbor size with Plot random effect
ACNE_seq<-seq(from=floor(min(FABdata$all.neighbors[FABdata$Species=="ACNE"],na.rm=T)+1),
              to=ceiling(max(FABdata$all.neighbors[FABdata$Species=="ACNE"],na.rm=T)),
              by=1)

ACNE_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~log(all.neighbors)+(1|Plot),
                           data=FABdata[FABdata$Species=="ACNE" & FABdata$all.neighbors!=0,]),
                           terms="all.neighbors [ACNE_seq]")

## and then we plot
ACNE_biomassRGR_lmer_plot<-ggplot(ACNE_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="ACNE",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.ag.ACNE,
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                  color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  theme_bw()+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(face="italic"))+
  ggtitle("Acer negundo")+guides(color=F) ## +
  # annotation_custom(ggplotGrob(ACNE_biomassRGR_lmer_inset),
  #                   xmin = 6000, xmax = 9000, 
  #                   ymin = -2, ymax = 0)


ACRU_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="ACRU",]),terms="all.neighbors")
ACRU_biomassRGR_lmer_plot<-ggplot(ACRU_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="ACRU",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="ACRU",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Acer rubrum")+
  guides(color=guide_legend(title="Treatment"))

BEPA_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="BEPA",]),terms="all.neighbors")
BEPA_biomassRGR_lmer_plot<-ggplot(BEPA_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="BEPA",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="BEPA",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=BEPA_Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Betula papyrifera")+
  guides(color=F)

QUAL_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="QUAL",]),terms="all.neighbors")
QUAL_biomassRGR_lmer_plot<-ggplot(QUAL_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="QUAL",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUAL",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Quercus alba")+
  guides(color=F)

QUEL_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="QUEL",]),terms="all.neighbors")
QUEL_biomassRGR_lmer_plot<-ggplot(QUEL_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="QUEL",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUEL",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Quercus ellipsoidalis")+
  guides(color=F)

QUMA_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="QUMA",]),terms="all.neighbors")
QUMA_biomassRGR_lmer_plot<-ggplot(QUMA_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="QUMA",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QUMA",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Quercus macrocarpa")+
  guides(color=F)

QURU_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~all.neighbors+(1|Plot),
                                     data=FABdata[FABdata$Species=="QURU",]),terms="all.neighbors")
QURU_biomassRGR_lmer_plot<-ggplot(QURU_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="QURU",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="QURU",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Quercus rubra")+
  guides(color=F)

TIAM_seq<-seq(from=floor(min(FABdata$all.neighbors[FABdata$Species=="TIAM"],na.rm=T)),
              to=ceiling(max(FABdata$all.neighbors[FABdata$Species=="TIAM"],na.rm=T)),
              by=1)

TIAM_biomassRGR_lmer<-ggpredict(lmer(RGR16.18~log(all.neighbors)+(1|Plot),
                                     data=FABdata[FABdata$Species=="TIAM" & FABdata$all.neighbors!=0,]),
                                terms="all.neighbors [TIAM_seq]")

TIAM_biomassRGR_lmer_plot<-ggplot(TIAM_biomassRGR_lmer,aes(x=x,y=predicted))+
  geom_point(data=FABdata[FABdata$Species=="TIAM",],
             aes(x=all.neighbors,y=RGR16.18),color="gray")+
  geom_point(data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species=="TIAM",],
             aes(x=Neighbor.biomass,y=RGR16.18plot,
                 color=Treatment),size=4)+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high),alpha=0.15)+
  geom_line(color="red",size=1)+
  xlab("Mean volume of neighbors")+
  ylab(expression(paste("Relative growth rate (",year^-1,")")))+
  coord_cartesian(ylim=c(-2,2.5),xlim=c(0,9000))+theme_bw()+
  theme(text=element_text(size=25),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.8,0.8),
        legend.title = element_text(size = 15), 
        legend.text = element_text(size = 15),
        plot.title = element_text(face="italic"),
        legend.background = element_rect(fill=alpha('blue', 0)))+
  ggtitle("Tilia americana")+
  guides(color=F)

pdf("../LTER_ASM_2018/biomassRGR.pdf",width=12,height=20,onefile = F)
egg::ggarrange(plots = list(ACNE_biomassRGR_lmer_plot,ACRU_biomassRGR_lmer_plot,
                            BEPA_biomassRGR_lmer_plot,QUAL_biomassRGR_lmer_plot,
                            QUEL_biomassRGR_lmer_plot,QUMA_biomassRGR_lmer_plot,
                            QURU_biomassRGR_lmer_plot,TIAM_biomassRGR_lmer_plot),
               nrow=4,ncol=2)
dev.off()

##############################
## compile statistical results

## for two-lines test
## download R script by Uri Simonsohn at:
## http://webstimate.org/twolines/
source("twolines.R")
## replace with other species codes as needed
twolines(RGR16.18~all.neighbors,data=FABdata[FABdata$Species=="TIAM",])

## everything from this point on will be written to a file
sink("biomass_stats.txt")
bl_sp<-c("ACNE","ACRU","BEPA","QUAL","QUEL","QUMA","QURU","TIAM")

print("######## RGR16.18~neighbor biomass+(1|Plot)")
for(i in bl_sp){
  print(i)
  lmer.obj<-lmer(RGR16.18~all.neighbors+(1|Plot),data=FABdata[FABdata$Species==i,])
  print(summary(lmer.obj))
  print(r.squaredGLMM(lmer.obj))
}

print("######## RGR16.18~neighbor biomass (no edge)+(1|Plot)")
for(i in bl_sp){
  print(i)
  lmer.obj<-lmer(RGR16.18~all.neighbors+(1|Plot),data=FABdata.noedge[FABdata.noedge$Species==i,])
  print(summary(lmer.obj))
  print(r.squaredGLMM(lmer.obj))
}

print("######## RGR16.18~log(neighbor biomass)+(1|Plot)")
for(i in bl_sp){
  print(i)
  lmer.obj<-lmer(RGR16.18~log(all.neighbors)+(1|Plot),
                   data=FABdata[FABdata$Species==i & FABdata$all.neighbors!=0,])
  print(summary(lmer.obj))
  print(r.squaredGLMM(lmer.obj))
}

print("######## RGR16.18~neighbor biomass+2016 size")
for(i in bl_sp){
  print(i)
  lmer.obj<-lmer(RGR16.18~all.neighbors+V_2016_proj+(1|Plot),data=FABdata[FABdata$Species==i,])
  print(summary(lmer.obj))
  print(r.squaredGLMM(lmer.obj))
}


print("######## RGR16.18~log(neighbor biomass)+2016 size")
for(i in bl_sp){
  print(i)
  lmer.obj<-lmer(RGR16.18~log(all.neighbors)+V_2016_proj+(1|Plot),
                 data=FABdata[FABdata$Species==i & FABdata$all.neighbors!=0,])
  print(summary(lmer.obj))
  print(r.squaredGLMM(lmer.obj))
}

print("######## PlotRGR16.18~neighbor biomass")
for(i in bl_sp){
  print(i)
  print(summary(lm(RGR16.18plot~Neighbor.biomass,data=FAB.sp.neighbor.ag[FAB.sp.neighbor.ag$Species==i,])))
}

print("######## Absolute growth~neighbor biomass+(1|Plot)")
for(i in bl_sp){
  print(i)
  lmer.obj<-rlmer(V_growth_proj~all.neighbors+(1|Plot),data=FABdata[FABdata$Species==i,])
  print(summary(lmer.obj))
  var.fix <- get_variance_fixed(lmer.obj)
  var.ran <- get_variance_random(lmer.obj)
  var.res <- get_variance_residual(lmer.obj)
  R2m = var.fix/(var.fix+var.ran+var.res)
  R2c = (var.fix+var.ran)/(var.fix+var.ran+var.res)
  print(c(R2m,R2c))
}

print("######## Absolute growth~neighbor biomass (no mortality)+(1|Plot)")
for(i in bl_sp){
  print(i)
  lmer.obj<-rlmer(V_growth_proj~all.neighbors+(1|Plot),data=FABdata[FABdata$Species==i & FABdata$V_2018_proj>0,])
  print(summary(lmer.obj))
  var.fix <- get_variance_fixed(lmer.obj)
  var.ran <- get_variance_random(lmer.obj)
  var.res <- get_variance_residual(lmer.obj)
  R2m = var.fix/(var.fix+var.ran+var.res)
  R2c = (var.fix+var.ran)/(var.fix+var.ran+var.res)
  print(c(R2m,R2c))
}

sink()