######################################################################################
# R-code to reproduce results presented in:
# Vicari et al. 
# "Unpacking multi-trophic herbivore-grass-endophyte interactions: feedbacks across
# different scales in vegetation responses to Soay sheep herbivory"
#
# In review at Science of Nature
#
# Contact: Gustaf.Granath@gmail.com
######################################################################################

# packages used in most sections, load first
library(ggplot2)
library(gplots)
library(nlme)
library(car)

# sheep population ####
# read data here
# extracted from http://soaysheep.biology.ed.ac.uk/population-ecology
# while island population size
dat <-structure(list(year = 1985:2010, population = c(1297.29401714987, 
                                                706.645692226542, 1038.5615392874, 1467.74713040718, 696.936455848648, 
                                                891.896055155931, 1461.62284284574, 953.867123498686, 1289.367919376, 
                                                1520.35438712779, 1178.38387893329, 1834.59755682004, 1755.68200383702, 
                                                1972.25398988932, 931.171783465357, 1407.2231117169, 1890.44433760135, 
                                                914.236634629299, 1577.64821756158, 1996.06029062359, 1372.95137446378, 
                                                1805.77792922527, 1543.06839876954, 1918.21835512466, 2224.86218019036, 
                                                2063.08202904369)), .Names = c("year", "population"), class = "data.frame", row.names = c(NA, 
                                                                                                                                          -26L))
ggplot(dat, aes(y=population, x =year)) +
  geom_path() +
  geom_point(data=dat[dat$year == 1991 |dat$year == 1992|dat$year == 1993,], aes(y=population, x =year), shape = 21, size=3) +
  xlab("Year") +
  ylab("Population size") +
  scale_y_continuous( limits=c(0, 2500), breaks=seq(0,3000,500)) +
  theme_bw(base_size = 12, base_family = "Helvetica") +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_line(color="white"), axis.line = element_line(colour = "black"),
        axis.text = element_text(size=12, color = "black"),
        axis.title = element_text(size=15))
ggsave("figure.pdf", height=3, width = 4)


# Standing crop and cum ANPP ####
dat <- read.csv("vicari_etal_prod_91-93.csv",sep=",")

tt<-aggregate(grass.dicots ~ Date+year+plot +site+treatment+exclosure, dat ,mean)
tt$Date<-as.Date(tt$Date, "%d/%m/%Y")

#analysis 1991 standing crop
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year=="1991")
tt2<-droplevels(tt2)
mod<-lme(log(grass.dicots)~poly(Date,2)*site*treatment,random= ~1|exclosure/plot,tt2)
summary(mod)
plot(mod)
anova(mod,type="marginal")

#analysis 1992 standing crop
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year=="1992")
tt2<-droplevels(tt2)
mod<-lme(log(grass.dicots)~poly(Date,2)*site*treatment,random= ~1|exclosure/plot,tt2)
summary(mod)
anova(mod,type="marginal")
plot(mod)

#analysis 1993  standing crop
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year=="1993")
tt2<-droplevels(tt2)
mod<-lme(log(grass.dicots)~poly(Date,2)*treatment,random= ~1|exclosure/plot,tt2)
summary(mod)
plot(mod)
anova(mod,type="marginal")

# between year analysis standing crop
# both sites, 1991-1992
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year=="1991" | year=="1992")
tt2<-droplevels(tt2)
tt2$ex.yr<-paste(tt2$exclosure,tt2$year,sep="")
mod<-lme(grass.dicots~factor(year)*treatment*site,random= ~1|ex.yr/plot,tt2,weights=varExp())
summary(mod)
plot(mod)
anova(mod,type="marginal")

  # only VM site
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  site=="VM")
tt2<-droplevels(tt2)
tt2$ex.yr<-paste(tt2$exclosure,tt2$year,sep="")
mod<-lme(grass.dicots~factor(year)*treatment,random= ~1|ex.yr/plot,tt2,weights=varPower())
summary(mod)
plot(mod)
anova(mod,type="marginal")

# cumulative ANPP
tt<-aggregate(cum.anpp ~ Date+year+plot +site+treatment+exclosure, dat ,mean)
tt$Date<-as.Date(tt$Date, "%d/%m/%Y")

# analyses 1991
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year=="1991")
tt2<-droplevels(tt2)
mod<-lme(cum.anpp~poly(Date,2)*site*treatment,random= ~1|exclosure/plot,tt2,
         control = lmeControl(msMaxIter = 1000,opt="optim"),weights=varExp() )
summary(mod)
plot(mod)
anova(mod,type="marginal")

# analyses 1992
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year=="1992")
tt2<-droplevels(tt2)
mod<-lme(log(cum.anpp+1)~poly(Date,2)*site*treatment,random= ~1|exclosure/plot,tt2)
summary(mod)
anova(mod,type="marginal")
plot(mod)

#analyses 1993
tt2<- subset(tt,  treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year=="1993")
tt2<-droplevels(tt2)
mod<-lme(cum.anpp~poly(Date,2)*treatment,random= ~1|plot,tt2)
summary(mod)
plot(mod)
anova(mod,type="marginal")

# between year analyses, all sites at last date of the season
tt2<- subset(dat,  treatment=="ungrazed" | treatment=="grazed")
tt2$Date2 <- strptime(tt2$Date, "%d/%m/%Y")$yday

tt3<- subset(tt2,  Date2 > 247 & year=="1993")
tt2<- subset(tt2,  Date2 > 256)
tt2<- rbind(tt2,tt3)
tt2<-droplevels(tt2)

ttt<-aggregate(cum.anpp~ year+ Date + exclosure + site+treatment,tt2,mean)
ave<-aggregate(cum.anpp~ year+site+treatment,tt2,mean) #ger means
ttt$plot3<-paste(ttt$exclosure, ttt$year,sep="")
mod<-lme(cum.anpp~year*treatment*site,random= ~1|plot3,ttt[ttt$year != "1993",],weights=varPower() )
# no big interaction so continue without it and include year 1993

mod<-lme(cum.anpp~factor(year)*treatment + site*treatment,random= ~1|plot3,ttt,weights=varPower() )
summary(mod)
plot(mod)
anova(mod,type="marginal")

# Plot figure 2 ####
se<- function (x) sd(x)/sqrt(length(x))

tt2<- subset(dat,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year==1991)
tt2<-droplevels(tt2)
ttt<-aggregate(grass.dicots ~ Date+plot +site+treatment,tt2,mean)

ttp<-aggregate(grass.dicots ~ Date+site+treatment,ttt,mean)
ttp$ses<-aggregate(grass.dicots ~ Date+site+treatment,ttt,se)$grass.dicots
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")

ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
day91.vb<-ttp.re.vb$Date
VB_G_91<-ttp.re.vb$grass.dicots.grazed_caged
VB_UG_91<-ttp.re.vb$grass.dicots.ungrazed
VB_UG_se_91<-ttp.re.vb$ses.ungrazed  
VB_G_se_91<-ttp.re.vb$ses.grazed_caged  		

ttp.re.f <-reshape(ttp[ttp$site=="SBF",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.f<-ttp.re.f[order(ttp.re.f$Date),]
day91.f<-ttp.re.f$Date
F_G_91<-ttp.re.f$grass.dicots.grazed_caged
F_UG_91<-ttp.re.f$grass.dicots.ungrazed
F_UG_se_91<-ttp.re.f$ses.ungrazed  
F_G_se_91<-ttp.re.f$ses.grazed_caged  		

#1992
tt2<- subset(dat,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year==1992)
tt2<-droplevels(tt2)
ttt<-aggregate(grass.dicots ~ Date+plot +site+treatment,tt2,mean)

ttp<-aggregate(grass.dicots ~ Date+site+treatment,ttt,mean)
ttp$ses<-aggregate(grass.dicots ~ Date+site+treatment,ttt,se)$grass.dicots
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")
ttp$Date<- ttp$Date-365-1
ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
day92.vb<-ttp.re.vb$Date
VB_G_92<-ttp.re.vb$grass.dicots.grazed_caged
VB_UG_92<-ttp.re.vb$grass.dicots.ungrazed
VB_UG_se_92<-ttp.re.vb$ses.ungrazed  
VB_G_se_92<-ttp.re.vb$ses.grazed_caged  		

ttp.re.f <-reshape(ttp[ttp$site=="SBF",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.f<-ttp.re.f[order(ttp.re.f$Date),]
day92.f<-ttp.re.f$Date
F_G_92<-ttp.re.f$grass.dicots.grazed_caged
F_UG_92<-ttp.re.f$grass.dicots.ungrazed
F_UG_se_92<-ttp.re.f$ses.ungrazed  
F_G_se_92<-ttp.re.f$ses.grazed_caged  		

#1993
tt2<- subset(dat,  treatment=="ungrazed" | treatment=="grazed_caged")
tt2<- subset(tt2,  year==1993)
tt2<-droplevels(tt2)
ttt<-aggregate(grass.dicots ~ Date+plot +site+treatment,tt2,mean)

ttp<-aggregate(grass.dicots ~ Date+site+treatment,ttt,mean)
ttp$ses<-aggregate(grass.dicots ~ Date+site+treatment,ttt,se)$grass.dicots
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")
ttp$Date<- ttp$Date-2*365-1
ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
day93.vb<-ttp.re.vb$Date
VB_G_93<-ttp.re.vb$grass.dicots.grazed_caged
VB_UG_93<-ttp.re.vb$grass.dicots.ungrazed
VB_UG_se_93<-ttp.re.vb$ses.ungrazed  
VB_G_se_93<-ttp.re.vb$ses.grazed_caged  		


# fig 2a - V. Meadow
plotCI(
  x=day91.vb ,
  y=VB_G_91,
  uiw=VB_G_se_91, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="blue", # fill colour for symbol
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  col="blue",
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  bty="l",
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  
  xlim=c(as.Date("1991/4/15"), as.Date("1991/10/1")), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(0,1250),
  #		yaxp=c(-4,8,6),
  ylab=expression(paste("Standing crop ( g",m^{-2},")",sep="")),
  cex.lab=1.3,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
)

Sys.setlocale(category = "LC_TIME", locale="English") 
axis.Date(1, at = seq(as.Date("1991/4/1"), as.Date("1991/10/1"), "month"),format="%b %Y")
Sys.setlocale(category = "LC_TIME", locale="") 

lines(x=day91.vb,y=VB_UG_91,lty=2,col="blue")
plotCI(
  x=day91.vb ,
  y=VB_UG_91,
  uiw=VB_UG_se_91, # error bar length
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  col="blue",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day92.vb ,
  y=VB_G_92,
  uiw=VB_G_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="red",
  col="red",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day92.vb,y=VB_UG_92,lty=2,col="red")
plotCI(
  x=day92.vb ,
  y=VB_UG_92,
  uiw=VB_UG_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="white",
  col="red",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day93.vb ,
  y=VB_G_93,
  uiw=VB_G_se_93, # error bar length
  
  pch=24, # symbol type 21 = filled circle
  pt.bg="black",
  col="black",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day93.vb,y=VB_UG_93,lty=2,col="black")
plotCI(
  x=day93.vb ,
  y=VB_UG_93,
  uiw=VB_UG_se_93, # error bar length 
  
  pch=24, # symbol type 21 = filled circle
  pt.bg="white",
  col="black",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plot.new()
legend(
  "topleft",
  #x=7889, # x coordinate of the top left of the legend
  #y=1520, # y coordinate of the top left of the legend
  box.lty=0, # line type to surround the legend box (0 for none)
  
  
  legend=c("1991 Grazed","1991 Ungrazed","1992 Grazed","1992 Ungrazed","1993 Grazed","1993 Ungrazed"), # sequence of text for the legend
  pch=c(21,21,22,22,24,24), # sequence of point types for the legend; -1 is a nonexistent point
  pt.bg=c("blue", "white", "red","white","black","white"), # sequence of fill colours for the points
  col=c("blue","blue","red","red","black","black"),
  pt.cex=c(1.3),
  lty=c(1,2,1,2,1,2), # sequence of line types for the legend
  lwd=1.1,
  xpd=TRUE,
  cex=1.3
)

#fig 2b - F-banks
plotCI(
  x=day91.f ,
  y=F_G_91,
  uiw=F_G_se_91, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="blue", # fill colour for symbol
  col="blue",
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  bty="l",
  xlim=c(as.Date("1991/4/15"), as.Date("1991/10/1")), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(0,1250),
  #		yaxp=c(-4,8,6),
  ylab=expression(paste("Standing crop ( g",m^{-2},")",sep="")),
  cex.lab=1.3,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
)
Sys.setlocale(category = "LC_TIME", locale="English") 
axis.Date(1, at = seq(as.Date("1991/4/1"), as.Date("1991/10/1"), "month"),format="%b %Y")
Sys.setlocale(category = "LC_TIME", locale="") 

lines(x=day91.f,y=F_UG_91,lty=2,col="blue")
plotCI(
  x=day91.f ,
  y=F_UG_91,
  uiw=F_UG_se_91, # error bar length 
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  col="blue",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day92.f ,
  y=F_G_92,
  uiw=F_G_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="red",
  col="red",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day92.f,y=F_UG_92,lty=2,col="red")
plotCI(
  x=day92.f ,
  y=F_UG_92,
  uiw=F_UG_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="white",
  col="red",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

legend("topleft",
       #x=7787, # x coordinate of the top left of the legend
       #y=1870, # y coordinate of the top left of the legend
       box.lty=0, # line type to surround the legend box (0 for none)
       
       legend=c("1991 Grazed","1992 Grazed","1991 Ungrazed","1992 Ungrazed"), # sequence of text for the legend
       pch=c(21,22,21,22), # sequence of point types for the legend; -1 is a nonexistent point
       pt.bg=c("black", "black","white","white"), # sequence of fill colours for the points
       pt.cex=c(1.1),
       lty=c(1,2,1,2,1,2), # sequence of line types for the legend
       xpd=TRUE
)

# Plot Figure 3 ####
se<- function (x) sd(x)/sqrt(length(x))

#1991
tt2<- subset(dat, treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year==1991)
tt2<-droplevels(tt2)
ttt<-aggregate(cum.anpp ~ Date+plot +site+treatment,tt2,mean)

ttp<-aggregate(cum.anpp ~ Date+site+treatment,ttt,mean)
ttp$ses<-aggregate(cum.anpp ~ Date+site+treatment,ttt,se)$cum.anpp
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")
ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
ttp.re.vb[1,3:4]<-0
day91.vb<-ttp.re.vb$Date
VB_G_91<-ttp.re.vb$cum.anpp.grazed
VB_UG_91<-ttp.re.vb$cum.anpp.ungrazed
VB_UG_se_91<-ttp.re.vb$ses.ungrazed  
VB_G_se_91<-ttp.re.vb$ses.grazed  		

ttp.re.f <-reshape(ttp[ttp$site=="SBF",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.f<-ttp.re.f[order(ttp.re.f$Date),]
ttp.re.f[1,3:4]<-0
day91.f<-ttp.re.f$Date
F_G_91<-ttp.re.f$cum.anpp.grazed
F_UG_91<-ttp.re.f$cum.anpp.ungrazed
F_UG_se_91<-ttp.re.f$ses.ungrazed  
F_G_se_91<-ttp.re.f$ses.grazed  

#1992
tt2<- subset(dat, treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year==1992)
tt2<-droplevels(tt2)
ttt<-aggregate(cum.anpp ~ Date+plot +site+treatment,tt2,mean)

ttp<-aggregate(cum.anpp ~ Date+site+treatment,ttt,mean)
ttp$ses<-aggregate(cum.anpp ~ Date+site+treatment,ttt,se)$cum.anpp
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")
ttp$Date<- ttp$Date-365-1
ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
ttp.re.vb[1,3:4]<-0
day92.vb<-ttp.re.vb$Date
VB_G_92<-ttp.re.vb$cum.anpp.grazed
VB_UG_92<-ttp.re.vb$cum.anpp.ungrazed
VB_UG_se_92<-ttp.re.vb$ses.ungrazed  
VB_G_se_92<-ttp.re.vb$ses.grazed  		

ttp.re.f <-reshape(ttp[ttp$site=="SBF",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.f<-ttp.re.f[order(ttp.re.f$Date),]
ttp.re.f[1,3:4]<-0
day92.f<-ttp.re.f$Date
F_G_92<-ttp.re.f$cum.anpp.grazed
F_UG_92<-ttp.re.f$cum.anpp.ungrazed
F_UG_se_92<-ttp.re.f$ses.ungrazed  
F_G_se_92<-ttp.re.f$ses.grazed  

#1993
tt2<- subset(dat, treatment=="ungrazed" | treatment=="grazed")
tt2<- subset(tt2,  year==1993)
tt2<-droplevels(tt2)

ttp<-aggregate(cum.anpp ~ Date+site+treatment,tt2,mean)
ttp$ses<-aggregate(cum.anpp ~ Date+site+treatment,tt2,se)$cum.anpp
ttp$Date<- as.Date(ttp$Date, "%d/%m/%Y")
ttp$Date<- ttp$Date-2*365-1
ttp.re.vb <-reshape(ttp[ttp$site=="VM",],direction="wide",timevar="treatment",idvar="Date")
ttp.re.vb<-ttp.re.vb[order(ttp.re.vb$Date),]
ttp.re.vb[1,3:4]<-0
day93.vb<-ttp.re.vb$Date
VB_G_93<-ttp.re.vb$cum.anpp.grazed
VB_UG_93<-ttp.re.vb$cum.anpp.ungrazed
VB_UG_se_93<-ttp.re.vb$ses.ungrazed  
VB_G_se_93<-ttp.re.vb$ses.grazed  		

plotCI(
  x=day91.vb ,
  y=VB_G_91,
  uiw=VB_G_se_91, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="blue", # fill colour for symbol
  col="blue",
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  bty="l",
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  
  xlim=c(as.Date("1991/4/1"), as.Date("1991/10/15")), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(0,1800),
  #		yaxp=c(-4,8,6),
  ylab=bquote("Cumulative ANPP ("*g*~m^-2*")"),
  cex.lab=1.3,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
  
)

Sys.setlocale(category = "LC_TIME", locale="English") 
axis.Date(1, at = seq(as.Date("1991/4/1"), as.Date("1991/10/1"), "month"),format="%b")
Sys.setlocale(category = "LC_TIME", locale="") 

lines(x=day91.vb,y=VB_UG_91,lty=2,col="blue")
plotCI(
  x=day91.vb ,
  y=VB_UG_91,
  uiw=VB_UG_se_91, # error bar length 
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  col="blue",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day92.vb ,
  y=VB_G_92,
  uiw=VB_G_se_92, # error bar length
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="red",
  col="red",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day92.vb,y=VB_UG_92,lty=2,col="red")
plotCI(
  x=day92.vb ,
  y=VB_UG_92,
  uiw=VB_UG_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="white",
  col="red",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day93.vb ,
  y=VB_G_93,
  uiw=VB_G_se_93, # error bar length 
  
  pch=24, # symbol type 21 = filled circle
  pt.bg="black",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day93.vb,y=VB_UG_93,lty=2)
plotCI(
  x=day93.vb ,
  y=VB_UG_93,
  uiw=VB_UG_se_93, # error bar length (default is to put this much above and below point)
  
  pch=24, # symbol type 21 = filled circle
  pt.bg="white",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)
plot.new()
legend(
  "topleft",
  #x=7889, # x coordinate of the top left of the legend
  #y=1520, # y coordinate of the top left of the legend
  box.lty=0, # line type to surround the legend box (0 for none)
  
  
  legend=c("1991 Ungrazed","1991 Grazed","1992 Ungrazed","1992 Grazed","1993 Ungrazed","1993 Grazed"), # sequence of text for the legend
  pch=c(21,21,22,22,24,24), # sequence of point types for the legend; -1 is a nonexistent point
  pt.bg=c( "white","blue", "white","red","white","black"), # sequence of fill colours for the points
  col=c("blue","blue","red","red","black","black"),
  pt.cex=c(1.3),
  lty=c(2,1,2,1,2,1), # sequence of line types for the legend
  lwd=1.1,
  xpd=TRUE,
  cex=1.3
)

#####fig 1b - F-banks
plotCI(
  x=day91.f ,
  y=F_G_91,
  uiw=F_G_se_91, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="blue", # fill colour for symbol
  col="blue",
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  bty="l",
  xlim=c(as.Date("1991/4/1"), as.Date("1991/10/15")), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(0,1800),
  #		yaxp=c(-4,8,6),
  ylab=bquote("Cumulative ANPP ("*g*~m^-2*")"),
  cex.lab=1.3,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
)
Sys.setlocale(category = "LC_TIME", locale="English") 
axis.Date(1, at = seq(as.Date("1991/4/1"), as.Date("1991/10/1"), "month"),format="%b")
Sys.setlocale(category = "LC_TIME", locale="") 

lines(x=day91.f,y=F_UG_91,lty=2,col="blue")
plotCI(
  x=day91.f ,
  y=F_UG_91,
  uiw=F_UG_se_91, # error bar length 
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  col="blue",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

plotCI(
  x=day92.f ,
  y=F_G_92,
  uiw=F_G_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="red",
  col="red",
  cex=1,
  lty=1,
  type="o",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

lines(x=day92.f,y=F_UG_92,lty=2,col="red")
plotCI(
  x=day92.f ,
  y=F_UG_92,
  uiw=F_UG_se_92, # error bar length 
  
  pch=22, # symbol type 21 = filled circle
  pt.bg="white",
  col="red",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
)

legend("topleft",
       #x=7787, # x coordinate of the top left of the legend
       #y=1870, # y coordinate of the top left of the legend
       box.lty=0, # line type to surround the legend box (0 for none)
       
       legend=c("1991 Grazed","1992 Grazed","1991 Ungrazed","1992 Ungrazed"), # sequence of text for the legend
       pch=c(21,22,21,22), # sequence of point types for the legend; -1 is a nonexistent point
       pt.bg=c("black", "black","white","white"), # sequence of fill colours for the points
       pt.cex=c(1.1),
       lty=c(1,2,1,2,1,2), # sequence of line types for the legend
       xpd=TRUE
)



# Relative growth rate - RGR ####
dat<-read.csv("vicari_etal_rgr.csv")
dat$rand <- interaction(dat$Exclosure,dat$Year)

mod <- lme(RGR ~ factor(Year)*Treatment * Site, random = ~ 1|rand, data=dat[dat$Year < 1993,])
summary(mod)
anova(mod,type="marginal")
plot(mod)
# no 3-way interaction effect
# we fit a model with all three years without the 3-way interaction (impossible as 1993 only has one site)
mod <- lme(RGR ~ factor(Year)*Treatment + Site, random = ~ 1|rand, data=dat)
summary(mod)
anova(mod,type="marginal")
plot(mod)

# Nutrient (N,P) application field experiment ####
library(tidyr)
library(dplyr)
library(ggplot2)

dat <- read.csv("vicari_etal_NP_exp.csv",sep=",")

# Type 3 test
options(contrasts = c("contr.sum", "contr.poly"))
# Run and show model results
mod<-(lm(log(biomass_no_moss_gm2_11Aug)~factor(P)*factor(N)*site, dat))
Anova(mod,type=3)
# Re-set contrasts
options(contrasts = c("contr.treatment", "contr.poly"))


# individual sites
mod<-(lm(log(biomass_no_moss_gm2_11Aug)~factor(P)*factor(N), dat,subset=site=="Village Meadows"))
summary.aov(mod)
mod<-(lm(log(biomass_no_moss_gm2_11Aug)~N*P, dat,subset=site=="St. Brianan Fanks"))
summary.aov(mod)
plot(mod)

# plot data
data_long <- gather(dat, date, biomass.g.sqm, c(biomass_no_moss_gm2_8May,
                                                biomass_no_moss_gm2_12June, 
                                                biomass_no_moss_gm2_11Aug), factor_key=TRUE)
levels(data_long$date)[levels(data_long$date)=="biomass_no_moss_gm2_8May"] <- "1992-05-08"
levels(data_long$date)[levels(data_long$date)=="biomass_no_moss_gm2_12June"] <- "1992-06-12"
levels(data_long$date)[levels(data_long$date)=="biomass_no_moss_gm2_11Aug"] <- "1992-08-11"

sum.dat <- data_long %>% select(site,Treatment,date,biomass.g.sqm) %>% group_by(site,Treatment,date) %>% 
  summarise_all(funs(mean(., na.rm=TRUE), sd = sd(.), samp=length(.)))
sum.dat <- sum.dat  %>% mutate(se = sd/sqrt(samp))
sum.dat <- sum.dat  %>% mutate(date = as.POSIXct(strptime(date, format = "%Y-%m-%d")))
sum.dat$low <- sum.dat$mean-sum.dat$se
sum.dat$high <- sum.dat$mean+sum.dat$se

# add noise for plotting so error bars are visual
dd <- 24*60*60
sum.dat$date <- sum.dat$date + rep(c(0,0,0,-dd*3,-dd*3,-dd*3, +dd*3,+dd*3,+dd*3, +dd*1.5,+dd*1.5,+dd*1.5),2)

ggplot(sum.dat, aes(x=date, y=mean, group = Treatment, shape = Treatment) ) +
  geom_point(size = 2) +
  geom_path() +
    geom_errorbar(data=sum.dat[!(is.na(sum.dat$se)),], aes(ymin=low, ymax=high), width=0.001)+
  ylab(bquote("Standing crop g (" * m^-2*")")) +
  facet_wrap(~site) +
  theme_bw()

##### F. rubra leaf demography ####
dat<-read.csv("vicari_etal_demogr_Frubra_field.csv",na.strings = "",sep=",")
library(nlme)
se <- function(x) sd(x)/sqrt(length(x))

# leaf births
new.sum<-aggregate(Cum_leaf_births~Site+Treatment+Quadrat_pair,dat,mean)
str(new.sum)

mod<-lme(Cum_leaf_births~Site*Treatment,random=~1|factor(Quadrat_pair),new.sum)
anova(mod)
plot(mod)


aggregate(Cum_leaf_births~Treatment,new.sum,mean)
aggregate(Cum_leaf_births~Treatment,new.sum,se)

xtabs(~Site+Treatment+Quadrat_pair,dat)

# leaf deaths
new.sum<-aggregate(Cum_leaf_deaths~Site+Treatment+Quadrat_pair,dat,mean)
mod<-lme(Cum_leaf_deaths~Site*Treatment,random=~1|factor(Quadrat_pair),new.sum)
anova(mod)
plot(mod)

aggregate(Cum_leaf_deaths~Treatment,new.sum,mean)
aggregate(Cum_leaf_deaths~Treatment,new.sum,se)

# life span
dat$life.ave<-apply(cbind(dat$Lifespan_leaf_1,dat$Lifespan_leaf_2),1,mean)
new.sum<-aggregate(life.ave~Site+Treatment+Quadrat_pair,dat,mean)

mod<-lme(life.ave~Site*Treatment,random=~1|factor(Quadrat_pair),new.sum)
anova(mod)
plot(mod)

aggregate(life.ave ~ Treatment,new.sum,mean)
aggregate(life.ave ~ Treatment,new.sum,se)

new.sum2<-aggregate(Lifespan_leaf_1~Site+Treatment+Quadrat_pair+leaf_1_status,dat,mean)
new.sum2<-cbind(new.sum2, aggregate(Lifespan_leaf_2~Site+Treatment+Quadrat_pair+leaf_2_status,dat,mean)[,4:5])
new.sum2$new.ave<-apply(cbind(new.sum2$Lifespan_leaf_1,new.sum2$Lifespan_leaf_2),1,mean)
new.sum.fin<-aggregate(new.ave~Treatment+leaf_1_status,new.sum2, function (x) c(mean(x),se(x)))

# F. rubra leaf demography - no of tillers ####

dat <- read.csv("vicari_etal_no_tillers_Frubra.csv",sep=",")
library(lme4)
se <- function (x) sd(x)/sqrt(length(x))

dat2<-reshape(dat,direction="long",varying=list(5:16),times=colnames(dat)[5:16],v.names="tiller")
dat2<-dat2[!(is.na(dat2$tiller)),]
dat2$time2<-(substr(dat2$time,2,11) )
dat2<-transform(dat2, Treatment=factor(Treatment), Plot.pair=factor(Plot.pair),time2=ordered(time2))

dat.agg<-aggregate(tiller~time2+Site+Plot.pair+Treatment, dat2, mean)
dat.agg$id<-paste(dat.agg$Plot.pair, dat.agg$Treatment,paste="")

xtabs(~time2 +Treatment, dat.agg)

# seasonal means
aggregate(tiller ~ Treatment, dat.agg, FUN=mean)

options(contrasts = c("contr.sum", "contr.poly"))
# Run and show model results
mod<-lmer(tiller~Treatment*Site*time2 +( 1|Plot.pair/id),data=dat.agg)
Anova(mod,type=3, test.statistic = "F")
# Re-set contrasts
options(contrasts = c("contr.treatment", "contr.poly"))

# Plot Figure 6
dat2<-transform(dat2, time=as.Date(time2,"%Y.%m.%d"))

mean2 <- aggregate(tiller ~ time2 + Treatment, dat.agg, FUN=mean)
mean2$ses <- aggregate(tiller ~ time2 + Treatment, dat.agg, FUN=se)[,3]

mean2<-transform(mean2, time=as.Date(time2,"%Y.%m.%d"))
png("FIG6_png.png", width = 580, height = 480, type="cairo")
with(mean2[mean2$Treatment=="grazed",], plotCI(
  x = time ,
  y = tiller,
  uiw=ses, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="black", # fill colour for symbol
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  bty="l",
  xlim=c(as.Date("1993/05/01"), as.Date("1993/10/15")), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(2,5),
  #		yaxp=c(-4,8,6),
  ylab="Number of leaves per tiller",
  cex.lab=1.3,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
))
Sys.setlocale(category = "LC_TIME", locale="English") 
axis.Date(1, at = seq(as.Date("1993/5/1"), as.Date("1993/10/15"), "month"),format="%b %Y")
Sys.setlocale(category = "LC_TIME", locale="") 

with(mean2[mean2$Treatment=="ungrazed",], lines(x=time,y=tiller,lty=2))

with(mean2[mean2$Treatment=="ungrazed",], plotCI(
  x=time ,
  y=tiller,
  uiw=ses, # error bar length 
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
))

legend("topleft",
       #x=7787, # x coordinate of the top left of the legend
       #y=1870, # y coordinate of the top left of the legend
       box.lty=0, # line type to surround the legend box (0 for none)
       
       legend=c("Grazed","Ungrazed"), # sequence of text for the legend
       cex=1.2,
       pch=c(21,21), # sequence of point types for the legend; -1 is a nonexistent point
       pt.bg=c("black","white"), # sequence of fill colours for the points
       pt.cex=c(1.2),
       lty=c(1,2),lwd=1.2, # sequence of line types for the legend
       xpd=TRUE
)
dev.off()



#### Leaf demography of R. acris and L. autumnalis ####
dat<-read.csv("vicari_etal_demogr_Racr_Laut_field.csv",na.strings = "",sep=",")
se<- function (x) sd(x)/sqrt(length(x))

# t-tests and means for L. autumnalis
with(dat[dat$species=="L_autumnalis",], t.test(Cum_leaves~treatment ,paired=TRUE,var.equal=TRUE))
with(dat[dat$species=="L_autumnalis",], t.test(leaf_mm~treatment ,paired=TRUE))
aggregate(Cum_leaves~treatment, dat[dat$species=="L_autumnalis",],mean )
aggregate(Cum_leaves~treatment, dat[dat$species=="L_autumnalis",],se )

# t-tests and means for R. acris
with(dat[dat$species=="R_acris",], t.test(Cum_leaves~treatment ,paired=TRUE,var.equal=FALSE))
with(dat[dat$species=="R_acris",], t.test(blade_mm~treatment ,paired=TRUE))
with(dat[dat$species=="R_acris",], t.test(petiole_mm~treatment ,paired=TRUE))
aggregate(petiole_mm~treatment, dat[dat$species=="R_acris",],mean )
aggregate(petiole_mm~treatment, dat[dat$species=="R_acris",],se )
aggregate(blade_mm~treatment, dat[dat$species=="R_acris",],mean )
aggregate(blade_mm~treatment, dat[dat$species=="R_acris",],se )


#### Endophytes ####

# hyphal load in the field
dat<-read.csv("vicari_etal_hyphal_load.csv",sep=",")
se<- function (x) sd(x)/sqrt(length(x))

summary.aov(lm(log(red.green_value)~grazing, dat,subset=time_period=="1992 amaj-jun"))

means <- aggregate(red.green_value~grazing+time_period, dat,mean)
means$ses <- aggregate(red.green_value~grazing+time_period, dat, se)$red.green_value

with(means[means$grazing=="yes",], plotCI(
  x=as.numeric(time_period)-0.02 ,
  y=red.green_value,
  uiw=ses, # error bar length (default is to put this much above and below point)
  
  pch=21, # symbol (plotting character) type: see help(pch); 24 = filled triangle pointing up
  pt.bg="black", # fill colour for symbol
  col="black",
  cex=1, # symbol size multiplier
  lty=1, # line type: see help(par)
  type="o", # p=points, l=lines, b=both, o=overplotted points/lines, etc.; see help(plot.default)
  bty="l",
  gap=0, # distance from symbol to error bar
  sfrac=0, # width of error bar as proportion of x plotting region (default 0.01)
  
  xlim=c(0.5,4.5), # x axis limits
  #		xaxp=c(1,12,11), # x-min tick mark, x-max tick mark, number of intervals between tick marks
  xlab="", # x axis label
  ylim=c(0,200),
  #		yaxp=c(-4,8,6),
  ylab="Proxy for hyphal load",
  cex.lab=1.4,
  las=1, # axis labels horizontal (default is 0 for always parallel to axis)
  #font.lab=2,
  xaxt = "n"# 1 plain, 2 bold, 3 italic, 4 bold italic, 5 symbol
  
))
axis(1,at=c(1,2,3,4), labels=c("Aug-Sep \n1991","May-June \n1992","Aug-Sep \n1992","Aug-Sep \n1993"),cex.axis=1.1, padj=0.5)

with(means[means$grazing=="no",], lines(x=as.numeric(time_period)+0.02 ,y=red.green_value,lty=2))

with(means[means$grazing=="no",], plotCI(
  x=as.numeric(time_period)+0.02 ,
  y=red.green_value,
  uiw=ses, # error bar length 
  
  pch=21, # symbol type 21 = filled circle
  pt.bg="white",
  col="black",
  cex=1,
  lty=1,
  type="p",
  
  gap=0,
  sfrac=0,
  
  add=TRUE # ADD this plot to the previous one
))

legend("topright",
       #x=7787, # x coordinate of the top left of the legend
       #y=1870, # y coordinate of the top left of the legend
       box.lty=0, # line type to surround the legend box (0 for none)
       
       legend=c("Grazed","Ungrazed"), # sequence of text for the legend
       cex=1.2,
       pch=c(21,21), # sequence of point types for the legend; -1 is a nonexistent point
       pt.bg=c("black","white"), # sequence of fill colours for the points
       pt.cex=c(1.2),
       lty=c(1,2),lwd=1.2, # sequence of line types for the legend
       xpd=TRUE
)


# Ergovaline levels-greenhous exp ####
dat <- read.csv("vicari_etal_ergovaline_exp.csv")

mod <- lme(ergovaline_ppb ~ordered(no_cuts_over_8_wks), random=~1|genotype,data=dat)
anova(mod)
summary(mod)
plot(mod)

mod <- lme(ergovaline_ppb ~factor(no_cuts_over_8_wks), random=~1|genotype,data=dat)
summary(mod)

##Fig 5b
dat.sum <- aggregate(ergovaline_ppb ~ no_cuts_over_8_wks, dat, function (x) c(mean(x), se(x)))

pp<-barplot(dat.sum$ergovaline_ppb[,1] , ylim=c(0,150),xlim=c(0,6),ylab="Ergovaline (ppb)",las=1,
            xlab="Cuts",cex.lab=1.4,axisname=TRUE,names.arg=c("0","1","2","4","8"))
abline(h=0)
arrows(pp,dat.sum$ergovaline_ppb[,1]-dat.sum$ergovaline_ppb[,2],  
       pp,dat.sum$ergovaline_ppb[,1]+dat.sum$ergovaline_ppb[,2],length=0)