## graphic function #####
mytheme_classic <- function (base_size = 20, base_family = "") 
{
    theme_bw(base_size = base_size, base_family = base_family) %+replace% 
        theme( axis.line = element_line(colour = "black"), 
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            strip.background = element_rect(colour = "black", 
                size = 1), legend.key = element_blank())
}
palette(c("dodgerblue4","#D55E00")) #56B4F9

## packages 
library(plyr)
library(ggplot2)



# Raw database HLC ####
mypath    <- "C:/Users/Mattia/Documents/LaSapienza/HLC/Submission/Data/"
mypathfig <- "C:/Users/Mattia/Documents/LaSapienza/HLC/Submission/Figure/"

dbhlc0 <- read.csv(paste(mypath,"DB HLC-2014.csv",sep=""),sep=";") #HLC sampling
summary(dbhlc0)
names(dbhlc0)

# Spot within site
dbhlc0$Site         <- dbhlc0$Spot
levels(dbhlc0$Site) <- c("Site-B","Site-B","Site-A","Site-A")
dbhlc0$Site         <- relevel(dbhlc0$Site,ref="Site-A") 

# time
dbhlc0$week   <- as.numeric(strftime(as.POSIXlt(dbhlc0$Date,format="%d/%m/%Y"),origin="2014-01-01",format="%W") )
dbhlc0$julian <- as.numeric(strftime(as.POSIXlt(dbhlc0$Date,format="%d/%m/%Y"),origin="2014-01-01",format="%j") )
dbhlc0$dayMon <- as.numeric(strftime(as.POSIXlt(dbhlc0$Date,format="%d/%m/%Y"),origin="2014-01-01",format="%d") )
dbhlc0$month  <- as.numeric(strftime(as.POSIXlt(dbhlc0$Date,format="%d/%m/%Y"),origin="2014-01-01",format="%m") )

table(dbhlc0$Technicians)
dbhlc1 <- droplevels(subset(dbhlc0,Technicians %in% c("Antonello","Giovanni")))


# Descriptive statistics HLC
sum(dbhlc1$Female_albopictus)
tapply(dbhlc1$Female_albopictus,dbhlc1$Site,mean)
tapply(dbhlc1$Female_albopictus,dbhlc1$Site,sd)/sqrt(tapply(dbhlc1$Female_albopictus,dbhlc1$Site,length))
  
tapply(dbhlc1$Female_albopictus,dbhlc1$Site,max)
max(dbhlc1$julian)-min(dbhlc1$julian)
length(dbhlc1$Female_albopictus[dbhlc1$Female_albopictus==0])
length(dbhlc1$Female_albopictus)

range(dbhlc1$Male_albopictus)
sum(dbhlc1$Male_albopictus)

# temperature #####
dbtempisthlc        <- read.csv(paste(mypath,"temperatura_istantanea_2m_2014.csv",sep=""),sep=";")

dbtempisthlc$julian <- as.numeric(strftime(as.POSIXlt(dbtempisthlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%j") )
dbtempisthlc$month  <- as.numeric(strftime(as.POSIXlt(dbtempisthlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%m") )
dbtempisthlc$day    <- factor(as.numeric(strftime(as.POSIXlt(dbtempisthlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%d") ))


dbtp2mhlc2 <- ddply(dbtempisthlc,c("month","julian"),summarize, 
                  temp = mean(temp,na.rm=T),
                  temp2 = sum(temp,na.rm=T))
#ggplot(dbtp2mhlc2)+geom_line(aes(x=julian,y=temp))+theme_bw()+facet_wrap(~month,scales = "free_x")

# lag
dbhlc1$tm2m   <- rep(NA,nrow(dbhlc1))
dbhlc1$tm2m7  <- rep(NA,nrow(dbhlc1))
dbhlc1$tm2m14 <- rep(NA,nrow(dbhlc1))
dbhlc1$tm2m21 <- rep(NA,nrow(dbhlc1))
dbhlc1$tm2m21early <- rep(NA,nrow(dbhlc1))
for(i in 1:nrow(dbhlc1)){
 a <- which(dbtp2mhlc2$julian==dbhlc1$julian[i])
     dbhlc1$tm2m[i]   <- dbtp2mhlc2[a,"temp"]
      dbhlc1$tm2m7[i]  <- mean(dbtp2mhlc2[(a-6):a,"temp"])
      dbhlc1$tm2m14[i] <- mean(dbtp2mhlc2[(a-13):a,"temp"])
     dbhlc1$tm2m21[i] <- mean(dbtp2mhlc2[(a-20):a,"temp"])
     dbhlc1$tm2m21early[i] <- mean(dbtp2mhlc2[(a-20-7):(a-7),"temp"])
 }

# wind #####
dbwindhlc        <- read.csv(paste(mypath,"Vento.csv",sep=""),sep=";") 
dbwindhlc$julian <- as.numeric(strftime(as.POSIXlt(dbwindhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%j") )
dbwindhlc$month  <- as.numeric(strftime(as.POSIXlt(dbwindhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%m") )
dbwindhlc$day    <- factor(as.numeric(strftime(as.POSIXlt(dbwindhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%d") ))

dbwind2 <- ddply(dbwindhlc,c("month","julian"),summarize, 
                  wind  = mean(vel_ms,na.rm=T),
                  wind2 = sum(vel_ms,na.rm=T))

dbhlc1$wind   <- rep(NA,nrow(dbhlc1))

for(i in 1:nrow(dbhlc1)){
 a <- which(dbwind2$julian==dbhlc1$julian[i])
     dbhlc1$wind[i]   <- mean(dbwind2[a,"wind"])
}

# Rainfall #####
dbprecdhlc        <- read.csv(paste(mypath,"precipitazione_oraria_2014.csv",sep=""),sep=";") #C:\Users\Mattia\Documents\LaSapienza\Attivitą 2014\HLC_Ovi\2014_HLC\HLC_dinamica_Orto_anatomia

dbprecdhlc$julian <- as.numeric(strftime(as.POSIXlt(dbprecdhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%j") )
dbprecdhlc$month  <- as.numeric(strftime(as.POSIXlt(dbprecdhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%m") )
dbprecdhlc$day    <- factor(as.numeric(strftime(as.POSIXlt(dbprecdhlc$date,format="%d/%m/%Y"),origin="2014-01-01",format="%d") ))

dbprecdhlc2 <- ddply(dbprecdhlc,c("month","julian"),summarize, 
                  rain  = mean(misura,na.rm=T) #media oraria
                  )

# sum mm rain
dbhlc1$mmrnH   <- rep(NA,nrow(dbhlc1))

for(i in 1:nrow(dbhlc1)){
 a <- which(dbprecdhlc2$julian==dbhlc1$julian[i])
     dbhlc1$mmrnH[i]   <- dbprecdhlc2[a,"rain"]}

dbhlc1$mmrn   <- rep(NA,nrow(dbhlc1))
dbhlc1$mmrn7  <- rep(NA,nrow(dbhlc1))
dbhlc1$mmrn14 <- rep(NA,nrow(dbhlc1))
dbhlc1$mmrn21 <- rep(NA,nrow(dbhlc1))
dbhlc1$mmrn21early <- rep(NA,nrow(dbhlc1))
for(i in 1:nrow(dbhlc1)){
 a <- which(dbprecdhlc2$julian==dbhlc1$julian[i])
     dbhlc1$mmrn[i]   <- dbprecdhlc2[a,"rain"]*24
     dbhlc1$mmrn7[i]  <- sum(dbprecdhlc2[(a-6):a,"rain"]*24)
     dbhlc1$mmrn14[i] <- sum(dbprecdhlc2[(a-13):a,"rain"]*24)
     dbhlc1$mmrn21[i] <- sum(dbprecdhlc2[(a-20):a,"rain"]*24)
     dbhlc1$mmrn21early[i] <- sum(dbprecdhlc2[(a-20-7):(a-7),"rain"]*24)
}

# Ovitraps #####
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") )

dbdis <- subset(dbovi,Date == "01/10/2014" & Site =="Anatomia")
library(raster)
dst <- pointDistance(dbdis[-7,c("long","lat")], longlat=TRUE)
# coerce to dist object
dst <- as.dist(dst)

range(dst)
mean(dst)
hist(sort(dst))


dbdis <- subset(dbovi,Date == "01/09/2014" & Site =="Orto")
library(raster)
dst <- pointDistance(dbdis[,c("long","lat")], longlat=TRUE)
# coerce to dist object
dst <- as.dist(dst)

range(dst)
mean(dst)
hist(sort(dst))
#plot(dbovi$julian,dbovi$NrEggs)

db1 <- droplevels(subset(dbovi, is.na(NrEggs)==FALSE))

levels(db1$Site) <- c("Site-B","Site-A")
db1$Site         <- relevel(db1$Site,ref="Site-A")


# descittive totale uova albopictus
 
sum(db1$NrEggs,na.rm=TRUE)
tapply(db1$NrEggs,db1$Site,mean,na.rm=TRUE)
tapply(db1$NrEggs,db1$Site,max,na.rm=TRUE)
tapply(db1$NrEggs,db1$Site,sd,na.rm=TRUE)/sqrt(tapply(db1$NrEggs,db1$Site,length))


tapply(db1$NrEggs[db1$NrEggs==0],db1$Site[db1$NrEggs==0],length)/tapply(db1$NrEggs,db1$Site,length)
round(length(db1$NrEggs[db1$NrEggs==0])/length(db1$NrEggs)*100,2)
length(db1$NrEggs[db1$NrEggs==0])

VMR = tapply(db1$NrEggs,db1$Site,var,na.rm=TRUE)/tapply(db1$NrEggs,db1$Site,mean,na.rm=TRUE)
VMR

idv = (tapply(db1$NrEggs[db1$Site=="Site-A"],droplevels(db1$OvitrapID[db1$Site=="Site-A"]),var,na.rm=TRUE))
Lidv = log(idv)
idm = tapply(db1$NrEggs[db1$Site=="Site-A"],droplevels(db1$OvitrapID[db1$Site=="Site-A"]),mean,na.rm=TRUE)
mvrm <- lm(Lidv~idm)
summary(mvrm)
plot(mvrm)

a = 5.014
b = 0.066
N = (1.96/ 0.1)^2 * a * mean(db1$NrEggs[db1$Site=="Site-A"])^(b-2)
N



idv = (tapply(db1$NrEggs[db1$Site=="Site-B"],droplevels(db1$OvitrapID[db1$Site=="Site-B"]),var,na.rm=TRUE))
Lidv = log(idv)
idm = tapply(db1$NrEggs[db1$Site=="Site-B"],droplevels(db1$OvitrapID[db1$Site=="Site-B"]),mean,na.rm=TRUE)
mvrm <- lm(Lidv~idm)
summary(mvrm)
plot(mvrm)

a = 6.188
b = 0.032
N = (1.96/ 0.25)^2 * a * mean(db1$NrEggs[db1$Site=="Site-B"])^(b-2)
N

tapply(db1$NrEggs,db1$Site,sd,na.rm=TRUE)/tapply(db1$NrEggs,db1$Site,mean,na.rm=TRUE)

### HLC and ovitrap relationship######

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))

#plot(ovitrapp$julian,ovitrapp$uova)
range(ovitrapp$julian)
day <- 189:304
ovitrapp$tempo


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

# Eggs lag ###

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)
}


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)
}

dayA <- ifelse(is.na(mdayA),0,1)
dayB <- ifelse(is.na(mdayB),0,1)


dbhlc$E.day    <- rep(NA,nrow(dbhlc))
dbhlc$E.week1  <- rep(NA,nrow(dbhlc))
dbhlc$E.week2 <- 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}
    tt <- sum(dayA[dd])
    dbhlc$E.week1[i] <- sum(mdayA[dd],na.rm=T)/tt
    tt2 <- sum(dayA[dd2])
    dbhlc$E.week2[i] <- sum(mdayA[dd2],na.rm=T)/tt2
    
  }else{
    if(dbhlc$julian[i] %in% ovitrappB$julian){dbhlc$E.day[i] <- mdayB[day == dbhlc$julian[i]]}
    else{dbhlc$E.day[i] <-NA}    
    tt <- sum(dayB[dd])
    dbhlc$E.week1[i] <- sum(mdayB[dd],na.rm=T)/tt
    tt2 <- sum(dayB[dd2])
    dbhlc$E.week2[i] <- sum(mdayB[dd2],na.rm=T)/tt2
  }
}


plot(dbhlc$E.day,dbhlc$E.week1)
plot(dbhlc$E.day,dbhlc$E.week2)
plot(dbhlc$E.week1,dbhlc$E.week2)


# merge HLC, Egg and climate dataset
db2        <- subset(dbhlc, is.na(E.day)==FALSE&is.na(E.week2)==FALSE )

db2$data   <- as.Date(db2$julian-1, origin=as.Date("2014-01-01"))

# center variable
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$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)
db2$windC   <- db2$wind -mean(db2$wind)


# bimodal pattern
library(mgcv)
gam1 <- gam(meanfem ~ s(julian), data=db2)
plot(gam1);summary(gam1)
gam2 <- gam(E.day ~ s(julian), data=db2)
plot(gam2);summary(gam2)
detach("package:mgcv")


# correlation ####
cor.test(db2$meanfem,db2$E.day)
cor.test(db2$meanfem,db2$E.week1)
cor.test(db2$meanfem,db2$E.week2)


#MODEL I #####
source("HighstatLibV10.R")
myv <- c("E.day","E.week1","E.week2","Site","rainF","tm2mC","windC")
corvif(db2[,myv])

library(nlme)
library(MuMIn)

first.modelI <- gls(meanfem ~ E.day+E.week1+E.week2+Site+rainF+tm2mC+I(tm2mC^2)+windC,
                   na.action = "na.fail",data=db2)
summary(first.modelI)
acf(resid(first.modelI,type="normalized"),na.action=na.pass)



full.modelI <- gls(meanfem ~ E.day+E.week1+E.week2+Site+rainF+tm2mC+I(tm2mC^2)+windC, 
                   correlation=corCAR1(form=~julian|Site), method = "ML",
                   na.action = "na.fail",data=db2)
summary(full.modelI)
acf(resid(full.modelI,type="normalized"),na.action=na.pass)
plot(resid(full.modelI,type="normalized"),pch=19)
plot(fitted(full.modelI),resid(full.modelI,type="normalized"),pch=19);min(fitted(full.modelI))
plot(db2$julian,resid(full.modelI,type="normalized"),pch=19)
plot(db2$wind,resid(full.modelI,type="normalized"),pch=19)
plot(db2$tm2mC,resid(full.modelI,type="normalized"),pch=19)
plot(db2$Site,resid(full.modelI,type="normalized"),pch=19)
plot(db2$E.day,resid(full.modelI,type="normalized"),pch=19)
plot(db2$E.week1,resid(full.modelI,type="normalized"),pch=19)




msubset <- expression(tm2mC | !`I(tm2mC^2)`) 
mod <- dredge(full.modelI,  rank = "AICc",subset=msubset)# Get the models 
nrow(mod)
importance(mod)
importance(get.models(mod,subset = delta < 2))

get.models(mod,1)

write.csv(as.data.frame(mod),file="Table_s1_ModelI.csv",quote=FALSE)
# Model average models with delta AICc < 4
modelavg_mod<-model.avg(mod, subset = delta < 2)
modelavg_mod
summary(modelavg_mod)


modelI <- gls(meanfem ~ E.day+windC,
            correlation=corCAR1(form=~julian|Site), method = "REML",
            na.action = "na.fail",data=db2)
summary(modelI)
acf(resid(modelI,type="normalized"),na.action=na.pass)
plot(resid(modelI,type="normalized"),pch=19)
plot(fitted(modelI),resid(modelI,type="normalized"),pch=19);min(fitted(modelI));abline(h=0)
plot(db2$julian,resid(modelI,type="normalized"),pch=19)
plot(db2$wind,resid(modelI,type="normalized"),pch=19)
plot(db2$tm2mC,resid(modelI,type="normalized"),pch=19)
plot(db2$Site,resid(modelI,type="normalized"),pch=19)
plot(db2$E.day,resid(modelI,type="normalized"),pch=19)
plot(db2$E.week1,resid(modelI,type="normalized"),pch=19)

#RMSE
error <- (resid(modelI,type="response"))^2
sqrt(sum(error)/length(error)) #RSME
cor.test(db2$meanfem,as.numeric(fitted(modelI)))
plot(db2$meanfem,as.numeric(fitted(modelI)))

detach("package:MuMIn")

library(AICcmodavg)

db.fig <- data.frame(E.day = seq(0,60,length=100),
                     windC = 0)

pr <- predictSE.gls(modelI, newdata = db.fig, se.fit = TRUE, print.matrix = FALSE)
db.fig$pred <-  pr$fit 
#XV <- model.matrix(~E.daymean,data=db.fig)
#db.fig$pred <-  XV %*% coef(x1.4) 
db.fig$predseup <-  pr$fit +1.96*pr$se.fit
db.fig$predselo <-  pr$fit -1.96*pr$se.fit

prraw0 <- predict(modelI,newdata=db2)
db2$pred <-  prraw0
prraw <- predictSE.gls(modelI, newdata = db2, se.fit = TRUE, print.matrix = FALSE)
plot(prraw0,prraw$fit)
db2$pred2 <-  prraw$fit 
db2$predseup2 <-  prraw$fit +1.96*prraw$se.fit
db2$predselo2 <-  prraw$fit -1.96*prraw$se.fit


Sys.setlocale("LC_TIME", "English")
db2$CommonDate <- as.Date((db2$julian-1),origin = "2014-1-1")
cd <- c(as.Date("2014-06-30"),db2$CommonDate)
db2$MonthDay <- as.Date(strftime(as.POSIXlt(db2$CommonDate,format="%Y-%m-%d"),origin="2014-01-01",format="%b-%d"),format="%b-%d")


# Figure 2 ##########
png(paste(mypathfig,"Figure2.png",sep=""),width=900, height=900,units="px",bg="white")


dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))

m <- rbind(c(1, 2), c(3,4))
layout(m)
#layout.show(3)
par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(db.fig$E.day,db.fig$pred, type = "l", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  xlim=c(0,65),lwd = 2, pch = 5, axes = T, main = " ")
mtext("Mean eggs/day", side = 1, line = 3, cex = 2, font = 2)

mtext("A", at=48,side = 2, line = 3, cex = 3, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2)
points(db2$E.day, db2$meanfem,col=db2$Site, cex = 2,  pch = 19)
text(3, 50, "Site-A", cex = 2, font = 1, adj = 0)
points(2, 50, col="dodgerblue4",cex = 2, pch = 19)
text(3, 45, "Site-B", cex = 2, font = 1, adj = 0)
points(2, 45,col="#D55E00", cex = 2, pch = 19)
text(as.Date("2014-07-05"), 40, expression(paste("Estimated")), cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex = 2, pch = 3)
lines(db.fig$E.day,db.fig$pred)
lines(db.fig$E.day,db.fig$predseup,lty=2)
lines(db.fig$E.day,db.fig$predselo,lty=2)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(db2$meanfem,db2$pred, type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  xlim=c(0,50),lwd = 2, pch = 1, axes = T, main = " ")
mtext("Observed HLC", side = 1, line = 3, cex = 2, font = 2)

mtext("B", at=48,side = 2, line = 3, cex = 3, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Fitted HLC", side = 2, line = 3, cex = 2, font = 2)
abline(0,1,lty=2)
text(2, 45, expression(paste("RMSE = 8.9")), cex = 2, font = 1, adj = 0)
text(2, 40, expression(paste("Pearson Corr = 0.53")), cex = 2, font = 1, adj = 0)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  lwd = 2, pch = 5, axes = T, main = " ")
par(las = 1)
mtext("C", at=50,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2,las=3)
points(dbA$CommonDate, dbA$meanfem, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 45, "Observed", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 45, cex =2, pch = 19)
text(as.Date("2014-07-05"), 40, "Estimated", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex =2, pch = 3)
points(dbA$CommonDate, dbA$pred, cex = 2, pch = 3)
segments(x0=dbA$CommonDate,x1=dbA$CommonDate, y0=dbA$predselo, y1=dbA$predseup, cex = 1, lwd = 2)
mtext("", side = 1, line = 3, cex = 2, font = 2)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  lwd = 2, pch = 5, axes = T, main = " ")
#axis(1, at = seq(202, 302,by=7), labels = seq(202, 302,by=7))
mtext("D", at=50,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2)

mtext("", side = 1, line = 3, cex = 2, font = 2)
axis(2, pos = 200 )
par(las = 0)
#mtext("Mean biting adults", side = 2, line = 3, cex = 1.5, font = 2)
points(dbB$CommonDate, dbB$meanfem, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, "Site-B", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 45, "Observed", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 45, cex = 2, lwd = 2, pch = 19)
text(as.Date("2014-07-05"), 40, "Estimated", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex = 2, lwd = 2, pch = 3)
#text(12.5, 0.98, expression(paste("N° Specimens = ", 5)), cex = 1.2, font = 1, adj = 0)
#points(dbB$CommonDate, dbB$mu, cex = 1, lwd = 1,  pch = 3)
#points(dbB$CommonDate, dbB$mu, cex = 1, lwd = 2)
#segments(x0=dbB$CommonDate,x1=dbB$CommonDate, y0=dbB$mul, y1=dbB$muh, cex = 1, lwd = 2)

points(dbB$CommonDate, dbB$pred, cex =2, lwd = 2, pch = 3)
segments(x0=dbB$CommonDate,x1=dbB$CommonDate, y0=dbB$predselo, y1=dbB$predseup, cex = 1, lwd = 2)

dev.off()

detach("package:AICcmodavg")


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

# correlation ####
cor.test(db2$meanfem,db2$E.day)
cor.test(db2$meanfem,db2$E.week1)
cor.test(db2$meanfem,db2$E.week2)


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)

# Figure 1 #######
Sys.setlocale("LC_TIME", "English")
db2$CommonDate <- as.Date((db2$julian-1),origin = "2014-1-1")
cd <- c(as.Date("2014-06-30"),db2$CommonDate)
db2$MonthDay <- as.Date(strftime(as.POSIXlt(db2$CommonDate,format="%Y-%m-%d"),origin="2014-01-01",format="%b-%d"),format="%b-%d")

dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))

# Figure 1 #######
png(paste(mypathfig,"Figure1.png",sep=""),width=900, height=900,units="px",bg="white")

par(mfrow=c(2,1))
par(cex.main = 1.5, mar = c(4, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 1.5, 
     ylim = c(0, 60),  lwd = 2, pch = 5, axes = T, main = " ")

par(las = 0)
mtext(expression(italic("Aedes albopictus")), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate, dbA$meanfem, cex = 2, pch = 19)
lines(dbA$CommonDate[order(dbA$CommonDate)], dbA$meanfem[order(dbA$CommonDate)], lty = 1, lwd = 1)
#segments(x0=dbA$CommonDate,x1=dbA$CommonDate,y0=dbA$selo,y1=dbA$seup)
text(as.Date("2014-07-05"), 60, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 55, expression(paste("Mean biting adults")), cex = 1.8, font = 1, adj = 0)
points(as.Date("2014-07-03"), 55, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, expression(paste("Mean eggs/day")), cex = 1.8, font = 1, adj = 0)
points(as.Date("2014-07-03"), 50, cex = 2, pch = 3)
points(dbA$CommonDate, dbA$E.day, cex = 2,  pch = 3)
lines(dbA$CommonDate[order(dbA$CommonDate)], dbA$E.day[order(dbA$CommonDate)], lty = 2, lwd = 1)

par(cex.main = 1.5, mar = c(4, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 1.5, 
     ylim = c(0, 60),  lwd = 2, pch = 5, axes = T, main = " ")
#axis(1, at = seq(202, 302,by=7), labels = seq(202, 302,by=7))

mtext("Year 2014", side = 1, line = 3, cex = 2, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext(expression(italic("Aedes albopictus")), side = 2, line = 3, cex = 2, font = 2)
points(dbB$CommonDate, dbB$meanfem, cex = 2,  pch = 19)
lines(dbB$CommonDate[order(dbB$CommonDate)], dbB$meanfem[order(dbB$CommonDate)], lty = 1, lwd = 1)
text(as.Date("2014-07-05"), 60, "Site-B", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 55, expression(paste("Mean biting adults")), cex = 1.8, font = 1, adj = 0)
points(as.Date("2014-07-03"), 55, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, expression(paste("Mean eggs/day")), cex = 1.8, font = 1, adj = 0)
points(as.Date("2014-07-03"), 50, cex = 2, pch = 3)
#text(12.5, 0.98, expression(paste("N° Specimens = ", 5)), cex = 1.2, font = 1, adj = 0)
points(dbB$CommonDate[order(dbB$CommonDate)], dbB$E.day[order(dbB$CommonDate)], cex = 2, lwd = 1,  pch = 3)
lines(dbB$CommonDate[order(dbB$CommonDate)], dbB$E.day[order(dbB$CommonDate)], lty = 2, lwd = 1)

dev.off()


library(MuMIn)
library(nlme)
vartmp <-  c("+ tm2mC ",
             "+ tm2m7C ",
             "+ tm2m14C ",
             "+ tm2m21C ",
             "+ tm2mC + I(tm2mC^2) ",
             "+ tm2m7C+ I(tm2m7C^2) ",
             "+ tm2m14C+ I(tm2m14C^2) ",
             "+ tm2m21C+ I(tm2m21C^2) ",
             "")
varwnd <-  c("+windC ", "")
varegg <-  c("+E.day +Site","+E.day + E.day:Site+Site","+E.day ","+Site","")
varain <-  c("+ rainF ","+ mmrn7C ","+ mmrn14C ","+ mmrn21C",
             "+ mmrn7C + I(mmrn7C^2) ","+ mmrn14C + I(mmrn14C^2) ",
             "+ mmrn21C + I(mmrn21C^2)", "")

dbestE <- data.frame(formula = rep(NA,length(varegg)*length(vartmp)*length(varain)*length(varwnd)),
                     AICc     = rep(NA,length(varegg)*length(vartmp)*length(varain)*length(varwnd)),
                     rsme     = rep(NA,length(varegg)*length(vartmp)*length(varain)*length(varwnd)))

x = 1
for(i in 1:length(vartmp)){
  for(z in 1:length(varwnd)){
    for(j in 1:length(varain)){
      for(q in 1:length(varegg)){
        f1  <- as.formula(paste("meanfem ~  1",varegg[q],vartmp[i],varain[j],varwnd[z],sep=""))
        b1 <- gls(f1, correlation = corCAR1(form=~julian|Site), 
                  #weights = varIdent(form=~1|monthF),
                  method="ML",
                  data=db2)
        dbestE$formula[x] <- paste("meanfem ~ 1 ",varegg[q],vartmp[i],varain[j],varwnd[z],sep="")
        dbestE$AICc[x]     <- AICc(b1)
        b2 <- gls(f1, correlation = corCAR1(form=~julian|Site), 
                  method="REML",                  data=db2)
        error <- (resid(b2,type="response"))^2
        dbestE$rsme[x]     <-sqrt(sum(error)/length(error))
        x <- x+1
      }
    }
  }
}


head(dbestE[order(dbestE$AICc),],5)
head(dbestE[order(dbestE$rsme),],5)
write.csv(dbestE[order(dbestE$AICc),],file="Table_s2_ModelII.csv",quote=FALSE)

source("HighstatLibV10.R")
myv <- c("E.day","mmrn14C","tm2mC","windC","Site")
corvif(db2[,myv])

modelII <- gls(meanfem ~ 1 +E.day + tm2mC + I(tm2mC^2) +mmrn14C +windC+Site,#+temp+I(temp^2)I(tempC^2)+
            correlation=corCAR1(form=~julian|Site), method = "REML",
            na.action = "na.fail",data=db2)
summary(modelII)
intervals(modelII)
detach("package:MuMIn")
acf(resid(modelII,type="normalized"),na.action=na.pass)
plot(resid(modelII,type="normalized"),pch=19)
plot(fitted(modelII),resid(modelII,type="normalized"),pch=19);min(fitted(modelII));abline(h=0)
plot(db2$julian,resid(modelII,type="normalized"),pch=19)
plot(db2$wind,resid(modelII,type="normalized"),pch=19)
plot(db2$tm2mC,resid(modelII,type="normalized"),pch=19)
plot(db2$Site,resid(modelII,type="normalized"),pch=19)
plot(db2$E.day,resid(modelII,type="normalized"),pch=19)
plot(db2$mmrn14C,resid(modelII,type="normalized"),pch=19)



error <- (resid(modelII,type="response"))^2
sqrt(sum(error)/length(error)) #RSME
#2*sqrt(sum(error)/length(error))/(max
cor.test(db2$meanfem,as.numeric(fitted(modelII)))
plot(db2$meanfem,as.numeric(fitted(modelII)))

library(AICcmodavg)

db.fig <- data.frame(E.day = seq(0,60,length=100),
                   tm2mC     = 0,
                   mmrn14C   = 0,
                   windC     = 0,
                   Site      = "Site-A")

pr <- predictSE.gls(modelII, newdata = db.fig, se.fit = TRUE, print.matrix = FALSE)
db.fig$pred <-  pr$fit 
#XV <- model.matrix(~E.daymean,data=db.fig)
#db.fig$pred <-  XV %*% coef(x1.4) 
db.fig$predseup <-  pr$fit +1.96*pr$se.fit
db.fig$predselo <-  pr$fit -1.96*pr$se.fit

prraw0 <- predict(modelII,newdata=db2)
db2$pred <-  prraw0
prraw <- predictSE.gls(modelII, newdata = db2, se.fit = TRUE, print.matrix = FALSE)
plot(prraw0,prraw$fit)
db2$pred2 <-  prraw$fit 
db2$predseup2 <-  prraw$fit +1.96*prraw$se.fit
db2$predselo2 <-  prraw$fit -1.96*prraw$se.fit


# FIgure 3 #######
png(paste(mypathfig, "Figure3.png",sep=""),width=900, height=900,units="px",bg="white")


dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))

m <- rbind(c(1, 2), c(3,4))
layout(m)
#layout.show(3)
par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(db.fig$E.day,db.fig$pred,group=db2$Site, type = "l", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  xlim=c(0,65),lwd = 2, pch = 5, axes = T, main = " ")
mtext("Mean eggs/day", side = 1, line = 3, cex = 2, font = 2)

mtext("A", at=48,side = 2, line = 3, cex = 3, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2)
points(db2$E.day, db2$meanfem,col=db2$Site, cex = 2, pch = 19)
text(3, 50, "Site-A", cex = 2, font = 1, adj = 0)
points(2, 50, col="dodgerblue4",cex = 2, pch = 19)
text(3, 45, "Site-B", cex = 2, font = 1, adj = 0)
points(2, 45,col="#D55E00", cex = 2, pch = 19)
text(as.Date("2014-07-05"), 40, expression(paste("Estimated")), cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex = 2, pch = 3)
lines(db.fig$E.day,db.fig$pred)
lines(db.fig$E.day,db.fig$predseup,lty=2)
lines(db.fig$E.day,db.fig$predselo,lty=2)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(db2$meanfem,db2$pred, type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  xlim=c(0,50),lwd = 2, pch = 1, axes = T, main = " ")
mtext("Observed HLC", side = 1, line = 3, cex = 2, font = 2)

mtext("B", at=48,side = 2, line = 3, cex = 3, font = 2)
axis(2, pos = 200 )
par(las = 0)
mtext("Fitted HLC", side = 2, line = 3, cex = 2, font = 2)
abline(0,1,lty=2)
text(2, 45, expression(paste("RMSE = 6.9")), cex = 2, font = 1, adj = 0)
text(2, 40, expression(paste("Pearson Corr = 0.76")), cex = 2, font = 1, adj = 0)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  lwd = 2, pch = 5, axes = T, main = " ")
par(las = 1)
mtext("C", at=50,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2,las=3)
points(dbA$CommonDate, dbA$meanfem, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 45, "Observed", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 45, cex =2, pch = 19)
text(as.Date("2014-07-05"), 40, "Estimated", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex = 2, pch = 3)
points(dbA$CommonDate, dbA$pred, cex = 2, pch = 3)
#points(dbA$CommonDate, dbA$mu, cex = 1, lwd = 2)
segments(x0=dbA$CommonDate,x1=dbA$CommonDate, y0=dbA$predselo, y1=dbA$predseup, cex = 1, lwd = 2)
mtext("", side = 1, line = 3, cex = 2, font = 2)

par(cex.main = 1.5, mar = c(5, 5, 2, 2) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 2, 
    ylim = c(0, 50),  lwd = 2, pch = 5, axes = T, main = " ")
#axis(1, at = seq(202, 302,by=7), labels = seq(202, 302,by=7))
mtext("D", at=50,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext("Mean biting adults", side = 2, line = 3, cex = 2, font = 2,las=3)

mtext("", side = 1, line = 3, cex = 2, font = 2)
axis(2, pos = 200 )
par(las = 0)
#mtext("Mean biting adults", side = 2, line = 3, cex = 1.5, font = 2)
points(dbB$CommonDate, dbB$meanfem, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 50, "Site-B", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 45, "Observed", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 45, cex = 2, pch = 19)
text(as.Date("2014-07-05"), 40, "Estimated", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 40, cex = 2, pch = 3)
#text(12.5, 0.98, expression(paste("N° Specimens = ", 5)), cex = 1.2, font = 1, adj = 0)
#points(dbB$CommonDate, dbB$mu, cex = 1, lwd = 1,  pch = 3)
#points(dbB$CommonDate, dbB$mu, cex = 1, lwd = 2)
#segments(x0=dbB$CommonDate,x1=dbB$CommonDate, y0=dbB$mul, y1=dbB$muh, cex = 1, lwd = 2)

points(dbB$CommonDate, dbB$pred, cex = 1, lwd = 2, pch = 3)
segments(x0=dbB$CommonDate,x1=dbB$CommonDate, y0=dbB$predselo, y1=dbB$predseup, cex = 1, lwd = 2)


dev.off()
# trap simulation






# r0 #######
library(AICcmodavg)
prraw <- predictSE.gls(modelII, newdata = db2, se.fit = TRUE, print.matrix = FALSE)
db2$pred <-  prraw$fit 
db2$predseup <-  prraw$fit +1.96*prraw$se.fit
db2$predselo <-  prraw$fit -1.96*prraw$se.fit

kVH.obs    = db2$meanfem
kVH.obsH   = db2$seup
kVH.obsL   = ifelse(db2$selo<0,0,db2$selo)

k = 0.09
xv =0.85
xh =0.65
g = 1/4.5
wv = 1/2.5
m = 4*(0.031 + 95820*exp(db2$tm2m-50.4))

db2$r0.obs_M    = sqrt(0.101*kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))

range(db2$r0.obs_M )

r0hv.obs_M    = 0.101*kVH.obs*xv*wv/(g*(wv+m))
r0hv.obsH_M   = 0.101*kVH.obsH*xv*wv/(g*(wv+m))
r0hv.obsL_M   = 0.101*kVH.obsL*xv*wv/(g*(wv+m))

r0vh.obs_M    = k*xh/m
r0vh.obsH_M   = k*xh/m
r0vh.obsL_M   = k*xh/m

db2$outB  <- (1 - ( (r0vh.obs_M + 1)/(r0vh.obs_M*(r0hv.obs_M + 1) ) ))*100
db2$outBH <- (1 - ( (r0vh.obsH_M + 1)/(r0vh.obsH_M*(r0hv.obsH_M + 1) ) ))*100
db2$outBL <- (1 - ( (r0vh.obsL_M + 1)/(r0vh.obsL_M*(r0hv.obsL_M + 1) ) ))*100

db2$outB <- ifelse(db2$outB<0,0,db2$outB)
db2$outBL <- ifelse(db2$outBL<0,0,db2$outBL)
db2$outBH <- ifelse(db2$outBH<0,0,db2$outBH)


Sys.setlocale("LC_TIME", "English")
db2$CommonDate <- as.Date((db2$julian-1),origin = "2014-1-1")
cd <- c(as.Date("2014-06-30"),db2$CommonDate)
db2$MonthDay <- as.Date(strftime(as.POSIXlt(db2$CommonDate,format="%Y-%m-%d"),origin="2014-01-01",format="%b-%d"),format="%b-%d")

dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))

png(paste(mypathfig,"Figure5.png",sep=""),width=900, height=900,units="px",bg="white")

ml <- rbind(c(1, 2), c(3, 4))
layout(ml)
#layout.show(4)

par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 100),  lwd = 2, pch = 5, axes = T, main = "Chikungunya risk by observed data")
mtext("A", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic("Outbreak Risk")), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$outB, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$outBL,y1=dbA$outBH,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$outB, cex =2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$outBL,y1=dbB$outBH,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=20,col="grey75",lty=2)
abline(h=40,col="grey75",lty=2)
abline(h=60,col="grey75",lty=2)
abline(h=80,col="grey75",lty=2)


kVH.obs    = db2$pred
kVH.obsH   = db2$predseup
kVH.obsL   = ifelse(db2$predselo>0,db2$predselo,0)

r0hv.obs_M    = 0.101*kVH.obs*xv*wv/(g*(wv+m))
r0hv.obsH_M   = 0.101*kVH.obsH*xv*wv/(g*(wv+m))
r0hv.obsL_M   = 0.101*kVH.obsL*xv*wv/(g*(wv+m))

r0vh.obs_M    = k*xh/m
r0vh.obsH_M   = k*xh/m
r0vh.obsL_M   = k*xh/m

db2$outB  <- (1 - ( (r0vh.obs_M + 1)/(r0vh.obs_M*(r0hv.obs_M + 1) ) ))*100
db2$outBH <- (1 - ( (r0vh.obsH_M + 1)/(r0vh.obsH_M*(r0hv.obsH_M + 1) ) ))*100
db2$outBL <- (1 - ( (r0vh.obsL_M + 1)/(r0vh.obsL_M*(r0hv.obsL_M + 1) ) ))*100

db2$outB <- ifelse(db2$outB<0,0,db2$outB)
db2$outBL <- ifelse(db2$outBL<0,0,db2$outBL)
db2$outBH <- ifelse(db2$outBH<0,0,db2$outBH)

dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 100),  lwd = 2, pch = 5, axes = T, main = "Chikungunya risk by Model II")
mtext("B", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic("Outbreak Risk")), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$outB, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$outBL,y1=dbA$outBH,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$outB, cex = 2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$outBL,y1=dbB$outBH,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=20,col="grey75",lty=2)
abline(h=40,col="grey75",lty=2)
abline(h=60,col="grey75",lty=2)
abline(h=80,col="grey75",lty=2)


# zika
k = 0.09
xv =0.5#0.0665
xh =0.5
g = 1/5.8
wv = 1/10.5
m = 4*(0.031 + 95820*exp(db2$tm2m-50.4))

kVH.obs    = db2$meanfem
kVH.obsH   = db2$seup
kVH.obsL   = ifelse(db2$selo<0,0,db2$selo)

db2$r0.obs_M    = sqrt(kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))
range(db2$r0.obs_M )

r0hv.obs_M    = 0.101*kVH.obs*xv*wv/(g*(wv+m))
r0hv.obsH_M   = 0.101*kVH.obsH*xv*wv/(g*(wv+m))
r0hv.obsL_M   = 0.101*kVH.obsL*xv*wv/(g*(wv+m))

r0vh.obs_M    = k*xh/m
r0vh.obsH_M   = k*xh/m
r0vh.obsL_M   = k*xh/m

db2$outB  <- (1 - ( (r0vh.obs_M + 1)/(r0vh.obs_M*(r0hv.obs_M + 1) ) ))*100
db2$outBH <- (1 - ( (r0vh.obsH_M + 1)/(r0vh.obsH_M*(r0hv.obsH_M + 1) ) ))*100
db2$outBL <- (1 - ( (r0vh.obsL_M + 1)/(r0vh.obsL_M*(r0hv.obsL_M + 1) ) ))*100

db2$outB <- ifelse(db2$outB<0,0,db2$outB)
db2$outBL <- ifelse(db2$outBL<0,0,db2$outBL)
db2$outBH <- ifelse(db2$outBH<0,0,db2$outBH)
dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 100),  lwd = 2, pch = 5, axes = T, main = "Zika risk by observed data")
mtext("C", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic("Outbreak Risk")), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$outB, cex =2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$outBL,y1=dbA$outBH,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$outB, cex = 2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$outBL,y1=dbB$outBH,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=20,col="grey75",lty=2)
abline(h=40,col="grey75",lty=2)
abline(h=60,col="grey75",lty=2)
abline(h=80,col="grey75",lty=2)

kVH.obs    = db2$pred
kVH.obsH   = db2$predseup
kVH.obsL   = db2$predselo

r0hv.obs_M    = 0.101*kVH.obs*xv*wv/(g*(wv+m))
r0hv.obsH_M   = 0.101*kVH.obsH*xv*wv/(g*(wv+m))
r0hv.obsL_M   = 0.101*kVH.obsL*xv*wv/(g*(wv+m))

r0vh.obs_M    = k*xh/m
r0vh.obsH_M   = k*xh/m
r0vh.obsL_M   = k*xh/m

db2$outB  <- (1 - ( (r0vh.obs_M + 1)/(r0vh.obs_M*(r0hv.obs_M + 1) ) ))*100
db2$outBH <- (1 - ( (r0vh.obsH_M + 1)/(r0vh.obsH_M*(r0hv.obsH_M + 1) ) ))*100
db2$outBL <- (1 - ( (r0vh.obsL_M + 1)/(r0vh.obsL_M*(r0hv.obsL_M + 1) ) ))*100

db2$outB <- ifelse(db2$outB<0,0,db2$outB)
db2$outBL <- ifelse(db2$outBL<0,0,db2$outBL)
db2$outBH <- ifelse(db2$outBH<0,0,db2$outBH)
dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 100),  lwd = 2, pch = 5, axes = T, main = "Zika risk by Model II")
mtext("D", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic("Outbreak Risk")), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$outB, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$outBL,y1=dbA$outBH,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$outB, cex = 2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$outBL,y1=dbB$outBH,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=20,col="grey75",lty=2)
abline(h=40,col="grey75",lty=2)
abline(h=60,col="grey75",lty=2)
abline(h=80,col="grey75",lty=2)
dev.off()




range(db2$tm2m)

#mean monthly climatic value ####
tapply(dbprecdhlc$misura, dbprecdhlc$month,mean,na.rm=T)
rr<-ddply(dbprecdhlc,c("month","julian"),summarise,r = mean(misura,na.rm=TRUE))
tapply(rr$r, rr$month,mean,na.rm=T)
dbprecdhlc2$rain14<-NA
for(i in 15:365){
     dbprecdhlc2$rain14[i] <- sum(dbprecdhlc2[(i-13):i,"rain"]*24)
 }
rr<-tapply(dbprecdhlc2$rain14, dbprecdhlc2$month,mean,na.rm=T)-mean(db2$mmrn14)
rmese <- as.numeric(rr[6:11])



tapply(dbwindhlc$vel_ms, dbwindhlc$month,mean,na.rm=T)
ww<-ddply(dbwindhlc,c("month","julian"),summarise,r = mean(vel_ms,na.rm=TRUE))
wmese <-as.numeric(tapply(ww$r, ww$month,mean,na.rm=T)[6:11]-mean(db2$wind))

tmese <- as.numeric(tapply(dbtempisthlc$temp, dbtempisthlc$month,mean,na.rm=T)[6:11]-mean(db2$tm2m))


dbr0 <- data.frame(E.day = rep(seq(0,100,length=100),6),
                   tm2mC     = rep(tmese,each=100),
                   mmrn14C   = rep(rmese,each=100),
                   windC     = rep(wmese,each=100),
                   Site      = "Site-A",
                   mese      = factor(rep(c("June","July","August","September","October","November"),each=100),
                                      levels = c("June","July","August","September","October","November")))



modelII <- gls(meanfem ~ E.day+ tm2mC+I(tm2mC^2)  + mmrn14C + windC  + Site,
          correlation=corCAR1(form=~julian|Site), 
          method = "REML",na.action = "na.fail",data=db2)
summary(modelII)
error  <- (resid(modelII,type="response"))^2
rmse.b <- sqrt(sum(error)/length(error)) 
rmse.b

b.pr   <- predictSE.gls(modelII, newdata = dbr0, se.fit = TRUE, print.matrix = FALSE)
dbr0$b.fit <-  b.pr$fit 
dbr0$b.seup  <-  b.pr$fit +1.96* 7.572908 #rmse.b #residual standard error
dbr0$b.selo  <-  b.pr$fit -1.96* 7.572908 #rmse.b #residual standard error


b.obs    = dbr0$b.fit
b.obsH   = dbr0$b.seup 
b.obsL   = ifelse(dbr0$b.selo<0,0,dbr0$b.selo)

# mean hlc generated by 0-100 eggs
plot(dbr0$mese,dbr0$b.fit)


k = 0.09
xv =0.85
xh =0.65
g = 1/4.5
wv = 1/2.5
m = 4*(0.031 + 95820*exp(dbr0$tm2mC+mean(db2$tm2m)-50.4))

dbr0$r0.b  = sqrt(0.101*b.obs*k*xv*xh*wv/(g*m*(wv+m)))
dbr0$r0.bh = sqrt(0.101*b.obsH*k*xv*xh*wv/(g*m*(wv+m)))
dbr0$r0.bl = sqrt(0.101*b.obsL*k*xv*xh*wv/(g*m*(wv+m)))

dbr0$val <- NA
jj <- dbr0[dbr0$mese == "June",]
dbr0$val[dbr0$mese == "June"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj <- dbr0[dbr0$mese == "July",]
dbr0$val[dbr0$mese == "July"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj <- dbr0[dbr0$mese == "August",]
dbr0$val[dbr0$mese == "August"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj <- dbr0[dbr0$mese == "September",]
dbr0$val[dbr0$mese == "September"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj <- dbr0[dbr0$mese == "October",]
dbr0$val[dbr0$mese == "October"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj <- dbr0[dbr0$mese == "November",]
dbr0$val[dbr0$mese == "November"] <- round( min(jj$E.day[which(jj$r0.bl>1)]),0)
jj$E.day[which(jj$r0.bl<1)]

tapply(dbr0$val,dbr0$mese,mean)

library(ggplot2)

head(dbr0)


library(grid)
png(paste(mypathfig,"Figure4.png",sep=""),width=900, height=900,units="px",bg="white")

lbx1 = bquote(paste(R[0], ">1 when Mean eggs/day >",28))
lbx2 = bquote(paste(R[0], ">1 when Mean eggs/day >",20))
lbx3 = bquote(paste(R[0], ">1 when Mean eggs/day >",20))
lbx4 = bquote(paste(R[0], ">1 when Mean eggs/day >",3))
lbx5 = bquote(paste(R[0], ">1 when Mean eggs/day >",12))
lbx6 = bquote(paste(R[0], ">1 when Mean eggs/day >",80))

my_grob = grobTree(textGrob(lbx1, x=0.5, y=0.95)) 

ggplot(dbr0,aes(x=E.day,y=r0.b))+
  geom_ribbon(aes(ymin = r0.bl, ymax = r0.bh), fill = "grey80")+
  #geom_line(aes(x=E.day,y=r0.bh),linetype="dashed")+
  #geom_line(aes(x=E.day,y=r0.bl),linetype="dashed")+
  annotation_custom(my_grob)+#annotate("text",x=40,y=2.9,size=5,label=lbx1)+
  mytheme_classic()+facet_wrap(~mese)+geom_line()+
  ylab(expression(R[0]))+xlab("Mean eggs/day")+geom_hline(yintercept=1)
dev.off()

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

# supplementary ##
# r0 #######
library(AICcmodavg)
prraw <- predictSE.gls(modelII, newdata = db2, se.fit = TRUE, print.matrix = FALSE)
db2$pred <-  prraw$fit 
db2$predseup <-  prraw$fit +1.96*prraw$se.fit
db2$predselo <-  prraw$fit -1.96*prraw$se.fit

kVH.obs    = db2$meanfem
kVH.obsH   = db2$seup
kVH.obsL   = ifelse(db2$selo<0,0,db2$selo)

k = 0.09
xv =0.85
xh =0.65
g = 1/4.5
wv = 1/2.5
m = 4*(0.031 + 95820*exp(db2$tm2m-50.4))

db2$r0.obs_M    = sqrt(0.101*kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(0.101*kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(0.101*kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))

Sys.setlocale("LC_TIME", "English")
db2$CommonDate <- as.Date((db2$julian-1),origin = "2014-1-1")
cd <- c(as.Date("2014-06-30"),db2$CommonDate)
db2$MonthDay <- as.Date(strftime(as.POSIXlt(db2$CommonDate,format="%Y-%m-%d"),origin="2014-01-01",format="%b-%d"),format="%b-%d")

dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))

png(paste(mypathfig,"Figure Supplementary 1.png",sep=""),width=900, height=900,units="px",bg="white")

ml <- rbind(c(1, 2), c(3, 4))
layout(ml)
#layout.show(4)

par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 3),  lwd = 2, pch = 5, axes = T, main = "Chikungunya Ro by observed data")
mtext("A", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic(R[0])), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$r0.obs_M, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$r0.obsL_M,y1=dbA$r0.obsH_M,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$r0.obs_M, cex =2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$r0.obsL_M,y1=dbB$r0.obsH_M,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=1,col="grey75",lty=2)
abline(h=2,col="grey75",lty=2)
abline(h=3,col="grey75",lty=2)
abline(h=4,col="grey75",lty=2)


kVH.obs    = db2$pred
kVH.obsH   = db2$predseup
kVH.obsL   = ifelse(db2$predselo>0,db2$predselo,0)

db2$r0.obs_M    = sqrt(0.101*kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(0.101*kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(0.101*kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))


dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 3),  lwd = 2, pch = 5, axes = T, main = "Chikungunya Ro by Model II")
mtext("B", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic(R[0])), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$r0.obs_M, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$r0.obsL_M,y1=dbA$r0.obsH_M,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$r0.obs_M, cex =2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$r0.obsL_M,y1=dbB$r0.obsH_M,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=1,col="grey75",lty=2)
abline(h=2,col="grey75",lty=2)
abline(h=3,col="grey75",lty=2)
abline(h=4,col="grey75",lty=2)


# zika
k = 0.09
xv =0.5#0.0665
xh =0.5
g = 1/5.8
wv = 1/10.5
m = 4*(0.031 + 95820*exp(db2$tm2m-50.4))

kVH.obs    = db2$meanfem
kVH.obsH   = db2$seup
kVH.obsL   = ifelse(db2$selo<0,0,db2$selo)

db2$r0.obs_M    = sqrt(0.101*kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(0.101*kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(0.101*kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))


dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 3),  lwd = 2, pch = 5, axes = T, main = "Zika Ro by observed data")
mtext("C", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic(R[0])), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$r0.obs_M, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$r0.obsL_M,y1=dbA$r0.obsH_M,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$r0.obs_M, cex =2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$r0.obsL_M,y1=dbB$r0.obsH_M,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=1,col="grey75",lty=2)
abline(h=2,col="grey75",lty=2)
abline(h=3,col="grey75",lty=2)
abline(h=4,col="grey75",lty=2)

kVH.obs    = db2$pred
kVH.obsH   = db2$predseup
kVH.obsL   = db2$predselo

db2$r0.obs_M    = sqrt(0.101*kVH.obs*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsH_M   = sqrt(0.101*kVH.obsH*k*xv*xh*wv/(g*m*(wv+m)))
db2$r0.obsL_M   = sqrt(0.101*kVH.obsL*k*xv*xh*wv/(g*m*(wv+m)))


dbA <- droplevels(subset(db2,Site == "Site-A"))
dbB <- droplevels(subset(db2,Site == "Site-B"))
par(cex.main = 1.5, mar = c(5, 5, 2, 1) + 0.1, mgp = c(1, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(cd, rep(-10, length(cd)), type = "p", ylab = "", xlab = " ", cex = 3, 
    ylim = c(0, 3),  lwd = 2, pch = 5, axes = T, main = "Zika Ro by Model II")
mtext("D", at=100,side = 2, line = 3, cex = 3, font = 2)
par(las = 0)
mtext(expression(italic(R[0])), side = 2, line = 3, cex = 2, font = 2)
points(dbA$CommonDate-0.5, dbA$r0.obs_M, cex = 2, pch = 19,col="dodgerblue4")
segments(x0=dbA$CommonDate-0.5,x1=dbA$CommonDate-0.5,y0=dbA$r0.obsL_M,y1=dbA$r0.obsH_M,col="dodgerblue4")
points(dbB$CommonDate+0.5, dbB$r0.obs_M, cex =2, pch = 19,col="#D55E00")
segments(x0=dbB$CommonDate+0.5,x1=dbB$CommonDate+0.5,y0=dbB$r0.obsL_M,y1=dbB$r0.obsH_M,col="#D55E00")
text(as.Date("2014-07-05"), 95, "Site-A", cex = 2, font = 1, adj = 0)
text(as.Date("2014-07-05"), 90, "Site-B", cex = 2, font = 1, adj = 0)
points(as.Date("2014-07-03"), 95, cex = 2, pch = 19,col="dodgerblue4")
points(as.Date("2014-07-03"), 90, cex = 2, pch = 19,col="#D55E00")
abline(h=1,col="grey75",lty=2)
abline(h=2,col="grey75",lty=2)
abline(h=3,col="grey75",lty=2)
abline(h=4,col="grey75",lty=2)
dev.off()