mu1 <- mu.tilde$muhat
mu1 <- t(matrix(mu1,m,n))

mu2 <- mu.hat$muhat
mu2 <- t(matrix(mu2,m,n))

Y <- t(matrix(y.sub,m,n))

pdf('ILI_MeanEst.pdf',width=10,height=7.5)

ids <- c(6,34,52,57) # Broward, Lake, Pinellas, Seminole
par(mfrow=c(2,2),mgp=c(2.4,1,0))

par(mar=c(3.5,3.5,1.5,0.5))
plot(1:365,mu2[,ids[1]],type='l',lty=1,lwd=1.5,xaxt='n',
     ylim=c(0,8e-5),main='Broward',xlab='Time',ylab='Incidence rate',
     cex=1.2,cex.lab=1.3,cex.axis=1.2,cex.main=1.3)
lines(1:365,mu1[,ids[1]],lty=2,lwd=2,cex=1.2)
points(1:365,Y[,ids[1]],cex=0.6)
axis(1,cex.axis=1.2,at=c(1+c(1,62,123,184,245,306, 367)),
     label=c('Jan','Mar','May','July','Sep','Nov','Jan'))
     
par(mar=c(3.5,3,1.5,1))
plot(1:365,mu2[,ids[2]],type='l',lty=1,lwd=1.5,xaxt='n',
     ylim=c(0,8e-5),main='Lake',xlab='Time',ylab='',
     cex=1.2,cex.lab=1.3,cex.axis=1.2,cex.main=1.3)
lines(1:365,mu1[,ids[2]],lty=2,lwd=2,cex=1.2)
points(1:365,Y[,ids[2]],cex=0.6)
axis(1,cex.axis=1.2,at=c(1+c(1,62,123,184,245,306, 367)),
     label=c('Jan','Mar','May','July','Sep','Nov','Jan'))
     
par(mar=c(3.5,3.5,1.5,0.5))
plot(1:365,mu2[,ids[3]],type='l',lty=1,lwd=1.5,xaxt='n',
     ylim=c(0,8e-5),main='Pinellas',xlab='Time',ylab='Incidence rate',
     cex=1.2,cex.lab=1.3,cex.axis=1.2,cex.main=1.3)
lines(1:365,mu1[,ids[3]],lty=2,lwd=2,cex=1.2)
points(1:365,Y[,ids[3]],cex=0.6)
axis(1,cex.axis=1.2,at=c(1+c(1,62,123,184,245,306, 367)),
     label=c('Jan','Mar','May','July','Sep','Nov','Jan'))
     
par(mar=c(3.5,3,1.5,1))
plot(1:365,mu2[,ids[4]],type='l',lty=1,lwd=1.5,xaxt='n',
     ylim=c(0,8e-5),main='Seminole',xlab='Time',ylab=' ',
     cex=1.2,cex.lab=1.3,cex.axis=1.2,cex.main=1.3)
lines(1:365,mu1[,ids[4]],lty=2,lwd=2,cex=1.2)
points(1:365,Y[,ids[4]],cex=0.6)
axis(1.2,at=c(1+c(1,62,123,184,245,306, 367)),
     label=c('Jan','Mar','May','July','Sep','Nov','Jan'))
     
graphics.off()