required files:

dat<-read.csv("Thames Estuary WeBS_2017.CSV")

required packages:

library(reshape)
library(ggplot2)
library(plyr)
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
AllData<-read.csv("LondonGatewayDataFINAL_2017.CSV")
#View(AllData)
AllData$Period<-factor(AllData$Period, levels = c("Pre-works","During","Post-works"))
AllData$Date<-as.Date(AllData$Date, "%d/%m/%Y")

library(reshape)

calculate peak counts - all species combined

TotalAssmb<-cast(AllData,Period+Winter+Month~Date,sum,value="Number")
df<-melt(TotalAssmb)
df<-rename(df,c("value"="PeakAssmb"))
df<-subset(df,PeakAssmb>0)

MonthPeakAssmb<-cast(df,Period+Winter~Month,max,value="PeakAssmb")
df2<-melt(MonthPeakAssmb)
WinterPeakAssmb<-rename(df2,c("value"="PeakAssmb"))
WinterPeakAssmb<-subset(WinterPeakAssmb,is.finite(PeakAssmb))
df3<-cast(WinterPeakAssmb,Period+Winter~.,max,value="PeakAssmb")
WinterPeakAssmb<-rename(df3,c("(all)"="PeakAssmb"))

View(WinterPeakAssmb)

KeySpp<-subset(AllData,Species=="Avocet"|Species=="Black-tailed godwit"|Species=="Dunlin")
SpYrPeaks<-cast(KeySpp,Period+Winter+Month+Species~Date,sum,value="Number")
df4<-melt(SpYrPeaks)
df4<-rename(df4,c("value"="Count"))
SpMnthPks<-cast(df4,Species+Period+Winter~Month,max,value="Count")
SpYrPeaks<-cast(df4,Species+Period~Winter,max,value="Count")
SpYrPeaks<-melt(SpYrPeaks)
SpYrPeaks<-rename(SpYrPeaks,c("value"="PeakCount"))
SpYrPeaks[13,3]<-0
SpYrPeaks<-subset(SpYrPeaks,is.finite(SpYrPeaks$PeakCount))
SpYrPeaks<-merge(SpYrPeaks,WinterPeakAssmb)
SpYrPeaks$Prop<-round(SpYrPeaks$PeakCount/SpYrPeaks$PeakAssmb*100,digits=2)

tableExport<-cast(SpYrPeaks,Species~Winter,value = "PeakCount")
tableExport2<-cast(SpYrPeaks,Species~Winter,value = "Prop")

write.csv(tableExport,file="SpeciesPeakCounts.csv")
write.csv(tableExport2,file="SpeciesProps.csv")

postworks<-subset(SpYrPeaks, SpYrPeaks$Period=="Post-works")
write.csv(postworks,file="postworks.csv")

historic data summary table

historic<-subset(AllData, Winter<2007)
histdata<-cast(historic, Species+Winter~Month,sum,value="Number")
histdata2<-melt(histdata)
histdata<-cast(histdata2,Species+Winter~.,max,value="value")
histdata2<-cast(histdata,Species~Winter,value="(all)")

write.csv(histdata2,file=“HistoricDataSummary.csv”)

WeBS data

WeBS<-read.csv("Thames Estuary WeBS_2017.CSV")
WeBS<-WeBS[,c(1:4)]

View(WeBS)

WeBSTotal<-cast(WeBS,Period+Winter~.,sum,value="PeakCount")
WeBSTotal<-rename(WeBSTotal,c("(all)"="Total"))
WeBSTotal$Period<-factor(WeBSTotal$Period, levels = c("Pre-works","During","Post-works"))
WeBSKeySpp<-subset(WeBS,Species=="Avocet"|Species=="Black-tailed godwit"|Species=="Dunlin")

plot graphs

library(ggplot2)

# Total assemblage trend
fig1a<-ggplot(WinterPeakAssmb,aes(x=Winter,y=PeakAssmb))+
  #coord_cartesian(ylim=c(0,17500))+
  scale_x_continuous(limits=c(1998,2016),breaks=seq(1998,2016,by=2))+
  geom_smooth(method="lm",colour="gray15",aes(linetype=Period))+
  geom_point(colour="gray15")+
  labs(y="Winter peak assemblage")+
  scale_linetype_manual(values=c("solid","dashed","dotted"))+
  theme(legend.position="none")+
  #theme(legend.position="bottom")+
  theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1))+
  theme(legend.title=element_text(size=11,face="bold"),legend.key.width=unit(1.5,"cm"))+
  theme(axis.title=element_text(size=14), axis.text=element_text(size=10,colour="black"))+
  theme(panel.background=element_blank(),panel.grid.major = element_line(colour="gray90"),
        panel.grid.minor = element_blank(),panel.border = element_rect(colour = "gray",fill=NA),strip.background=element_rect(colour="gray"))
fig1a
## Warning in qt((1 - level)/2, df): NaNs produced

stat test: population trend post-works?

data<-subset(WinterPeakAssmb,Period=="Post-works")
attach(data)
model1<-lm(PeakAssmb~Winter,data=data)
summary(model1) # apparent increase in total assemblage since 2012/13 is not significant (p=0.712)
## 
## Call:
## lm(formula = PeakAssmb ~ Winter, data = data)
## 
## Residuals:
##       9      10      11      12      13 
##   506.4   240.5   557.6 -3862.3  2557.8 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -693818.2  1728011.7  -0.402    0.715
## Winter          347.9      858.0   0.405    0.712
## 
## Residual standard error: 2713 on 3 degrees of freedom
## Multiple R-squared:  0.05196,    Adjusted R-squared:  -0.2641 
## F-statistic: 0.1644 on 1 and 3 DF,  p-value: 0.7123
par(mfcol=c(2,2))
plot(model1) #model assumptions not seriously violated 

par(mfcol=c(1,1))
acf(residuals(model1)) #no problem with non-independence

detach(data)
fig1b<-ggplot(WeBSTotal,aes(x=Winter,y=Total))+
  geom_smooth(method="lm",colour="gray15",aes(linetype=Period))+
  geom_point(colour="gray15")+
  scale_x_continuous(limits=c(1998,2016),breaks=seq(1998,2016,by=2))+
  scale_y_continuous(breaks=seq(120000,220000,by=20000))+
  coord_cartesian(ylim=c(115000,222000))+
  scale_linetype_manual(values=c("solid","dashed","dotted"))+
  labs(y="Thames Estuary total assemblage")+
  theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1))+
  theme(legend.position="none",legend.title=element_text(size=11,face="bold"),legend.key.width=unit(1.5,"cm"))+
  theme(axis.title=element_text(size=14), axis.text=element_text(size=10,colour="black"))+
  theme(panel.background=element_blank(),panel.grid.major = element_line(colour="gray90"),
        panel.grid.minor = element_blank(),panel.border = element_rect(colour = "gray",fill=NA),strip.background=element_rect(colour="gray"))
fig1b
## Warning in qt((1 - level)/2, df): NaNs produced

stat test: WeBS population trend post-works?

data2<-subset(WeBSTotal,Period=="Post-works")
attach(data2)
model1b<-lm(Total~Winter)
summary(model1b) #(p=0.185)
## 
## Call:
## lm(formula = Total ~ Winter)
## 
## Residuals:
##      1      2      3      4 
## -13054  18734   1693  -7374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30586521   15316000   1.997    0.184
## Winter        -15118       7607  -1.987    0.185
## 
## Residual standard error: 17010 on 2 degrees of freedom
## Multiple R-squared:  0.6639, Adjusted R-squared:  0.4958 
## F-statistic:  3.95 on 1 and 2 DF,  p-value: 0.1852
par(mfcol=c(2,2))
plot(model1b) #model assumptions OK (except leverage) 

par(mfcol=c(1,1))
acf(residuals(model1b)) #no problem with non-independence

detach(data2)

Key species peak counts

fig2a<-ggplot(SpYrPeaks,aes(x=Winter, y=PeakCount))+
  facet_grid(~Species)+
  scale_y_continuous(breaks=seq(0,12000,by=2000))+
  coord_cartesian(ylim=c(0,12000))+
  scale_x_continuous(limits=c(1998,2016),breaks=seq(1998,2016,by=2))+
  labs(y="Winter peak count")+
  geom_smooth(method="lm",colour="gray15",aes(linetype=Period))+
  geom_point(colour="gray15")+
  scale_linetype_manual(values=c("solid","dashed","dotted"))+
  #theme(legend.position="bottom")+
  theme(legend.position= "none")+
  theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1))+
  theme(legend.text=element_text(size=18),legend.title=element_text(size=24,face="bold"),legend.key.width=unit(1.5,"cm"))+
  theme(axis.title=element_text(size=24), axis.text=element_text(size=18,colour="black"),strip.text=element_text(size=18,face="bold"))+
  theme(panel.background=element_blank(),panel.grid.major = element_line(colour="gray90"),
        panel.grid.minor = element_blank(),panel.border = element_rect(colour = "gray",fill=NA),strip.background=element_rect(colour="gray"))

fig2a
## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

stats: tests for population trends post-works

AVdata<-subset(SpYrPeaks,Period=="Post-works" & Species=="Avocet")
attach(AVdata)
model2<-lm(PeakCount~Winter)
summary(model2) # no significant trend in peak counts since 2012/13 (p=0.284)
## 
## Call:
## lm(formula = PeakCount ~ Winter)
## 
## Residuals:
##      1      2      3      4      5 
##  679.6 -447.5 -637.6 -100.7  506.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 555199.2   424857.9   1.307    0.282
## Winter        -274.9      211.0  -1.303    0.284
## 
## Residual standard error: 667.1 on 3 degrees of freedom
## Multiple R-squared:  0.3615, Adjusted R-squared:  0.1486 
## F-statistic: 1.698 on 1 and 3 DF,  p-value: 0.2835
par(mfcol=c(2,2))
plot(model2) #model assumptions OK 

par(mfcol=c(1,1))
acf(residuals(model2)) #no problem with non-independence

detach(AVdata)
BWdata<-subset(SpYrPeaks,Period=="Post-works" & Species=="Black-tailed godwit")
attach(BWdata)
model3<-lm(PeakCount~Winter)
summary(model3) # no significant change in peak counts since 2012/13 (p=0.24)
## 
## Call:
## lm(formula = PeakCount ~ Winter)
## 
## Residuals:
##       1       2       3       4       5 
##  -263.0   649.3   187.6 -1271.1   697.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  22986.6   596059.8   0.039    0.972
## Winter         -10.3      296.0  -0.035    0.974
## 
## Residual standard error: 935.9 on 3 degrees of freedom
## Multiple R-squared:  0.0004036,  Adjusted R-squared:  -0.3328 
## F-statistic: 0.001211 on 1 and 3 DF,  p-value: 0.9744
par(mfcol=c(2,2))
plot(model3) #model assumptions OK 

par(mfcol=c(1,1))
acf(residuals(model3)) #no problem with non-independence

detach(BWdata)
DNdata<-subset(SpYrPeaks,Period=="Post-works" & Species=="Dunlin")
attach(DNdata)
model4<-lm(PeakCount~Winter)
summary(model4) # no significant change in peak counts since 2012/13 (p=0.48)
## 
## Call:
## lm(formula = PeakCount ~ Winter)
## 
## Residuals:
##       1       2       3       4       5 
##   365.0  -912.2  2152.6 -3028.6  1423.2 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  446949.6  1507002.6   0.297    0.786
## Winter         -219.8      748.3  -0.294    0.788
## 
## Residual standard error: 2366 on 3 degrees of freedom
## Multiple R-squared:  0.02796,    Adjusted R-squared:  -0.2961 
## F-statistic: 0.08629 on 1 and 3 DF,  p-value: 0.7881
par(mfcol=c(2,2))
plot(model4) #model assumptions OK 

par(mfcol=c(1,1))
acf(residuals(model4)) #no problem with non-independence

detach(DNdata)

comparison between pre- and post-works counts

library(plyr)

compData<-subset(WinterPeakAssmb,Period!=“During”) compData<-subset(compData,Winter!=2007) compData<-subset(compData,Winter!=2008) comp<-ddply(compData, c(“Period”), summarise, N = sum(!is.na(PeakAssmb)), mean = mean(PeakAssmb,na.rm=TRUE), sd = sd(PeakAssmb,na.rm=TRUE), se = sd / sqrt(N))

fig3<-ggplot(comp,aes(x=Period,y=mean,fill=Period))+ geom_bar(stat=“identity”)+ geom_errorbar(aes(ymin=mean-se,ymax=mean+se),width=0.2, size=0.75, colour=“black”)+ labs(y=“Mean peak assemblage”)+ scale_fill_manual(values=c(“gray30”,“gray55”))+ theme(legend.position=“none”)+ theme(axis.title=element_text(size=14), axis.text=element_text(size=10,colour=“black”))+ theme(panel.background=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.border = element_rect(colour = “gray”,fill=NA),strip.background=element_rect(colour=“gray”))

stats: compare 4-yr mean peak assemblages

attach(compData) model5<-lm(PeakAssmb~Period) summary(model5) #no significant difference in 4-yr mean total assemblage (5-yr post-works) (p=0.56) par(mfcol=c(2,2)) plot(model5) #model assumptions OK par(mfcol=c(1,1)) acf(residuals(model5)) #no problem with non-independence detach(compData)

SpCompData<-subset(SpYrPeaks,Period!=“During”) SpCompData<-subset(SpCompData,Winter!=2007) SpCompData<-subset(SpCompData,Winter!=2008) SpComp<-ddply(SpCompData, c(“Species”,“Period”), summarise, N = sum(!is.na(PeakCount)), mean = mean(PeakCount,na.rm=TRUE), sd = sd(PeakCount,na.rm=TRUE), se = sd / sqrt(N))

fig4<-ggplot(SpComp,aes(Species,mean,fill=Period))+ geom_bar(stat=“identity”,position=“dodge”)+ geom_errorbar(aes(ymin=mean-se,ymax=mean+se), size=0.5, colour=“black”,position=position_dodge(0.9),width=0.2)+ labs(y=“Mean peak count”)+ scale_y_continuous(limits = c(0,8500))+ scale_fill_manual(values=c(“gray30”,“gray55”))+ theme(legend.position=“bottom”,legend.title=element_text(size=11,face=“bold”))+ theme(axis.title=element_text(size=14), axis.text=element_text(size=10,colour=“black”))+ theme(panel.background=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.border = element_rect(colour = “gray”,fill=NA),strip.background=element_rect(colour=“gray”))

stats compare key species’ 4-yr mean peaks and contributions to assemblage

SpYrPeaks<-subset(SpYrPeaks,Winter!=2007 & Winter!=2008) AVdata2<-subset(SpYrPeaks,Period!=“During” & Species==“Avocet”) attach(AVdata2) model6<-lm(PeakCount~Period) summary(model6) #significant increase in AV winter peak counts(p=0.012) par(mfcol=c(2,2)) plot(model6) #model assumptions OK par(mfcol=c(1,1)) acf(residuals(model6)) #no problem with non-independence

model7<-lm(Prop~Period) summary(model7) #significant increase in AV contribution to assemblage (p=0.018) par(mfcol=c(2,2)) plot(model7) #model assumptions OK par(mfcol=c(1,1)) acf(residuals(model7))#no problem with non-independence

detach(AVdata2)

BWdata2<-subset(SpYrPeaks,Period!=“During” & Species==“Black-tailed godwit”) attach(BWdata2) model8<-lm(PeakCount~Period) summary(model8) # significant increase in BW winter peak counts (p=0.002) par(mfcol=c(2,2)) plot(model8) # model assumptions OK par(mfcol=c(1,1)) acf(residuals(model8)) #no problem with non-independence

model9<-lm(Prop~Period) summary(model9) # significant increase in BW contribution to assemblage (p=0.0001) par(mfcol=c(2,2)) plot(model9) #model assumptions ok par(mfcol=c(1,1)) acf(residuals(model9))#no problem with non-independence

detach(BWdata2)

DNdata2<-subset(SpYrPeaks,Period!=“During” & Species==“Dunlin”) attach(DNdata2) model10<-lm(PeakCount~Period) summary(model10) # no significant change in DN winter peak counts (p=0.09) par(mfcol=c(2,2)) plot(model10) #model assumptions ok par(mfcol=c(1,1)) acf(residuals(model10)) #no problem with non-independence

model11<-lm(Prop~Period) summary(model11) # significant decrease in DN contribution to assemblage (p=0.03) par(mfcol=c(2,2)) plot(model11) #model assumptions not seriously violated par(mfcol=c(1,1)) acf(residuals(model11))#no problem with non-independence

detach(DNdata2)

figure: compare contributions to assemblage

SpComp2<-ddply(SpCompData, c(“Species”,“Period”), summarise, N = sum(!is.na(Prop)), mean = mean(Prop,na.rm=TRUE), sd = sd(Prop,na.rm=TRUE), se = sd / sqrt(N))

fig5<-ggplot(SpComp2,aes(Species,mean,fill=Period))+ geom_bar(stat=“identity”,position=“dodge”)+ geom_errorbar(aes(ymin=mean-se,ymax=mean+se), size=0.5, colour=“black”,position=position_dodge(0.9),width=0.2)+ labs(y=“Mean contribution to assemblage (%)”)+ scale_y_continuous(limits = c(0,100),breaks=seq(0,100,by=25))+ scale_fill_manual(values=c(“gray30”,“gray55”))+ theme(legend.position =“bottom”,legend.title=element_text(size=11,face=“bold”))+ theme(axis.title=element_text(size=14), axis.text=element_text(size=10,colour=“black”))+ theme(panel.background=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.border = element_rect(colour = “gray”,fill=NA),strip.background=element_rect(colour=“gray”))

export figures

ggsave(fig1a,file=“1a-Total assemblage trend.png”,width=140,height=100,units=“mm”,dpi=300) ggsave(fig1b,file=“1b-WeBS assemblage trend.png”,width=140,height=100,units=“mm”,dpi=300) ggsave(fig2a,file=“2a-Species winter peaks trends.png”,width=300,height=140,units=“mm”,dpi=300) ggsave(fig2b,file=“2b-WeBS species trends.png”,width=300,height=140,units=“mm”,dpi=300) ggsave(fig3,file=“3-Mean total assemblage comparison.png”,width=80,height=100,units=“mm”,dpi=300) ggsave(fig4,file=“4-Species mean peak comparision.png”,width=140,height=100,units=“mm”,dpi=300) ggsave(fig5,file=“5-Species mean contribution comparision.png”,width=140,height=120,units=“mm”,dpi=300)