## packages 
library(plyr)
library(ggplot2)
# run primary code up to Ovitraps section#####
mypath    <- "C:/Users/Mattia/Documents/LaSapienza/HLC/Submission/Data/"

dbovi <- read.csv(paste(mypath,"2014_ovitrap_Orto vs Antomia.csv",sep=""),sep=";") 
summary(dbovi)
levels(dbovi$DaysInField)[4] <- "4gg"

dbovi$DaysInFieldN <- as.numeric(dbovi$DaysInField)+1
dbovi$leftover     <- dbovi$EndH2O
dbovi$Evapday      <- (dbovi$StartH20 - dbovi$EndH2O)
table(dbovi$OvitrapID,dbovi$Site)



# 
dbovi$julian <- as.numeric(strftime(strptime(dbovi$Date,format="%d/%m/%Y"),origin="01/01/14",format="%j") )
dbovi$week   <- as.numeric(strftime(strptime(dbovi$Date,format="%d/%m/%Y"),origin="01/01/14",format="%W") )

#plot(dbovi$julian,dbovi$NrEggs)
nsim = 1000
corval <- array(NA,c(nsim,2,15))
modval <- array(NA,c(nsim,3,15))
eggval <- array(NA,c(nsim,73,15))
for(o in 1:15){
  
  for(s in 1:nsim){

idchose <- droplevels(sample(unique(dbovi$OvitrapID),o,replace = TRUE))

levels(dbovi$Site) <- c("Site-B","Site-A")

db1 <- droplevels(subset(dbovi, is.na(NrEggs)==FALSE & OvitrapID %in% idchose))

#db1$Site         <- relevel(db1$Site,ref="Site-A")
dbhlc <- ddply(dbhlc1,c("Site","month","Date","julian"),summarize,
               totfem  = sum(Female_albopictus),
               meanfem = mean(Female_albopictus),
               length  = length(Female_albopictus),
               std     = sd(Female_albopictus),
               seup    =( meanfem + 1.96*std/sqrt(length)),
               selo    =( meanfem - 1.96*std/sqrt(length)),
               wind    = mean(wind),
               mmrn    = mean(mmrn),
               mmrn7   = mean(mmrn7),    
               mmrn14  = mean(mmrn14),
               mmrn21  = mean(mmrn21),
               tm2m    = mean(tm2m),
               tm2m7   = mean(tm2m7),
               tm2m14  = mean(tm2m14),
               tm2m21  = mean(tm2m21)
               )  

range(dbhlc$julian)


#db1$NrEggs1 <- log10(db1$NrEggs+1)
ovitrapp <- ddply(db1,c("Site","julian","Date"),summarize,
                  uova     = mean(NrEggs,na.rm=TRUE),
                  std      = sd(NrEggs),
                  length   = length(NrEggs),
                  seup     =( uova + 1.96*std/sqrt(length)),
                  selo     =( uova - 1.96*std/sqrt(length)),
                  totuova  = sum(NrEggs,na.rm=TRUE),
                  tempo    = mean(DaysInFieldN),
                  Evapday  = mean(Evapday,na.rm=TRUE),
                  leftover = mean(leftover,na.rm=TRUE),
                  ntrap    = sum(is.na(NrEggs)==FALSE))

day <- 189:304




# Eggs lag ###
if(sum(db1$Site=="Site-A")>0){
  ovitrappA <- (subset(ovitrapp, Site =="Site-A"))

eg1A <- rep(ovitrappA$uova[1],each=ovitrappA$tempo[1])
d1A  <- c((ovitrappA$julian[1]-ovitrappA$tempo[1]+1):ovitrappA$julian[1])
for(i in 2:length(ovitrappA$uova)){
  eg0A <- rep(ovitrappA$uova[i]/ovitrappA$tempo[i],each=ovitrappA$tempo[i])
  d0A  <- c((ovitrappA$julian[i]-ovitrappA$tempo[i]+1):ovitrappA$julian[i])
  eg1A <- c(eg1A,eg0A)
  d1A  <- c(d1A,d0A)}


mdayA <- rep(NA,length(day))
for(i in 1:length(day)){
  if(sum(d1A==day[i])==1){
    mdayA[i] <- eg1A[d1A==day[i]]}
  else( mdayA[i] <-NA)
}
dayA <- ifelse(is.na(mdayA),0,1)
}else{
  ovitrappA <- data.frame(julian = -1)
}

if(sum(db1$Site=="Site-B")>0){

ovitrappB <- (subset(ovitrapp, Site =="Site-B"))

eg1B <- rep(ovitrappB$uova[1],each=ovitrappB$tempo[1])
d1B  <- c((ovitrappB$julian[1]-ovitrappB$tempo[1]+1):ovitrappB$julian[1])
for(i in 2:length(ovitrappB$uova)){
  eg0B <- rep(ovitrappB$uova[i]/ovitrappB$tempo[i],each=ovitrappB$tempo[i])
  d0B  <- c((ovitrappB$julian[i]-ovitrappB$tempo[i]+1):ovitrappB$julian[i])
  eg1B <- c(eg1B,eg0B)
  d1B  <- c(d1B,d0B)}

mdayB <- rep(NA,length(day))
for(i in 1:length(day)){
  if(sum(d1B==day[i])==1){
    mdayB[i] <- eg1B[d1B==day[i]]}
  else( mdayB[i] <-NA)
}
dayB <- ifelse(is.na(mdayB),0,1)
}else{
  ovitrappB <- data.frame(julian = -1)
}


dbhlc$E.day    <- rep(NA,nrow(dbhlc))
for(i in 1:nrow(dbhlc)){
  dd <- day %in% c(dbhlc$julian[i]:(dbhlc$julian[i]-6))
  dd2 <- day %in% c((dbhlc$julian[i]-7):(dbhlc$julian[i]-13))
  if(dbhlc$Site[i] == "Site-A"){
    if(dbhlc$julian[i] %in% ovitrappA$julian){dbhlc$E.day[i] <- mdayA[day == dbhlc$julian[i]]}
    else{dbhlc$E.day[i] <-NA}
  }else{
    if(dbhlc$julian[i] %in% ovitrappB$julian){dbhlc$E.day[i] <- mdayB[day == dbhlc$julian[i]]}
    else{dbhlc$E.day[i] <-NA}    
   }
}



# merge HLC, Egg and climate dataset

library(nlme)
# Model II ######
db2        <- subset(dbhlc, is.na(E.day)==FALSE )

# correlation ####
corval[ s , 1 , o] <- cor.test(db2$meanfem,db2$E.day)$estimate
corval[ s , 2 , o] <- cor.test(db2$meanfem,db2$E.day)$p.value

#eggval[ s ,1:length(db2$E.day)  , o] <- db2$E.day

db2$data    <- as.Date(db2$julian-1, origin=as.Date("2014-01-01"))
db2$monthF  <- factor(db2$month)
db2$rainF   <- factor(db2$mmrn>0)
db2$tm2mC   <- db2$tm2m -mean(db2$tm2m)
db2$tm2m7C  <- db2$tm2m7 -mean(db2$tm2m7)
db2$tm2m14C <- db2$tm2m14 -mean(db2$tm2m14)
db2$tm2m21C <- db2$tm2m21 -mean(db2$tm2m21)
db2$windC   <- db2$wind -mean(db2$wind)
db2$mmrnC   <- db2$mmrn -mean(db2$mmrn)
db2$mmrn7C  <- db2$mmrn7 -mean(db2$mmrn7)
db2$mmrn14C <- db2$mmrn14 -mean(db2$mmrn14)
db2$mmrn21C <- db2$mmrn21 -mean(db2$mmrn21)

if(length(unique(db2$Site)) == 2){
modelII <- gls(meanfem ~ 1 +E.day + tm2mC + I(tm2mC^2) +mmrn14C +windC+Site,
            correlation=corCAR1(form=~julian|Site), method = "REML",
            na.action = "na.fail",data=db2)}
else{
  modelII <- gls(meanfem ~ 1 +E.day + tm2mC + I(tm2mC^2) +mmrn14C +windC,
            correlation=corCAR1(form=~julian), method = "REML",
            na.action = "na.fail",data=db2)}

modval[ s , 1:2 , o] <- summary(modelII)$tTable[2,c(1,4)]


error <- (resid(modelII,type="response"))^2
modval[ s , 3 , o] <- sqrt(sum(error)/length(error)) #RSME

}}

save(modval,file = "modval_sample_replacment_v2.RData")
save(corval,file = "corval_sample_replacment_v2.RData")


#load(file = "modval_sample_replacment.RData")
#load(file = "corval_sample_replacment.RData")

boxplot(modval[,1,1])
boxplot(modval[,1,2])
boxplot(modval[,1,3])

meanEff <- apply(modval[,1,], 2,median)
lowEff <- apply(modval[,1,], 2,quantile,0.025)
highEff <- apply(modval[,1,], 2,quantile,0.975)
plot(1:15,meanEff,ylim=c(-0.1,0.5),pch=19,
     xlab="N° of Ovitrap", ylab = "Mean number/eggs/day Parameter")
segments(x0 = 1:15, y0=lowEff,x1 = 1:15,y1=highEff)
lines(1:15,lowEff)
lines(1:15,highEff)



meanpv <- apply(modval[,2,], 2,mean)
lowpv <- apply(modval[,2,], 2,quantile,0.025)
highpv <- apply(modval[,2,], 2,quantile,0.975)
plot(1:15,meanpv,ylim=c(0.,0.5))
segments(x0 = 1:15, y0=lowpv,x1 = 1:15,y1=highpv)
lines(1:15,highpv)



meanpv2 <- apply(modval[,2,]<0.05, 2,mean)
plot(1:15,meanpv2*100,ylim=c(0.,100),pch=19,
     xlab="N° of Ovitrap", ylab = "% Significative Parameter")



meanerr <- apply(modval[,3,], 2,mean)
lowerr <- apply(modval[,3,], 2,quantile,0.025)
higherr <- apply(modval[,3,], 2,quantile,0.975)
plot(1:15,meanerr,ylim=c(6,10),pch=19,
     xlab="N° of Ovitrap", ylab = "RMSE")
segments(x0 = 1:15, y0=lowerr,x1 = 1:15,y1=higherr)
lines(1:15,lowerr)
lines(1:15,higherr)



meancor <- apply(corval[,1,], 2,mean)
lowcor <- apply(corval[,1,], 2,quantile,0.025)
highcor <- apply(corval[,1,], 2,quantile,0.975)
plot(1:15,meancor,ylim=c(0,1),pch=19,
     xlab="N° of Ovitrap", ylab = "Pearson Corr. Estimate")
segments(x0 = 1:15, y0=lowcor,x1 = 1:15,y1=highcor)
lines(1:15,lowcor)
lines(1:15,highcor)




meancorp <- apply(corval[,2,]<0.05, 2,mean)
lowcorp <- apply(corval[,2,], 2,quantile,0.025)
highcorp <- apply(corval[,2,], 2,quantile,0.975)
plot(1:15,meancorp,ylim=c(0,1),pch=19,
     xlab="N° of Ovitrap", ylab = "% Significative Corr. Estimater")
lines(1:15,lowcorp)
lines(1:15,highcorp)



#png("FigureSampleRep_v3.png",width=600, height=600,units="px",bg="white")
png("FigureSampleRep_v3.png",res = 600,width=18, height=16,units="cm",bg="white")
m <- rbind(c(1, 2), c(3,4))
layout(m)
par(cex.main = 1.5, mar = c(2, 4, 0.5, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(1:15,meanEff,ylim=c(0,0.5),pch=19,
     xlab="", ylab = "")
segments(x0 = 1:15, y0=lowEff,x1 = 1:15,y1=highEff)
mtext("N° of ovitraps", side = 1, line = 3, cex = 1, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Mean number/eggs/day parameter", side = 2, line = 3, cex = 1, font = 2)

plot(1:15,meanpv2*100,ylim=c(0.,100),pch=19,
     xlab="", ylab = "")
mtext("N° of ovitraps", side = 1, line = 3, cex = 1, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("% Significant parameter (p<0.05)", side = 2, line = 3, cex = 1, font = 2)
abline(h=80,lty=2)

plot(1:15,meancor,ylim=c(0,1),pch=19,
     xlab="", ylab = "")
segments(x0 = 1:15, y0=lowcor,x1 = 1:15,y1=highcor)
mtext("N° of ovitraps", side = 1, line = 3, cex = 1, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Pearson corr. estimate", side = 2, line = 3, cex = 1, font = 2)

plot(1:15,meanerr,ylim=c(6,10),pch=19,
     xlab="", ylab = "")
segments(x0 = 1:15, y0=lowerr,x1 = 1:15,y1=higherr)
mtext("N° of ovitraps", side = 1, line = 3, cex = 1, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("RMSE", side = 2, line = 3, cex = 1, font = 2)

dev.off()