# install.packages("plotrix")
# install.packages("pbkrtest")
# install.packages("afex")

library(plotrix)
library(sciplot)
library(ggplot2)
library(boot)
library(mgcv)
library(lme4)
library(MuMIn)
# library(pbkrtest)
# library(afex)

d <- read.csv("data_parasites.csv", fileEncoding="latin1")
names(d)
names(d)[8:10] <- c("Peso", "Talla", "P_bull")



pdf("Rplots_def.pdf")


def.par <- par(no.readonly=TRUE)



# Generar columna con 1/0
d$NI <- NA
d$NI <- ifelse(d$P_bull > 0, 1, 0)



my.aggre <- function(x){
	
	# Estimate prevalence (Ninfected/Ntotal)
	d1 <- aggregate(x$NI, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), function(i) sum(i)/length(i))
	names(d1)[4] <- "prev"
	
	# Aggregate bird density by LAT and sampling time
	d0 <- aggregate(x[, 17:21], list(LAT=x$LAT, SAMPL=x$season), mean)	
	d001 <- merge(d0, d1)
	
# Estimate parasite intensity (i.e. mean number of parasites in each parasited host)
	d01 <- droplevels(x[x$NI==1,  ])
		
	d12 <- aggregate(d01$P_bull, list(LAT=d01$LAT, SITE=d01$Lugar, SAMPL=d01$season), function(i) sum(i)/length(i))
	names(d12)[4] <- "intens"
		
	d002 <- merge(d001, d12)
		
	d13 <- aggregate(x$Talla, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), function(i) length(i))
	
	names(d13)[4] <- "N"
	d003 <- merge(d002, d13)
	
	d14 <- aggregate(d01$Talla, list(LAT=d01$LAT, SITE=d01$Lugar, SAMPL=d01$season), function(i) length(i))
	
	names(d14)[4] <- "Ninf"
	d15 <- merge(d003, d14)
	
# Estimate mean parasite load
	d16 <- aggregate(x$P_bull, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), mean)	
	names(d16)[4] <- "m"
	d17 <- merge(d15, d16)  

# Estimate parasite load var	
	d18 <- aggregate(x$P_bull, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season),  function(i) sd(i, na.rm=TRUE)^2 )	
	names(d18)[4] <- "s2"
	
	d19 <- merge(d17, d18)

# Corrected moment estimate of k from binomial negative distribution
# (s2 = m + m2/k)

my.k <- function(y) {
	m <- mean(y)
	s2 <- var(y)
	n <- length(y)
return( ( m^2 - (s2/n) ) / (s2-m) )
	
}

my.vmr <- function(y){
	m <- mean(y)
	s2 <- var(y)
return(s2/m)
}

my.I <- function(y){
	m <- mean(y)
	s2 <- var(y)
	n <- length(y)
return( (s2*(n-1))/m )
}
	d20 <- aggregate(x$P_bull, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), function(i) my.k(i))
names(d20)[4] <- "k"
	
	d21 <- merge(d19, d20)
		

	d22 <- aggregate(x$P_bull, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), function(i) my.vmr(i))
names(d22)[4] <- "vmr"

	d23 <- merge(d21, d22)

	d24 <- aggregate(x$P_bull, list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), function(i) my.I(i))
names(d24)[4] <- "I"
	
	d25 <- merge(d23, d24)
	d26 <- aggregate(x[, 24:27], list(LAT=x$LAT, SITE=x$Lugar, SAMPL=x$season), unique)

names(d26)[4:7] <- c("peop", "dogs", "cow", "total")

	return(merge(d25, d26)  )
	
}


# Mesoscale variation of temporal aggregation of parasite load

# All
d1 <- my.aggre(x = d)


# Split dataset by size modes ( > 15 mm < 15 mm)
d$SIZE <- NA
d$SIZE <- ifelse(d$Talla<15, 1, 2)
d$SIZE <- as.factor(d$SIZE)

par(mfrow=c(2, 2), las=1)
hist(d$Talla, main="Host size structure", xlab="Body size (mm)", ylab="Counts", bty="l")
plot(d$P_bull ~ d$Talla, xlab="Body size (mm)", ylab=expression("Nº of parasites ind"^-1), bty="l", main="Parasite burden")

atab <- .95
lac <- -3.6
lbd <- -25
cx <- 1.2
mtext(expression(bold("A")), 2, line=lac, cex=cx, outer=TRUE, las=1, at=atab)
mtext(expression(bold("B")), 2, line=lbd, cex=cx, outer=TRUE, las=1, at=atab)

par(def.par)



d2 <- my.aggre(droplevels(d[d$SIZE==1,  ]))
d2$SIZE <- as.factor(rep(1, nrow(d2)) ) 
d3 <- my.aggre(droplevels(d[d$SIZE==2,  ]))  
d3$SIZE <- as.factor(rep(2, nrow(d3)))

d4 <- rbind(d2, d3)


my.taylor <- function(dat, main=NULL, legend=TRUE)	{
	mods <- by(dat, -dat$LAT, function(i) lm(log(s2) ~ log(m), data=i))

print(lapply(mods, summary))

b <- rev(unlist(lapply(mods, function(i) coef(i)[2] )))
a <- rev(unlist(lapply(mods, function(i) coef(i)[1] )))

my.reg <- function(dat, id) {
    m1 <- lm(log10(s2) ~ log10(m), subset=id, data=dat)
    coef(m1)[1] + coef(m1)[2]
}

nR <- 999
l.m <- by(dat, -dat$LAT, function(i) boot(i, statistic=my.reg, R=nR))

b.m <- unlist(lapply(lapply(l.m, "[[", "t"), mean))
b.ci <- lapply(l.m, function(i) boot.ci(i, type="bca")$bca[4:5]  )
d.b <- do.call(rbind.data.frame, b.ci)
d.b$b.m <- b.m
names(d.b)[1:2] <- c("lci", "uci")

mm <- tapply(dat$densave, -dat$LAT, mean)
d.b$mm <- mm

bg <- c("red", "green4", "blue", "purple")[as.factor(-dat$LAT)]
s1 <- c("Cheuque", "Curiñanco", "Calfuco", "Chaihuín")
pcs <- c(16, 21, 15, 22)
plotCI(d.b$mm, d.b$b.m, ui=d.b$uci, li=d.b$lci, bty="l", xlab=expression("Bird abundance"), ylab= expression(italic("I")[10]*" "%+-%" CI"), pch=pcs, ylim=c(-1, 4), main=main, las=1, cex=1.2)
abline(h=0, lty=2)
abline(h=1, lty=3)
if(legend==TRUE) {legend("topright", legend=s1, pch=pcs, bty="n", lty=rep(1, 4))}
}



my.t.plot.1 <- function(Y, dat, legend=FALSE, ylab=NULL, main=NULL, ylim=NULL, ncol=2){
	d1 <- dat
	s.name <- c("Cheuque", "Curiñanco", "Calfuco", "Chaihuín")
	s.time <- c("2014", "2014", "2014", "2014", "2015", "2015", "2015", "2015")
	s.time2 <- c("Su", "Au", "Wi", "Sp", "Su", "Au", "Wi", "Sp")
	lineplot.CI(d1$SAMPL, Y, -d1$LAT, fixed=TRUE, x.leg="topright", 	leg.lab=s.name, bty="l", bg="grey", xaxt="n", xlab="", ylab=ylab, legend=legend, las=1, main=main, ylim=ylim, ncol=ncol)
	axis(1, at=1:8, labels=s.time2, cex.axis=.9)
	axis(1, at=1:8, labels=s.time, line=1.5, lwd=0)
	}
	
my.case4 <- function(four=TRUE, cex=cx) {
	atab <- .95
	atcd <- .45
	lac <- -3.5
	lbd <- -23
mtext(expression(bold("A")), 2, line=lac, cex=cx, outer=TRUE, las=1, at=atab)
mtext(expression(bold("B")), 2, line=lbd, cex=cx, outer=TRUE, las=1, at=atab)
mtext(expression(bold("C")), 2, line=lac, cex=cx, outer=TRUE, las=1, at=atcd)
if(four==TRUE) mtext(expression(bold("D")), 2, line=lbd, cex=cex, outer=TRUE, las=1, at=atcd)
	}
	
my.case2 <- function(cex=cx) {
	atab <- .95
	atcd <- .45
	lac <- -2.5
	lbd <- -23
mtext(expression(bold("A")), 2, line=lac, cex=cx, outer=TRUE, las=1, at=atab)
mtext(expression(bold("B")), 2, line=lac, cex=cx, outer=TRUE, las=1, at=atcd)
	}
	
# -----------------------------	
# Models
# -----------------------------
# Prev
# All bird summed. Include size as factor

# d$sDens <- scale(d$densave)
# d$sSize <- scale(d$Talla)
# ctrl <- glmerControl(optimizer="bobyqa")
# m01 <- glmer(NI ~  sDens*sSize + (1 + sDens | Lugar), data=d, family=binomial)

# tt <- getME(m01,"theta")
# ll <- getME(m01,"lower")
# coef(m01)
# min(tt[ll==0])
# plot(m01)
# summary(m01)
# predict(m01)
# plot(predict(m01) ~ sDens, data=d)
# d$pp <- predict(m01)
# a <- exp(0.5141010)
# b1 <- exp(1.5702218)
# ggplot(data=d, aes(x=sDens, y=pp, color=sSize))  + geom_point() + scale_color_gradient(low = "green", high = "red") +  geom_abline(intercept = 0.5141010, slope = 1.5702218, color = "green")


m1 <- lme(prev ~ densave*SIZE, random=list(SITE = pdDiag(~SAMPL)), correlation=corAR1(form= ~1| SITE), data=d4)

summary(m1)



# Obtaining approximate F-tests on the basis of Kenward-Roger approach
# Full
m1 <- lmer(prev ~ densave + SIZE + densave:SIZE + (1|SITE/SAMPL), data=d4 )

print(plot(m1, main="Prevalence vs. bird*size"))
acf(residuals(m1, type="deviance"), main="Prevalence vs. bird*size \nFull model")


print(summary(m1))
print(r.squaredGLMM(m1))


write.csv(summary(m1)[10], "res_LMM_prev.csv")


Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

Mode(d$Talla[d$Talla<20]  )
Mode(d$Talla[d$Talla>20]  )

mean(d3$prev)
sd(d3$prev)


diff(tapply(d4$prev, d4$SIZE, mean))

cof <- fixef(m1)


# Bird density

par(mfcol=c(2, 2), mgp=c(2.85, 1, 0))

my.t.plot.1(d2$densave, d2, ylab=expression("Bird abundance (ind km"^-1*")"), legend=TRUE, ylim=c(0, 1200), main="Bird abundance (parasite exposure)", ncol=2)

# Molecrab density
my.t.plot.1(d1$N, d1, ylab=expression("Molecrab density (ind 0.006 m"^-3*")"), legend=FALSE, ylim=c(0, 200), main="Molecrab abundance", ncol=2)


my.case2()

par(def.par)


min(tapply(d1$N, list(d1$SITE, d1$SAMPL), unique))



par(mfrow=c(2, 2), mgp=c(2.85, 1, 0))
my.t.plot.1(d2$prev, d2, ylab=expression("Prevalence"), main="Small hosts (<15 mm)", ylim=c(0, 1), legend=TRUE)

my.t.plot.1(d3$prev, d3, ylab=expression("Prevalence"), main="Large hosts (>15 mm)", ylim=c(0, 1))

gr <- rev(gray(1:2/2))
plot(d4$prev ~ d4$densave, bty="l", las=1, xlab=expression("Bird abundance (ind km"^-1*")"), ylab="Prevalence", pch=21, bg=gr[d4$SIZE], ylim=c(0, 1))
legend("bottomright", legend=c("Large", "Small"), pch=21, pt.bg=rev(gr[unique(d4$SIZE)]), cex=1.1, ncol=1, lty=c(2, 1), bty="n")

abline(a = cof[1], b = cof[2], lwd=1, lty=1)
abline(a = cof[1] + cof[3], b = cof[2] + cof[4], lwd=1, lty=2)


my.case4(four=FALSE)

par(def.par)





# Aggregation

# m2 <- lme(log10(vmr) ~ densave*SIZE, random=list(SITE = pdDiag(~SAMPL)), correlation=corAR1(form= ~1|SITE), data=d4)


m02 <- lmer(vmr ~ densave + SIZE + densave:SIZE + (1|SITE/SAMPL), data=d4 )

print(plot(m02, main="Variance-to-mean ratio vs. bird*size"))

m2 <- lmer(log10(vmr) ~ densave + SIZE + densave:SIZE + (1|SITE/SAMPL), data=d4 )

print(plot(m2, main="log10 variance-to-mean ratio vs. bird*size"))

acf(residuals(m2, type="deviance"), main="log10 variance-to-mean ratio vs. bird*size")


print(summary(m2))
print(r.squaredGLMM(m2))


write.csv(summary(m2)[10], "res_LMM_mvr.csv")


print(summary(m2))
print(r.squaredGLMM(m2))



write.table(c(capture.output(r.squaredGLMM(m1)), capture.output(r.squaredGLMM(m2))), "res_R2.txt")



par(mfrow=c(2, 2), mgp=c(2.85, 1, 0))

my.t.plot.1(d2$vmr, d2, ylab=expression("Variance-mean ratio (s"^2*"m"^-1*")"), main="Small hosts (<15 mm)", legend=TRUE, ylim=c(0, 7))
abline(h=1, lty=3)

my.t.plot.1(d3$vmr, d3, ylab=expression("Variance-mean ratio (s"^2*"m"^-1*")"), main="Large hosts (>15 mm)", ylim=c(0, 7))
abline(h=1, lty=3)

# plot(d4$vmr ~ d4$densave, bty="l", las=1, xlab=expression("Bird density (mean ind km"^-2*")"), ylab=expression("Variance-mean ratio (s"^2*"m"^-1*")"), pch=21, bg=gr[d4$SIZE])
# legend("topright", legend=c("Small", "Large"), pch=21, pt.bg=gr[unique(d4$SIZE)], cex=1.1, ncol=2)

my.taylor(dat=d2, legend=FALSE)
my.taylor(dat=d3, legend=FALSE)


my.case4()

par(def.par)


# Appendix

# How does relate bird count with human-related disturbance
plot(densave ~ total, data=d1, las=1, pch=21, bg="pink", bty="l", xlab=expression("Counts (ind site"^-1*"survey"^-1* ")"), ylab=expression("Bird abundance (ind km"^-1*")" ), cex=1.3)

m0 <- lme(log(densave+1) ~ total, random=list(SITE = pdDiag(~SAMPL)), correlation=corAR1(form= ~1| SITE), data=d1)


summary(m0)
print(r.squaredGLMM(m0))



# Lagged responses 


# N vs. bird density

par(mfcol=c(4, 4), las=1, mar=c(4, 4, 3, 0)+.1 )

by(d2, -d2$LAT, function(i) ccf(i$N, i$densave, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="") )
mtext("Small", outer=TRUE, line=-2.7, at=0.17)

by(d3, -d3$LAT, function(i) ccf(i$N, i$densave, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )
mtext("Large", outer=TRUE, line=-2.7, at=0.41)

par(las=1, mar=c(4, 0, 2, 1) )
plot.new()
legend("topleft", legend="Cheuque", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Curiñanco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Calfuco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Chaihuín", bty="n", cex=1.5)

mtext(expression(bold("Molecrab count and bird density")), 3, outer=TRUE, line=-1.5, at=0.3)

mtext("CCF", 2, outer=TRUE, line=-2, las=0)

mtext("Lag", 1, outer=TRUE, line=-1.2, at=0.3)


par(def.par)





# Bird density vs. prevalence

par(mfcol=c(4, 4), las=1, mar=c(4, 4, 2, 1) )

by(d2, -d2$LAT, function(i) ccf(i$densave, i$prev, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")) 

mtext("Small", outer=TRUE, line=-2.7, at=0.17)

by(d3, -d3$LAT, function(i) ccf(i$densave, i$prev, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )

mtext("Large", outer=TRUE, line=-2.7, at=0.41)

par(las=1, mar=c(4, 0, 2, 1) )
plot.new()
legend("topleft", legend="Cheuque", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Curiñanco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Calfuco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Chaihuín", bty="n", cex=1.5)


mtext(expression(bold("Bird density and prevalence")), 3, outer=TRUE, line=-1.5, at=0.3)

mtext("CCF", 2, outer=TRUE, line=-2, las=0)

mtext("Lag", 1, outer=TRUE, line=-1.2, at=0.3)



par(def.par)



# Bird density vs. aggregation
par(mfcol=c(4, 4), las=1, mar=c(4, 4, 2, 1) )
by(d2, -d2$LAT, function(i) ccf(i$densave, i$vmr, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )

mtext("Small", outer=TRUE, line=-2.7, at=0.17)

by(d3, -d3$LAT, function(i) ccf(i$densave, i$vmr, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )

mtext("Large", outer=TRUE, line=-2.7, at=0.41)

par(las=1, mar=c(4, 0, 2, 1) )
plot.new()
legend("topleft", legend="Cheuque", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Curiñanco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Calfuco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Chaihuín", bty="n", cex=1.5)

mtext(expression(bold("Bird density and VMR")), 3, outer=TRUE, line=-1.5, at=0.3)

mtext("CCF", 2, outer=TRUE, line=-2, las=0)

mtext("Lag", 1, outer=TRUE, line=-1.2, at=0.3)


par(def.par)



# N vs. prevalence

par(mfcol=c(4, 4), las=1, mar=c(4, 4, 2, 1) )
by(d2, -d2$LAT, function(i) ccf(i$N, i$prev, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="") )

mtext("Small", outer=TRUE, line=-2.7, at=0.17)

by(d3, -d3$LAT, function(i) ccf(i$N, i$prev, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )

mtext("Large", outer=TRUE, line=-2.7, at=0.41)

par(las=1, mar=c(4, 0, 2, 1) )
plot.new()
legend("topleft", legend="Cheuque", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Curiñanco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Calfuco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Chaihuín", bty="n", cex=1.5)

mtext(expression(bold("Molecrab count and prevalence")), 3, outer=TRUE, line=-1.5, at=0.3)

mtext("CCF", 2, outer=TRUE, line=-2, las=0)

mtext("Lag", 1, outer=TRUE, line=-1.2, at=0.3)


par(def.par)


# N vs. aggregation

par(mfcol=c(4, 4), las=1, mar=c(4, 4, 2, 1) )

by(d2, -d2$LAT, function(i) ccf(i$N, i$vmr, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="") )

mtext("Small", outer=TRUE, line=-2.7, at=0.17)

by(d3, -d3$LAT, function(i) ccf(i$N, i$vmr, na.action=na.pass, xlim=c(0, 6), ylim=c(-1, 1), bty="l", main="", xlab="", ylab="")   )

mtext("Large", outer=TRUE, line=-2.7, at=0.41)

par(las=1, mar=c(4, 0, 2, 1) )
plot.new()
legend("topleft", legend="Cheuque", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Curiñanco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Calfuco", bty="n", cex=1.5)
plot.new()
legend("topleft", legend="Chaihuín", bty="n", cex=1.5)

mtext(expression(bold("Molecrab count and aggregation")), 3, outer=TRUE, line=-1.5, at=0.3)

mtext("CCF", 2, outer=TRUE, line=-2, las=0)

mtext("Lag", 1, outer=TRUE, line=-1.2, at=0.3)


par(def.par)

names(d1)

plot(s2 ~ m, data=d3)

range(d3$m)

d3$ran <- ifelse(d3$m < (median(d3$m)), "A", "B")


library(lattice)


xyplot(s2 ~ m | ran, data=d3, xlab="Mean")

mm1 <- by(d3, d3$ran, function(i) lm(log10(s2) ~ log10(m), data=i)  )

lapply(mm1, summary)

dev.off()


rm(list=ls(all=TRUE))