rm(list=ls())

library(MASS)
library(lme4)
library(lmerTest)
library(languageR)
library(Hmisc)
library(multcomp)
library(car)
library(nlme)
library(car)
library(pgirmess)
library(lsmeans)
library(FactoMineR)
library(factoextra)
library(ade4)
library(psych)
library(effects)
library(lme4)
library(ggplot2)
library(arm)
library(sjPlot)

###########################################################################################################################
################################################ CHLORPYRIFOS #############################################################
###########################################################################################################################
top <- read.table("//smbhome.uscs.susx.ac.uk/ct430/Desktop/Toxic load paper/15 top pesticides.txt", 
                  header=TRUE, sep="\t") 
top

year=top$Years
chlor=top$chlorpyrifos

m <- lm(chlor ~ year)
summary(m)

plot(residuals(m))
qqnorm(residuals(m))
qqline(residuals(m))
hist(residuals(m))

plot(year, chlor, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "chlorpyrifos",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)
lines(year, chlor, col="dark red", lty=1, lwd=3)
abline(h=mean(chlor), col="darkgrey", lty=2, lwd=1)


###########################################################################################################################
################################################ Pendimethalin ###########################################################
###########################################################################################################################


pendi=top$pendimethalin


m1 <- lm(pendi ~ year)
summary(m1)

plot(residuals(m1))
qqnorm(residuals(m1))
qqline(residuals(m1))
hist(residuals(m1))

plot(year, pendi, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "pendimethalin",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, pendi, col="dark orange", lty=1, lwd=3)
abline(m1)

###########################################################################################################################
################################################ Ethoprop ###########################################################
###########################################################################################################################
etho=top$ethoprop


m2 <- lm(etho ~ year)
summary(m2)

plot(residuals(m2))
qqnorm(residuals(m2))
qqline(residuals(m2))
hist(residuals(m2))

plot(year, etho, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "ethoprop",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, etho, col="dark orange", lty=1, lwd=3)
abline(m2)

###########################################################################################################################
################################################ Chlormequat ###########################################################
###########################################################################################################################

chlorM=top$chlormequat


m3 <- lm(chlorM ~ year)
summary(m3)

plot(residuals(m3))
qqnorm(residuals(m3))
qqline(residuals(m3))
hist(residuals(m3))

plot(year, chlorM, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "chlormequat",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, chlorM, col="dark orange", lty=1, lwd=3)
abline(m3)

###########################################################################################################################
################################################ Mancozeb ###########################################################
###########################################################################################################################

manco=top$mancozeb


m4 <- lm(manco ~ year)
summary(m4)

plot(residuals(m4))
qqnorm(residuals(m4))
qqline(residuals(m4))
hist(residuals(m4))

plot(year, manco, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "mancozeb",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, manco, col="dark orange", lty=1, lwd=3)
abline(m4)

###########################################################################################################################
################################################ Oxamyl ###########################################################
###########################################################################################################################

oxa=top$oxamyl


m5 <- lm(oxa ~ year)
summary(m5)

plot(residuals(m5))
qqnorm(residuals(m5))
qqline(residuals(m5))
hist(residuals(m5))

plot(year, oxa, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "oxamyl",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, oxa, col="dark orange", lty=1, lwd=3)
abline(m5)

###########################################################################################################################
################################################ FOSTHIAZATE ##############################################################
###########################################################################################################################

fosth=top$fosthiazate


m6 <- lm(fosth ~ year)
summary(m6)

plot(residuals(m6))
qqnorm(residuals(m6))
qqline(residuals(m6))
hist(residuals(m6))

plot(year, fosth,
     ylab="potential corn bunting killed", 
     xlab = "year", 
     cex.lab=1.25,
     main = "fosthiazate",
     pch=".",
     col="darkslategray3",
     xaxt="none")
axis(1, seq(1998,2016,1), las=2, cex=0.25)
lines(year, fosth, col="darkslategray3", lty=1, lwd=3)
abline(m6)

###########################################################################################################################
################################################ GLYPHOSATE ###############################################################
###########################################################################################################################

glyph=top$glyphosate


m7 <- lm(glyph ~ year)
summary(m7)

plot(residuals(m7))
qqnorm(residuals(m7))
qqline(residuals(m7))
hist(residuals(m7))

plot(year, glyph,
     ylab="potential corn bunting killed", 
     xlab = "year", 
     cex.lab=1.25,
     main = "glyphosate",
     pch=".",
     col="darkorchid",
     xaxt="none")
axis(1, seq(1990,2016,1), las=2, cex=0.25)
lines(year, glyph, col="darkorchid", lty=1, lwd=3)
abline(m7)

###########################################################################################################################
################################################ Diquat ###########################################################
###########################################################################################################################

diq=top$diquat


m8 <- lm(diq ~ year)
summary(m8)

plot(residuals(m8))
qqnorm(residuals(m8))
qqline(residuals(m8))
hist(residuals(m8))

plot(year, diq, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "diquat",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, diq, col="dark orange", lty=1, lwd=3)
abline(m8)

###########################################################################################################################
################################################ MCPA ###########################################################
###########################################################################################################################

mcpa=top$MCPA


m9 <- lm(mcpa ~ year)
summary(m9)

plot(residuals(m9))
qqnorm(residuals(m9))
qqline(residuals(m9))
hist(residuals(m9))

plot(year, mcpa, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "MCPA",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, mcpa, col="dark orange", lty=1, lwd=3)
abline(m9)

###########################################################################################################################
################################################ Chlorothalonil ###########################################################
###########################################################################################################################

chlorT=top$chlorothalonil


m10 <- lm(chlorT ~ year)
summary(m10)

plot(residuals(m10))
qqnorm(residuals(m10))
qqline(residuals(m10))
hist(residuals(m10))

plot(year, chlorT, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "chlorothalonil",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, chlorT, col="dark orange", lty=1, lwd=3)
abline(m10)

###########################################################################################################################
################################################ Mecoprop-P ###########################################################
###########################################################################################################################

meco=top$mecoprop.p

m11 <- lm(meco ~ year)
summary(m11)

plot(residuals(m11))
qqnorm(residuals(m11))
qqline(residuals(m11))
hist(residuals(m11))

plot(year, meco, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "mecoprop-P",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, meco, col="dark orange", lty=1, lwd=3)
abline(m11)

###########################################################################################################################
################################################ Cyproconazole ###########################################################
###########################################################################################################################

cypro=top$cyproconazole


m12 <- lm(cypro ~ year)
summary(m12)

plot(residuals(m12))
qqnorm(residuals(m12))
qqline(residuals(m12))
hist(residuals(m12))

plot(year, cypro, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "cyproconazole",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, cypro, col="dark orange", lty=1, lwd=3)
abline(m12)

###########################################################################################################################
################################################ Metaldehyde ###########################################################
###########################################################################################################################

metal=top$metaldehyde

m13 <- lm(metal ~ year)
summary(m13)

plot(residuals(m13))
qqnorm(residuals(m13))
qqline(residuals(m13))
hist(residuals(m13))

plot(year, metal, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "metaldehyde",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, metal, col="dark orange", lty=1, lwd=3)
abline(m13)

###########################################################################################################################
################################################ Cyprodinil ###########################################################
###########################################################################################################################


cyproD=top$cyprodinil

m14 <- lm(cyproD ~ year)
summary(m14)

plot(residuals(m14))
qqnorm(residuals(m14))
qqline(residuals(m14))
hist(residuals(m14))

plot(year, cyproD, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "cyprodinil",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, cyproD, col="dark orange", lty=1, lwd=3)
abline(m14)

###########################################################################################################################
############################################## TOTAL PESTICIDES ###########################################################
###########################################################################################################################
pesticide<- read.table("//smbhome.uscs.susx.ac.uk/ct430/Desktop/Toxic load paper/pesticides.txt", 
                       header=TRUE, sep="\t") 
pesticide

year=pesticide$years
herbicide=pesticide$HERBICIDES
fungicide=pesticide$FUNGICIDES
other=pesticide$OTHER
insecticide=pesticide$INSECTICIDES
total=pesticide$TOTAL

par(mfrow=c(2,2))

############# INSECTICIDES #############
m15 <- lm (insecticide ~ year, data = pesticide)
summary(m15)

plot(residuals(m15))
qqnorm(residuals(m15))
qqline(residuals(m15))
hist(residuals(m15))

plot(year, insecticide, type = "o",
     ylim = c(8e+11,1.65e+13),
     ylab="", 
     xlab = "year", pch=".",
     cex.lab=1.30,
     cex.main=1.25,
     col="orange",
     lty=1,
     lwd=2,
     main = "insecticides",
     xaxt="none", yaxt="none")
axis(1, seq(1990,2015,5), las=2, cex=0.15)
axis(2,seq(8e+11,1.65e+13,1.57e+13), cex=0.15)
abline(m15)

############# HERBICIDES #############
m16 <- lm(herbicide ~ year, data=pesticide)
summary(m16)

plot(residuals(m16))
qqnorm(residuals(m16))
qqline(residuals(m16))
hist(residuals(m16))

plot(year, herbicide,type = "o", 
                    ylim = c(4.5e+11,1.1e+12),
                   ylab="", 
                   xlab = "year",
                   cex.lab=1.30,
                   pch=".",
                   main = "herbicides",
                   cex.main=1.25,
                   col="red",
                   lty=1,
                   lwd=2, yaxt="none", xaxt="none")
axis(1, seq(1990,2015,5), las=2, cex=0.15)
axis(2, seq(4.5e+11,1.1e+12,6.5e+11), cex=0.15)
abline(m16)

############# FUNGICIDES #############
m17 <- lm(fungicide ~ year, data=pesticide)
summary(m17)

plot(residuals(m17))
qqnorm(residuals(m17))
qqline(residuals(m17))
hist(residuals(m17))

plot(year, fungicide,type = "o", 
                    ylim = c(2.5e+11,5.3e+11),
                    ylab="", 
                    xlab = "year",
                    cex.lab=1.30,
                    cex.main=1.25,
                    pch=".",
                    main = "fungicides", 
                    col="blue",
                    lty=1,
                    lwd=2, xaxt="none", yaxt="none")
axis(1, seq(1990,2015,5), las=2, cex=0.15)
axis(2,seq(2.5e+11,5.3e+11,2.8e+11), cex=0.15)
abline(m17)

############# OTHER PESTICIDES #############
m18 <- lm(other ~ year, data=pesticide)
summary(m18)

plot(residuals(m18))
qqnorm(residuals(m18))
qqline(residuals(m18))
hist(residuals(m18))

plot(year, other,type = "o", 
                    ylim = c(1.5e+11,5.7e+11),
                    ylab="", 
                    xlab = "year",
                    cex.lab=1.30,
                    cex.main=1.25,
                    pch=".",
                    main = "other pesticides",
                    col="green",
                    lty=1,
                    lwd=2, xaxt="none", yaxt="none")
axis(1, seq(1990,2015,5), las=2, cex=0.15)
axis(2,seq(1.5e+11,5.7e+11, 4.2e+11), cex=0.15)
abline(m18)


par(mfrow=c(1,1))
############# TOTAL ALL #############
m19 <- lm(total ~ year, data=pesticide)
summary(m19)

plot(residuals(m19))
qqnorm(residuals(m19))
qqline(residuals(m19))
hist(residuals(m19))

plot(year, total, ylab="potential corn bunting killed", 
     xlab = "year", 
     pch=".", 
     main= "All pesticides",
     xaxt="none",
     cex.lab=1.25)
axis(1, seq(1990,2016,1), las=2, cex=0.25)

lines(year, total, col="palevioletred3", lty=1, lwd=3)
abline(m19)


##########################################################################################################################
############################################## WEIGHT AND AREA ###########################################################
##########################################################################################################################
weight.area<- read.table("//smbhome.uscs.susx.ac.uk/ct430/Desktop/Toxic load paper/weigh and area applied.txt", 
                         header=TRUE, sep="\t") 
weight.area

year=weight.area$Year
weight=weight.area$Weight.applied..kg.
area=weight.area$Area.treated..ha.

m20 <- lm (weight ~ year, data = weight.area)
summary(m20)

plot(residuals(m20))
qqnorm(residuals(m20))
qqline(residuals(m20))
hist(residuals(m20))

m21 <- lm (area ~ year, data = weight.area)
summary(m21)

plot(residuals(m21))
qqnorm(residuals(m21))
qqline(residuals(m21))
hist(residuals(m21))


par(mfrow=c(1,1))

############# WEIGHT #############



plot(year, weight,type = "o",
                ylim = c(16321668, 80274553),
                    ylab="", 
                    xlab = "year", 
                cex.lab=1.25,
                    pch=".",
                    col="grey40",
                    lty=1,
                    lwd=3, xaxt="none")
axis(1, seq(1990,2016,1), las=2, cex=0.25)
abline(m20, col="black", lty=2, lwd=1)

############# AREA #############
points(year, area, 
       pch=".",
       col="grey80")
lines(year, area, col="grey80", lty=1, lwd=3)
abline(m21, col="grey70", lty=2, lwd=1)


# legend #
legend(1990,8.1e+07, legend = c("weight applied (kg)", "area treated (ha)"),
       col = c("grey40", "grey80"),
       lty=c(1,1),
       lwd = c(3,3),
       ncol=1, cex = 0.75)


########################################################################################################################################