###########seasonal multiple linear regression model################
library(tidyverse)
library("showtext")
font_files()
showtext.auto(enable = TRUE)
library(fpp2)
library(ggpubr)
library(MASS)
library(qcc)
data<-readxl::read_excel('HFMD_st_dataset1.xlsx')

regions <- as.character(unique(data$province))
## Build a list that stores data for 31 cities 
CC <- lapply(regions, function(regions) data[data$province==regions,])
names(CC) <- regions

result<-data.frame()
for (i in 1:31) {
  TM<-as.data.frame(CC[i])
  colnames(TM) <- names(data)
  TM<-TM %>% group_by(year) %>% mutate(prop=cases/sum(cases))
  TM<-as.data.frame(TM)
  qcc<-qcc.overdispersion.test(TM$cases, type = 'poisson')
  
  #TM<-TM %>% mutate(lag.value1 = dplyr::lag(prop, n = 1, default = NA))
  #TM<-TM %>% mutate(lag.value1 = dplyr::lag(incidence, n = 1, default = NA))
  b=cos(2*pi*TM$time/12)
  c=sin(2*pi*TM$time/12)
  d=cos(4*pi*TM$time/12)
  e=sin(4*pi*TM$time/12)
  #fit <- lm(prop~b+c+d+e, data = TM)
  fit <- glm.nb(cases~b+c+d+e+offset(log(pop/100000))+year,data = TM)
  #acf(fit$residuals)
  b1<-fit$coefficients[2]
  c1<-fit$coefficients[3]
  d1<-fit$coefficients[4]
  e1<-fit$coefficients[5]
  annual<-sqrt(b1^2+c1^2)
  semi<-sqrt(d1^2+ e1^2)
  ANR<-annual/(annual+semi)
  SEMIR<-semi/(annual+semi)
  #R2<-summary(fit)$r.squared
  #adjR2<-summary(fit)$adj.r.squared
  if(b1>0&c1<0){
    peak<-12*(atan(c1/b1)+2*pi)/pi/2
  }else if(b1>0&c1>0){
    peak<-12*(atan(c1/b1))/pi/2
  }else if(b1<0&c1>0){
    peak<-12*(atan(c1/b1)+pi)/pi/2
  }else{
    peak<-12*(atan(c1/b1)+pi)/pi/2
  }
  
  if(annual>semi){
    number<-1
  }else{
    number<-2
  }
  
  if(d1>0&e1<0){
    semipeak<-12*(atan(e1/d1)+2*pi)/pi/2
  }else if(d1>0&e1>0){
    semipeak<-12*(atan(e1/d1))/pi/2
  }else if(d1<0&e1>0){
    semipeak<-12*(atan(e1/d1)+pi)/pi/2
  }else{
    semipeak<-12*(atan(e1/d1)+pi)/pi/2
  }
  
  
  
  result[i,1]<-TM[1,'province']
  result[i,2]<-peak
  result[i,3]<-annual
  result[i,4]<-semi
  result[i,5]<-ANR
  result[i,6]<-SEMIR
  #result[i,7]<-R2
  #result[i,8]<-adjR2
  result[i,7]<-number
  result[i,8]<-semipeak
  result[i,9]<-qcc[3]
  names(result)[1]<-'province'
  names(result)[2]<-'peak'
  names(result)[3]<-'annual'
  names(result)[4]<-'semi'
  names(result)[5]<-'ANR'
  names(result)[6]<-'SEMIR'
  #names(result)[7]<-'R2'
  #names(result)[8]<-'adjR2'
  names(result)[7]<-'number'
  names(result)[8]<-'semipeak'
  #TM<-TM[-1,]
  
  res<-data.frame(year=TM$year,month=TM$month,time=TM$time,residuals=fit$residuals)
  predict<-data.frame(year=TM$year,month=TM$month,time=TM$time,fit=fit$fitted.values,pop=TM$pop,TM$cases)
  names(predict)[6]<-'cases'
  #predict$fit[predict$fit<0]<-0
  predict$inc<-predict$fit/predict$pop*100000
  R2<-1 - (sum((predict$fit- predict$cases)^2) / sum((predict$cases - mean(predict$cases))^2))
  #r2<-cor.test(TM$cases,predict$fit)
  r2<-round(R2,2)
  result[i,9]<-r2
  names(result)[9]<-'R2'
  #r2<-round(R2,2)  
  fitness<-ggplot(TM,aes(time,incidence))+
    geom_line()+
    geom_line(data=predict,aes(time,inc),color='red')+
    theme_bw()+
    theme(panel.grid = element_blank(),
          plot.title = element_text(hjust = 0.5,size = rel(1.8)))+
    labs(x='Month',y='Monthly HFMD incidence',
         title = TM[1,'province'])+
    scale_x_continuous(limits = c(0,120),breaks = seq(0,120,12))+
    scale_y_continuous(limits = c(0,150),breaks = seq(0,150,30))+
    theme(axis.title = element_text(size = rel(1.4)),
          axis.text = element_text(size = rel(1.4)))+
    annotate('text',x=108,y=136,
             label=bquote ('r' ~'='~ .(r2)),size=5)
  
  resp<-ggplot(res,aes(time,residuals))+
    geom_point(shape=1)+
    theme_bw()+
    theme(panel.grid = element_blank(),
          plot.title = element_text(hjust = 0.5,size = rel(1.8)))+
    labs(x='Month',y='Residuals',
         title = TM[1,'province'])+
    scale_x_continuous(limits = c(0,120),breaks = seq(0,120,12))+
    scale_y_continuous(limits = c(-10,10),breaks = seq(-10,10,5))+
    theme(axis.title = element_text(size = rel(1.8)),
          axis.text = element_text(size = rel(1.8))
    )
  #hz<-ggarrange(fitness,resp,ncol=1,nrow = 2,hjust = 0)
  file_name = paste("fit", TM[1,'province'], ".pdf", sep="")
  ggsave(fitness,file=file_name,width = 6,height = 4)
}

meancase<-data %>% group_by(province,year) %>% summarise(case=sum(incidence)) %>% mutate(mean_ind=mean(case))
meancase<-meancase %>% filter(year==2010)
meancase<-meancase[,c(1,4)]
result<-left_join(result,id)
result<-left_join(result,X3region)
result<-left_join(result,meancase)

fit<-summary(lm(peak~lat,data=result))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
#####fig 3A#####
ggplot(result,aes(peak,lat))+
  geom_point(aes(color=region,size=mean_ind),shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_x_continuous(limits = c(5,8))+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Peak time')+
  annotate('text',x=7.6,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=7.6,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='peak_all.pdf',width = 5,height = 5)

fit<-summary(lm(annual~lat,data=result))

r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
#####fig 3B#####
ggplot(result,aes(annual*100,lat))+
  geom_point(aes(color=region,size=mean_ind),shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0,250),breaks = seq(0,250,50))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Annual amplitude (%)')+
  annotate('text',x=226,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=226,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='annual_all.pdf',width = 5,height = 5)


fit<-summary(lm(semi~lat,data=result))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
#####fig 3C#####
ggplot(result,aes(semi*100,lat))+
  geom_point(aes(color=region,size=mean_ind),shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(40,101),breaks = seq(40,100,10))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Semi-annual amplitude (%)')+
  annotate('text',x=90,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=90,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='semi_all.pdf',width = 5,height = 5)

fit<-summary(lm(SEMIR~lat,data=result))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
#####fig 3D#####
ggplot(result,aes(SEMIR,lat))+
  geom_point(aes(color=region,size=mean_ind),shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0.1,0.74),breaks = seq(0.1,0.7,0.2))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Ratio semiannual component')+
  annotate('text',x=0.66,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=0.66,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='SEMIR_all.pdf',width = 5,height = 5)


summary(lm(peak~lat,data=result))
summary(lm(annual~lat,data=result))
summary(lm(semi~lat,data=result))
summary(lm(SEMIR~lat,data=result))

##############East###############
b<-filter(result,region==1)
summary(lm(peak~lat,data=b))
summary(lm(annual~lat,data=b))
summary(lm(semi~lat,data=b))
summary(lm(SEMIR~lat,data=b))

fit<-summary(lm(peak~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(peak,lat))+
  geom_point(aes(size=mean_ind),color='red',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_x_continuous(limits = c(5,8))+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Peak time')+
  annotate('text',x=7.6,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=7.6,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='peak_eastern.pdf',width = 5,height = 5)

fit<-summary(lm(annual~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(annual*100,lat))+
  geom_point(aes(size=mean_ind),color='red',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0,250),breaks = seq(0,250,50))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Annual amplitude (%)')+
  annotate('text',x=226,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=226,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='annual_eastern.pdf',width = 5,height = 5)

fit<-summary(lm(semi~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(semi*100,lat))+
  geom_point(aes(size=mean_ind),color='red',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(40,101),breaks = seq(40,101,10))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Semi-annual amplitude (%)')+
  annotate('text',x=90,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=90,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='semi_eastern.pdf',width = 5,height = 5)

fit<-summary(lm(SEMIR~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(SEMIR,lat))+
  geom_point(aes(size=mean_ind),color='red',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0.1,0.74),breaks = seq(0.1,0.7,0.2))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Ratio semiannual component')+
  annotate('text',x=0.66,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=0.66,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='SEMIR_eastern.pdf',width = 5,height = 5)

##############central###############
b<-filter(result,region==2)
summary(lm(peak~lat,data=b))
summary(lm(annual~lat,data=b))
summary(lm(semi~lat,data=b))
summary(lm(SEMIR~lat,data=b))

fit<-summary(lm(peak~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(peak,lat))+
  geom_point(aes(size=mean_ind),color='green',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_x_continuous(limits = c(5,8))+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Peak time')+
  annotate('text',x=7.6,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=7.6,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='peak_central.pdf',width = 5,height = 5)

fit<-summary(lm(annual~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(annual*100,lat))+
  geom_point(aes(size=mean_ind),color='green',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0,250),breaks = seq(0,250,50))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Annual amplitude (%)')+
  annotate('text',x=226,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=226,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='annual_central.pdf',width = 5,height = 5)

fit<-summary(lm(semi~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(semi*100,lat))+
  geom_point(aes(size=mean_ind),color='green',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(40,101),breaks = seq(40,101,10))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Semi-annual amplitude (%)')+
  annotate('text',x=90,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=90,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='semi_central.pdf',width = 5,height = 5)

fit<-summary(lm(SEMIR~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(SEMIR,lat))+
  geom_point(aes(size=mean_ind),color='green',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0.1,0.74),breaks = seq(0.1,0.7,0.2))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Ratio semiannual component')+
  annotate('text',x=0.66,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=0.66,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='SEMIR_central.pdf',width = 5,height = 5)


##############west###############
b<-filter(result,region==3)
summary(lm(peak~lat,data=b))
summary(lm(annual~lat,data=b))
summary(lm(semi~lat,data=b))
summary(lm(SEMIR~lat,data=b))

fit<-summary(lm(peak~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(peak,lat))+
  geom_point(aes(size=mean_ind),color='blue',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_x_continuous(limits = c(5,8),breaks = seq(5,8,1))+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Peak time')+
  annotate('text',x=7.6,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=7.6,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='peak_western.pdf',width = 5,height = 5)

fit<-summary(lm(annual~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(annual*100,lat))+
  geom_point(aes(size=mean_ind),color='blue',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0,250),breaks = seq(0,250,50))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Annual amplitude (%)')+
  annotate('text',x=226,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=226,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='annual_western.pdf',width = 5,height = 5)

fit<-summary(lm(semi~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(semi*100,lat))+
  geom_point(aes(size=mean_ind),color='blue',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(40,101),breaks = seq(40,101,10))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Semi-annual amplitude (%)')+
  annotate('text',x=90,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=90,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='semi_western.pdf',width = 5,height = 5)

fit<-summary(lm(SEMIR~lat,data=b))
r2<-round(fit$r.squared,2)
p<-round(fit$coefficients['lat','Pr(>|t|)'],3)
ggplot(b,aes(SEMIR,lat))+
  geom_point(aes(size=mean_ind),color='blue',shape=1)+
  theme_bw()+
  stat_smooth(method = lm,se=F,color='black')+
  theme(plot.title =element_text(hjust = 0.5,size=16),
        axis.text = element_text(size = rel(1.5)),
        axis.title = element_text(size = rel(1.5)),
        legend.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none')+
  scale_y_continuous(limits = c(20,46),breaks = seq(20,45,5))+
  scale_x_continuous(limits = c(0.1,0.74),breaks = seq(0.1,0.7,0.2))+
  scale_size_continuous(range = c(2,16))+
  labs(y='Latitude',x='Ratio semiannual component')+
  annotate('text',x=0.66,y=22,
           label=bquote (R^2 ~'='~ .(r2)),size=5)+
  annotate('text',x=0.66,y=20,
           label=bquote (italic('P') ~'='~ .(p)),size=5)
ggsave(file='SEMIR_western.pdf',width = 5,height = 5)


range(result$R2)
median(result$R2)






############sp+dlnm################################################################
###loading of the package###
library(INLA)

packages <- c("data.table", "tidyverse", "sf", "sp", "spdep",
              "dlnm", "tsModel", "hydroGOF","RColorBrewer", 
              "geofacet", "ggpubr", "ggthemes","splines")
lapply(packages, library, character.only = TRUE)
####Map file import and adjustment####
map <- sf::st_read('newshp.shp')
names(map)[10]<-'province'
map$ID<-seq(1,31)
china <- poly2nb(map)
nbs<-poly2nb(map, queen = T, row.names = map$ID)
map[11,10]<-'Inner Mongolia'
map[4,10]<-'Tibet'
####Data import adjustments####
data<-readxl::read_excel('HFMD_st_dataset1.xls')
g.file <- "china.graph"
if (!file.exists(g.file)) nb2INLA(g.file, nbs)
data<-left_join(data,map[,c(10:12)])
data<-data[with(data,order(ID,year,month)),]
data$tcases<-as.numeric(data$tcases)
data$tincidence<-as.numeric(data$tincidence)
data$PM2_5<-log(data$PM2_5)
data$SO2<-log(data$SO2)
data$O3_8h<-log(data$O3_8h)
###Determine the max lag###
nlag<-18

# creating lagged variables
# define the monthly mean climate variable lag term matrix
lag_mt <- tsModel::Lag(data$mean_tem, group = data$ID, k = 0:nlag)
lag_rhu <- tsModel::Lag(data$rhu, group = data$ID, k = 0:nlag)
lag_win <- tsModel::Lag(data$win, group = data$ID, k = 0:nlag)
lag_pre <- tsModel::Lag(data$pre, group = data$ID, k = 0:nlag)
lag_ndvi <- tsModel::Lag(data$ndvi, group = data$ID, k = 0:nlag)
lag_ssd <- tsModel::Lag(data$ssd, group = data$ID, k = 0:nlag)
lag_prs <- tsModel::Lag(data$prs, group = data$ID, k = 0:nlag)
lag_SO2 <- tsModel::Lag(data$SO2, group = data$ID, k = 0:nlag)
lag_NO2 <- tsModel::Lag(data$NO2, group = data$ID, k = 0:nlag)
lag_PM10 <- tsModel::Lag(data$PM10, group = data$ID, k = 0:nlag)
lag_CO <- tsModel::Lag(data$CO, group = data$ID, k = 0:nlag)
lag_O3_8h <- tsModel::Lag(data$O3_8h, group = data$ID, k = 0:nlag)
lag_PM2_5 <- tsModel::Lag(data$PM2_5, group = data$ID, k = 0:nlag)
# cross basis matrix
# Set the lag section
lagknot <- equalknots(0:nlag, 2)
varknot<-equalknots(data$mean_tem, 2)
var <- lag_mt
basis_mt <- crossbasis(var,
                       argvar = list(fun = "ns", knots = varknot),
                       arglag = list(fun = "ns", knots = lagknot))


var <- lag_rhu
basis_rhu <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$rhu, 2)),
                        arglag = list(fun = "ns", knots = lagknot))

var <- lag_win
basis_win <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$win, 2)),
                        arglag = list(fun = "ns", knots = lagknot))

var <- lag_ndvi
basis_ndvi <- crossbasis(var,
                         argvar = list(fun = "ns", knots = equalknots(data$ndvi, 2)),
                         arglag = list(fun = "ns", knots = lagknot))
var <- lag_ssd
basis_ssd <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$ssd, 2)),
                        arglag = list(fun = "ns", knots = lagknot))
var <- lag_prs
basis_prs <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$prs, 2)),
                        arglag = list(fun = "ns", knots = lagknot))
var <- lag_NO2
basis_NO2 <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$NO2, 2)),
                        arglag = list(fun = "ns", knots = lagknot))

var <- lag_CO
basis_CO<- crossbasis(var,
                      argvar = list(fun = "ns", knots = equalknots(data$CO, 2)),
                      arglag = list(fun = "ns", knots = lagknot))

var <- lag_SO2
basis_SO2<- crossbasis(var,
                       argvar = list(fun = "ns", knots = equalknots(data$SO2, 2)),
                       arglag = list(fun = "ns", knots = lagknot))
var <- lag_O3_8h
basis_O3_8h <- crossbasis(var,
                          argvar = list(fun = "ns", knots = equalknots(data$O3_8h, 2)),
                          arglag = list(fun = "ns", knots = lagknot))
var <- lag_PM10
basis_PM10 <- crossbasis(var,
                         argvar = list(fun = "ns", knots = equalknots(data$PM10, 2)),
                         arglag = list(fun = "ns", knots = lagknot))

var <- lag_pre
basis_pre <- crossbasis(var,
                        argvar = list(fun = "ns", knots = equalknots(data$pre, 2)),
                        arglag = list(fun = "ns", knots = lagknot))

var <- lag_PM2_5
basis_PM2_5 <- crossbasis(var,
                          argvar = list(fun = "ns", knots = equalknots(data$PM2_5, 2)),
                          arglag = list(fun = "ns", knots = lagknot))

# assign unique column names to cross-basis matrix for inla() model
# note: not necessary for glm(), gam() or glm.nb() models
colnames(basis_mt) = paste0("basis_mt.", colnames(basis_mt))
colnames(basis_rhu) = paste0("basis_rhu.", colnames(basis_rhu))
colnames(basis_win) = paste0("basis_win.", colnames(basis_win))
colnames(basis_ndvi) = paste0("basis_ndvi.", colnames(basis_ndvi))
colnames(basis_ssd) = paste0("basis_ssd.", colnames(basis_ssd))
colnames(basis_prs) = paste0("basis_prs.", colnames(basis_prs))
colnames(basis_SO2) = paste0("basis_SO2.", colnames(basis_SO2))
colnames(basis_NO2) = paste0("basis_NO2.", colnames(basis_NO2))
colnames(basis_CO) = paste0("basis_CO.", colnames(basis_CO))
colnames(basis_PM10) = paste0("basis_rhu.", colnames(basis_PM10))
colnames(basis_O3_8h) = paste0("basis_O3_8h.", colnames(basis_O3_8h))
colnames(basis_pre) = paste0("basis_pre.", colnames(basis_pre))
colnames(basis_PM2_5) = paste0("basis_PM2_5.", colnames(basis_PM2_5))


#########################################################################################################################################
# set up data and priors for INLA model
# set data for models
Y  <- data$tcases # response variable
E  <- data$popu/10^5 # model offset so that response is equivalent to an incidence rate per 100,000 people
# random variable
T1 <- data$month # for random effect to account for annual cycle (seasonality)
data$year_index <- data$year-2003
T2 <- data$year_index # for random effect to account for inter-annual variability
S <- data$ID # for country interaction with month and country spatial random effect
# Other covariables
ndvi<-data$ndvi
#healthper<-data$healthper
holiday<-data$holiday
birth<-data$birth
beds<-data$beds
mt<-data$mean_tem
rhu<-data$rhu
win<-data$win
ssd<-data$ssd
lncar<-data$lncar
lnpopden<-data$lnpopden
lnrgdp<-data$lnrgdp
prs<-data$prs
city<-data$city
PM10<-data$PM10
SO2<-data$SO2
NO2<-data$NO2
CO<-data$CO
O3_8h<-data$O3_8h
beds<-data$beds
healthper<-data$healthper
sixthfive<-data$sixtyfive
pre<-data$pre
lnrgdp<-data$lnrgdp
PM2_5<-data$PM2_5
region<-data$region

# create dataframe for model testing
df <- data.frame(Y, E, T1, T2, S, ndvi,holiday,birth,beds,mt,rhu,win,ssd,lncar,lnpopden,lnrgdp,prs,NO2,SO2,CO,PM10,O3_8h,pre,sixthfive,beds,healthper,PM2_5)

# define priors
precision.prior <- list(prec = list(prior = "pc.prec", param = c(0.5, 0.01)))

# inla model function

# include formula and set defaults for data, family (to allow other prob dist models e.g. Poisson) and config (to allow for sampling)
mymodel <- function(formula, data = df, family = "nbinomial", config = FALSE)
  
{
  model <- inla(formula = formula, data = data, family = family, offset = log(E),
                control.inla = list(strategy = 'adaptive',cmin=0), 
                control.compute = list(dic = TRUE, config = config, 
                                       cpo = TRUE, return.marginals = FALSE,
                                       waic=TRUE),
                control.fixed = list(correlation.matrix = TRUE,
                                     mean.intercept=0, 
                                     prec.intercept=1, # tight precision
                                     mean=0, 
                                     prec=1),
                control.predictor = list(link = 1, compute = TRUE), 
                verbose = FALSE)
  model <- inla.rerun(model)
  return(model)
}




# fixed effects priors
control.fixed1 = list(mean.intercept=0, # prior mean for intercept
                      prec.intercept=0.05, # prior precision for intercept (wide sd bc. n binomial; intercept likely v small)
                      mean=0, # prior mean for fixed effects
                      prec=1)  # prior precision for (scaled) fixed effects

# hyperpriors for iid, ar1, bym2
hyper.iid = list(theta = list(prior="pc.prec", param=c(1, 0.01)))

# for bym or bym2
hyper.bym = list(theta1 = list(prior="pc.prec", param=c(1, 0.01)),
                 theta2 = list(prior="pc.prec", param=c(1, 0.01)))
hyper.bym2 = list(theta1 = list(prior="pc.prec", param=c(1, 0.01)),
                  theta2 = list(prior="pc", param=c(0.5, 0.5)))

# ar1
hyper.ar1 = list(theta1 = list(prior='pc.prec', param=c(1, 0.01)), 
                 rho = list(prior='pc.cor0', param = c(0.75, 0.5)))

# random walk hyperparameters on smoothing
hyper1.rw = list(prec = list(prior='pc.prec', param=c(0.1, 0.01))) # strictest; sd constrained to be low
hyper2.rw = list(prec = list(prior='pc.prec', param=c(0.5, 0.01))) # medium
hyper3.rw = list(prec = list(prior='pc.prec', param=c(1, 0.01))) # weakest; sd can be quite wide 

####The formula for the basic model####
baseformula1 <- Y ~ 1 + f(T1, replicate = S, model = "rw1", cyclic = TRUE, constr = TRUE,
                          scale.model = TRUE,  hyper = hyper2.rw) +
  f(S, model = "bym2", graph = "china.graph",replicate = T2,constr = TRUE,
    scale.model = TRUE, hyper = hyper.bym2)+lnrgdp+lnpopden+healthper+city+sixthfive+beds+mt+rhu+win+pre+ssd+SO2+O3_8h+PM2_5

####Lagging model formulas####
finalformula <- Y ~ 1 + f(T1, replicate = S, model = "rw1", cyclic = TRUE, constr = TRUE,
                          scale.model = TRUE,  hyper = hyper2.rw) +
  f(S, model = "bym2", graph = "china.graph",replicate = T2,constr = TRUE,
    scale.model = TRUE, hyper = hyper.bym2)+basis_PM2_5+basis_win+basis_pre+basis_ssd+
  basis_mt+basis_rhu+lnrgdp+lnpopden+healthper+city+sixthfive+beds+basis_SO2+basis_O3_8h


####The basic regression model####
base<-mymodel(Y ~ 1+lnrgdp+lnpopden+healthper+city+sixthfive+beds+mt+rhu+win+pre+ssd+SO2+O3_8h+PM2_5)
####dlnm####
gam<-mymodel(Y ~ 1+basis_PM2_5+basis_win+basis_pre+basis_ssd+
               basis_mt+basis_rhu+lnrgdp+lnpopden+healthper+city+sixthfive+beds+basis_SO2+basis_O3_8h)
####sp####
sp<-mymodel(baseformula1, df)
####sp+dlnm####
spdlnm<-mymodel(finalformula, df)
####result####
summary(base)
summary(gam)
summary(sp)
summary(spdlnm)
###dic####
base$dic$dic
gam$dic$dic
sp$dic$dic
spdlnm$dic$dic
#####other result#####
rr<-exp(sp$summary.fixed)
round(-mean(log(base$cpo$cpo), na.rm = T), 3)
round(-mean(log(gam$cpo$cpo), na.rm = T), 3)
round(-mean(log(sp$cpo$cpo), na.rm = T), 3)
round(-mean(log(spdlnm$cpo$cpo), na.rm = T), 3)
relation1 <- lm(data$cases~spdlnm1$summary.fitted.values$mean)
relation2 <- lm(data$cases~spdlnm2$summary.fitted.values$mean)
relation3 <- lm(data$cases~spdlnm3$summary.fitted.values$mean)
summary(relation1)
summary(relation2)
summary(relation3)

postpred<-data.frame(data[,c(1:3,29)],spdlnm$summary.fitted.values$mean)
names(postpred)[5]<-'pred'

postpred$pred<-round(postpred$pred,0)
postpred$indicence<-postpred$pred/postpred$pop*100000

writexl::write_xlsx(postpred,'postpred.xlsx')

###################Random effects fs6-7##############
model<-spdlnm3
month_effects <- data.table(cbind(rep(unique(data$ID), each = 12),
                                  model$summary.random$T1))
names(month_effects)[1:2] <- c("ID", "Month")

month_effects<-left_join(month_effects,data[!duplicated(data[,c(1,33)]),][,c(1,33)])
ggplot(month_effects,aes(Month,mean))+
  geom_line()+
  geom_ribbon(aes(ymax=`0.975quant`,ymin=`0.025quant`),fill='grey',alpha=0.4)+
  facet_wrap(~month_effects$province, ncol = 5)+
  labs(y='Contribution to log(HFMD incidence)')+
  scale_x_continuous(breaks = seq(1,12,1))+
  theme_bw()+theme(panel.grid = element_blank())
ggsave("Fig S8_month_effect.pdf", height = 20, width = 30, units = "cm")  

write.csv(month_effects,file = "Random effects.csv",row.names = F)
#########################################################################################################################################
## make maps of spatial random effects per year
# extract posterior mean estimates for combined unstructured and structured random effects
n_country <- length(unique(data$ID))
nyear <- length(unique(data$year))

space <- data.table(model$summary.random$S)
space$year <- rep(min(data$year):max(data$year), each = 2*n_country)
space$re <- rep(c(rep(1,n_country),rep(2,n_country)),nyear)
space <- space[space$re == 2,]
space$ID <- rep(unique(data$ID), nyear)
mn <-min(space$mean)
mx <-max(space$mean)

write.csv(space,file = "Fig S9_spatial random effects per yeari.csv",row.names = F)

# Add the map geometry to the space dataframe
province<-readxl::read_excel('D:/DESKTOP/毕业论文/data/data/province_wushi.xlsx')
mapping <- sf::st_read('F:/region.shp')
#mapping <-left_join(province[,c('省份','province')],mapping)
names(mapping)[10]<-'province'
mapping[13,10]<-'Inner Mongolia'
mapping[5,10]<-'Tibet'
mapping<-left_join(mapping,data[!duplicated(data[,c(1,32)]),][,c(1,32)])

space <- left_join(map, space)

space_effects <- ggplot() + 
  geom_sf(data = space, aes(fill = mean), lwd = 0, color = NA) +
  scale_fill_distiller(palette = "RdBu", direction = -1, 
                       limits = c(min(mn,-mx),max(mx,-mn))) +
  labs(fill = "Contribution to \n log(HFMD incidence)") +
  theme_void() +
  facet_wrap(~space$year, ncol = 5)

ggsave("Fig S9_year_spatial_effect.pdf", height = 20, width = 30, units = "cm")

##############figure 7A-D##########

coef <- model$summary.fixed$mean
vcov <- model$misc$lincomb.derived.covariance.matrix
# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_mt", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_mt, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$mean_tem), 0)) 

par(mar = c(5, 4, 2, 4))
plot(predt,"overall",col=2,xlim=c(-25,35),ylim=c(-1.,4),axes=F,ann=F,lwd=2.5,cex.axis=5)
#lines(predt,"overall",ci="area",col=2,lwd=2.5)
axis(1,at=10*(-3:4),cex.axis=1.2)
axis(2,at=c(0:4*2),cex.axis=1.2)
title(xlab= expression(paste("Temperature (",degree,"C)")),cex.lab=1.2)
mtext("Relative risk",side = 2,at=4, line =2.5,cex=1.2)
par(new=T)
hist(data$mean_tem,xlim=c(-25,35),ylim=c(0,1600),axes=F,ann=F,col="wheat",breaks=20)
axis(4,at=(0:2)*200,cex.axis=1.2)
mtext("Frequency",side =4,at=200, line =2.5,cex=1.2,las=3)

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_pm", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_pm, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$pm), 0)) 

range(data$pm)
par(mar = c(5, 4, 2, 4))
plot(predt,"overall",col=2,xlim=c(0,90),ylim=c(0,2),axes=F,ann=F,lwd=2.5,cex.axis=5)
axis(1,at=10*(0:9),cex.axis=1.2)
axis(2,at=c(0:4*1),cex.axis=1.2)
title(xlab= "pm",cex.lab=1.2)
mtext("Relative risk",side = 2,at=2, line =2.5,cex=1.2)
par(new=T)
hist(data$pm,xlim=c(0,90),ylim=c(0,2000),axes=F,ann=F,col="wheat",breaks=20)
axis(4,at=(0:3)*200,cex.axis=1.2)

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_win", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_win, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$win), 0)) 

range(data$win)
par(mar = c(5, 4, 2, 4))
plot(predt,"overall",col=2,xlim=c(0,5),ylim=c(-1.,2),axes=F,ann=F,lwd=2.5,cex.axis=5)
axis(1,at=1*(0:9),cex.axis=1.2)
axis(2,at=c(0:4*1),cex.axis=1.2)
title(xlab= "Wind speed (m/s)",cex.lab=1.2)
mtext("Relative risk",side = 2,at=1, line =2.5,cex=1.2)
par(new=T)
hist(data$win,xlim=c(0,5),ylim=c(0,2000),axes=F,ann=F,col="wheat",breaks=20)
axis(4,at=(0:3)*200,cex.axis=1.2)

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_ndvi", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_ndvi, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$ndvi), 0)) 

range(data$ndvi)
par(mar = c(5, 4, 2, 4))
plot(predt,"overall",col=2,xlim=c(0,1),ylim=c(-1.,2),axes=F,ann=F,lwd=2.5,cex.axis=5)
axis(1,at=1*(0:9),cex.axis=1.2)
axis(2,at=c(0:4*1),cex.axis=1.2)
title(xlab= "NVDI",cex.lab=1.2)
mtext("Relative risk",side = 2,at=1, line =2.5,cex=1.2)
par(new=T)
hist(data$win,xlim=c(0,5),ylim=c(0,2000),axes=F,ann=F,col="wheat",breaks=20)
axis(4,at=(0:3)*200,cex.axis=1.2)

########figure 7I-L#######


coef <- model$summary.fixed$mean
vcov <- model$misc$lincomb.derived.covariance.matrix


# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_pm", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_pm, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$pm), 0)) 
y <- predt$predvar
x <- seq(0, nlag, 0.25)
z <- t(predt$matRRfit)

pal <- rev(brewer.pal(11, "PRGn"))
levels <- pretty(z, 20)
col1 <- colorRampPalette(pal[1:6])
col2 <- colorRampPalette(pal[6:11])
cols <- c(col1(sum(levels < 1)), col2(sum(levels > 1)))

filled.contour(x,y,z,
               xlab = "Lag", ylab = expression(paste("pm (",degree,"C)")), main = "",
               col = cols,levels = levels,
               plot.axes = { axis(1, at = 0:nlag, c(0:nlag)) 
                 axis(2)})

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_rhu", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_rhu, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$rhu), 0)) 
y <- predt$predvar
x <- seq(0, nlag, 0.25)
z <- t(predt$matRRfit)

pal <- rev(brewer.pal(11, "PRGn"))
levels <- pretty(z, 20)
col1 <- colorRampPalette(pal[1:6])
col2 <- colorRampPalette(pal[6:11])
cols <- c(col1(sum(levels < 1)), col2(sum(levels > 1)))

filled.contour(x,y,z,
               xlab = "Lag", ylab = "Relative humidity (%)", main = "",
               col = cols,levels = levels,
               plot.axes = { axis(1, at = 0:nlag, c(0:nlag)) 
                 axis(2)})


# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_ndvi", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_ndvi, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$ndvi), 0)) 
y <- predt$predvar
x <- seq(0, nlag, 0.25)
z <- t(predt$matRRfit)

pal <- rev(brewer.pal(11, "PRGn"))
levels <- pretty(z, 20)
col1 <- colorRampPalette(pal[1:6])
col2 <- colorRampPalette(pal[6:11])
cols <- c(col1(sum(levels < 1)), col2(sum(levels > 1)))

filled.contour(x,y,z,
               xlab = "Lag", ylab = "NDVI", main = "",
               col = cols,levels = levels,
               plot.axes = { axis(1, at = 0:nlag, c(0:nlag)) 
                 axis(2)})


# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_win", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_win, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$win), 0)) 
y <- predt$predvar
x <- seq(0, nlag, 0.25)
z <- t(predt$matRRfit)

pal <- rev(brewer.pal(11, "PRGn"))
levels <- pretty(z, 20)
col1 <- colorRampPalette(pal[1:6])
col2 <- colorRampPalette(pal[6:11])
cols <- c(col1(sum(levels < 1)), col2(sum(levels > 1)))

filled.contour(x,y,z,
               xlab = "Lag", ylab = "Wind speed (m/s)", main = "",
               col = cols,levels = levels,
               plot.axes = { axis(1, at = 0:nlag, c(0:nlag)) 
                 axis(2)})





######figure 7E-H#######

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_mt", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_mt, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$mean_tem), 0)) 
# get exposures values
vars <- predt$predvar

# obtain relative risk (RR) fit and upper and lower confidence limits for all exposure variables
rr <- predt$matRRfit
rr.lci <- predt$matRRlow
rr.uci <- predt$matRRhigh

# set relative risk range 
r1 <- min(range(rr, rr.lci, rr.uci))
r2 <- max(range(rr, rr.lci, rr.uci))

summary(data$mean_tem)
# get selected exposure variable positions
mn <- which(round(vars, 2) == round(quantile(data$mean_tem,0.05),0))
mx <- which(round(vars, 2) == round(quantile(data$mean_tem,0.95),0))

# define colours
col1 <- brewer.pal(11, "RdBu")[9]
tcol1 <- do.call(rgb, c(as.list(col2rgb(col1)), alpha = 255/4, max = 255))

col2 <- brewer.pal(11, "RdBu")[3]
tcol2 <- do.call(rgb, c(as.list(col2rgb(col2)), alpha = 255/4, max = 255))

col3 <- brewer.pal(11, "RdBu")[1]
tcol3 <- do.call(rgb, c(as.list(col2rgb(col3)), alpha = 255/4, max = 255))

# define x values (lag, by lag)
lagbylag <- seq(0, nlag, 0.25)

# cool
plot(lagbylag, rr[mn,], col = col1, type = "l", lwd = 1, 
     xlab = "Lag", ylab = "Relative risk", main = "", 
     ylim = range(r1, r2*1.11), frame.plot = T, axes = F)
axis(1, at = 0:nlag, labels = 0:nlag)
axis(2)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mn,], rev(rr.uci[mn,]))
polygon(xx, yy, col = tcol1, border = tcol1)
# warm
lines(lagbylag, rr[mx,], col = col2, lwd = 1)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mx,], rev(rr.uci[mx,]))
polygon(xx, yy, col = tcol2, border = tcol2)
abline(h = 1, lty = 3)

legend("topleft",
       legend = c(expression(paste("Tem =",' -7 ',degree,"C")),
                  expression(paste("Tem =",' 28 ',degree,"C"))),
       col = c(col1, col2), 
       lwd = 2, lty = 1, bty = "n", 
       y.intersp = 1.5, horiz = F)

#mtext(side = 2, at = r2*1.25, text = "b", las = 2, cex = 1.2, line = 2)


# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_rhu", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_rhu, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$rhu), 0)) 
# get exposures values
vars <- predt$predvar

# obtain relative risk (RR) fit and upper and lower confidence limits for all exposure variables
rr <- predt$matRRfit
rr.lci <- predt$matRRlow
rr.uci <- predt$matRRhigh

# set relative risk range 
r1 <- min(range(rr, rr.lci, rr.uci))
r2 <- max(range(rr, rr.lci, rr.uci))

summary(data$rhu)
# get selected exposure variable positions
mn <- which(round(vars, 2) == round(quantile(data$rhu,0.05),0))
mx <- which(round(vars, 2) == round(quantile(data$rhu,0.95),0))

# define colours
col1 <- brewer.pal(11, "RdBu")[9]
tcol1 <- do.call(rgb, c(as.list(col2rgb(col1)), alpha = 255/4, max = 255))

col2 <- brewer.pal(11, "RdBu")[3]
tcol2 <- do.call(rgb, c(as.list(col2rgb(col2)), alpha = 255/4, max = 255))


# define x values (lag, by lag)
lagbylag <- seq(0, nlag, 0.25)

# cool
plot(lagbylag, rr[mn,], col = col1, type = "l", lwd = 1, 
     xlab = "Lag", ylab = "Relative risk", main = "", 
     ylim = range(r1, r2*1.11), frame.plot = T, axes = F)
axis(1, at = 0:nlag, labels = 0:nlag)
axis(2)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mn,], rev(rr.uci[mn,]))
polygon(xx, yy, col = tcol1, border = tcol1)
# warm
lines(lagbylag, rr[mx,], col = col2, lwd = 1)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mx,], rev(rr.uci[mx,]))
polygon(xx, yy, col = tcol2, border = tcol2)
abline(h = 1, lty = 3)


legend("topleft",
       legend = c(paste0("Rhu = ",vars[mn],"%"),
                  paste0("Rhu = ", vars[mx],"%")),
       col = c(col1, col2), 
       lwd = 2, lty = 1, bty = "n", 
       y.intersp = 1.5, horiz = F)

# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_ndvi", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_ndvi, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$ndvi), 0)) 
# get exposures values
vars <- predt$predvar

# obtain relative risk (RR) fit and upper and lower confidence limits for all exposure variables
rr <- predt$matRRfit
rr.lci <- predt$matRRlow
rr.uci <- predt$matRRhigh

# set relative risk range 
r1 <- min(range(rr, rr.lci, rr.uci))
r2 <- max(range(rr, rr.lci, rr.uci))

summary(data$ndvi)
round(quantile(data$ndvi,0.025),1)
# get selected exposure variable positions
mn <- which(round(vars, 2) == round(quantile(data$ndvi,0.05),1))
mx <- which(round(vars, 2) == round(quantile(data$ndvi,0.95),1))


# define colours
col1 <- brewer.pal(11, "RdBu")[9]
tcol1 <- do.call(rgb, c(as.list(col2rgb(col1)), alpha = 255/4, max = 255))

col2 <- brewer.pal(11, "RdBu")[3]
tcol2 <- do.call(rgb, c(as.list(col2rgb(col2)), alpha = 255/4, max = 255))


# define x values (lag, by lag)
lagbylag <- seq(0, nlag, 0.25)

# cool
plot(lagbylag, rr[mn,], col = col1, type = "l", lwd = 1, 
     xlab = "Lag", ylab = "Relative risk", main = "", 
     ylim = range(r1, r2*1.11), frame.plot = T, axes = F)
axis(1, at = 0:nlag, labels = 0:nlag)
axis(2)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mn,], rev(rr.uci[mn,]))
polygon(xx, yy, col = tcol1, border = tcol1)
# warm
lines(lagbylag, rr[mx,], col = col2, lwd = 1)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mx,], rev(rr.uci[mx,]))
polygon(xx, yy, col = tcol2, border = tcol2)
abline(h = 1, lty = 3)



legend("topleft",
       legend = c(paste0("NDVI = ",vars[mn]),
                  paste0("NDVI = ", vars[mx])),
       col = c(col1, col2, col3), 
       lwd = 2, lty = 1, bty = "n", 
       y.intersp = 1.5, horiz = F)


# find position of the terms associated with Tmin crossbasis
indt <- grep("basis_win", model$names.fixed)
# extract predictions from the Tmin DLNM centred on overall mean Tmin (19 deg C)
predt <- crosspred(basis_win, coef = coef[indt], vcov=vcov[indt,indt],
                   model.link = "log", bylag = 0.25, cen = round(median(data$win), 0)) 
# get exposures values
vars <- predt$predvar

# obtain relative risk (RR) fit and upper and lower confidence limits for all exposure variables
rr <- predt$matRRfit
rr.lci <- predt$matRRlow
rr.uci <- predt$matRRhigh

# set relative risk range 
r1 <- min(range(rr, rr.lci, rr.uci))
r2 <- max(range(rr, rr.lci, rr.uci))

summary(data$win)
round(quantile(data$win,0.025),1)
# get selected exposure variable positions
mn <- which(round(vars, 2) == round(quantile(data$win,0.05),1))
mx <- which(round(vars, 2) == round(quantile(data$win,0.95),1))


# define colours
col1 <- brewer.pal(11, "RdBu")[9]
tcol1 <- do.call(rgb, c(as.list(col2rgb(col1)), alpha = 255/4, max = 255))

col2 <- brewer.pal(11, "RdBu")[3]
tcol2 <- do.call(rgb, c(as.list(col2rgb(col2)), alpha = 255/4, max = 255))


# define x values (lag, by lag)
lagbylag <- seq(0, nlag, 0.25)

# cool
plot(lagbylag, rr[mn,], col = col1, type = "l", lwd = 1, 
     xlab = "Lag", ylab = "Relative risk", main = "", 
     ylim = range(r1, r2*1.11), frame.plot = T, axes = F)
axis(1, at = 0:nlag, labels = 0:nlag)
axis(2)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mn,], rev(rr.uci[mn,]))
polygon(xx, yy, col = tcol1, border = tcol1)
# warm
lines(lagbylag, rr[mx,], col = col2, lwd = 1)
xx <- c(lagbylag, rev(lagbylag))
yy <- c(rr.lci[mx,], rev(rr.uci[mx,]))
polygon(xx, yy, col = tcol2, border = tcol2)
abline(h = 1, lty = 3)



legend("topleft",
       legend = c(paste0("Wind speed = ",vars[mn],"m/s"),
                  paste0("Wind speed = ", vars[mx],"m/s")),
       col = c(col1, col2, col3), 
       lwd = 2, lty = 1, bty = "n", 
       y.intersp = 1.5, horiz = F)




########figure5-6#############
# !/usr/bin/env Rscript


for(yyear in 10:10){
  for(mmonth in 1:12){
    
    #inla.setOption(num.threads = "4:1")
    
    
    # Script to run INLA models in cross validation prediction mode
    
    # Step 0: load packages pre-processed data and functions
    # Step 1: rerun the selected model (fitted with config = TRUE for sampling) 
    # Step 2: produce cross-validated posterior predictive samples leaving out one year and one month at a time
    
    # Step 0: load packages and pre-processed data
    #source("01.output/6.preds_test/01-historic-load-all.R")
    
    # Step 2: produce cross-validated posterior predictive samples leaving out one year and one month at a time
    
    # define number of samples
    s <- 1000
    
    # replace dengue data in testing period with NA for out of sample prediction
    casestopred <- data$cases # response variable
    idx.pred <- which(data$year_index == yyear & data$month == mmonth)
    casestopred[idx.pred] <- NA # replace cases in year and month of interest to NA
    mpred <- length(idx.pred)
    
    # set response variable and year indicator
    df$Y <- casestopred
    
    mymodel <- function(formula, data = df, family = "nbinomial", config = TRUE)
      
    {
      model <- inla(formula = formula, data = data, family = family, offset = log(E),
                    control.inla = list(strategy = 'adaptive',cmin=0), 
                    control.compute = list(dic = TRUE, config = config, 
                                           cpo = TRUE, return.marginals = FALSE,
                                           waic=TRUE),
                    control.fixed = list(correlation.matrix = TRUE,
                                         mean.intercept=0, 
                                         prec.intercept=1, # tight precision
                                         mean=0, 
                                         prec=1),
                    control.predictor = list(link = 1, compute = TRUE), 
                    verbose = FALSE)
      model <- inla.rerun(model)
      return(model)
    }
    
    finalformula <- Y ~ 1 + f(T1, replicate = S, model = "rw1", cyclic = TRUE, constr = TRUE,
                              scale.model = TRUE,  hyper = hyper2.rw) +
      f(S, model = "bym2", graph = "china.graph",replicate = T2,constr = TRUE,
        scale.model = TRUE, hyper = hyper.bym2)+
      f(T2,replicate = S, model = "rw1", constr = TRUE,
        scale.model = TRUE,  hyper = hyper2.rw)+ 
      nvdi+lncar+beds+holiday+birth+lnpopden+
      basis_mt+basis_rhu+basis_win  
    
    model <- mymodel(finalformula, df)
    
    xx <- inla.posterior.sample(s, model)
    xx.s <- inla.posterior.sample.eval(function(...) c(theta[1], Predictor[idx.pred]), xx)
    y.pred <- matrix(NA, mpred, s)
    for(s.idx in 1:s) {
      xx.sample <- xx.s[, s.idx]
      y.pred[, s.idx] <- rnbinom(mpred, mu = exp(xx.sample[-1]), size = xx.sample[1])
    }
    preds <- list(year = 2008 + yyear, month = mmonth, idx.pred = idx.pred, 
                  mean = apply(y.pred, 1, mean), 
                  median = apply(y.pred, 1, median),
                  lci = apply(y.pred, 1, quantile, probs = c(0.025)),
                  uci = apply(y.pred, 1, quantile, probs = c(0.975)))
    save(preds, file = paste0("output/preds/preds_all_",2008 + yyear, "_", mmonth, ".RData"))
  }
}



##############plot#########################

pred<-data.frame()
for(yyear in 1:10){
  for(mmonth in 1:12){
    load(paste0("output/preds/preds_all_",2008 + yyear, "_", mmonth, ".RData"))
    preds<-data.frame(preds)
    pred<-rbind(pred,preds)   
    
  }
}

pred<-pred[with(pred,order(idx.pred)),]


figdata<-data.frame(data$ID,data$province,data$cases,data$pop,pred)
names(figdata)[1]<-'ID'
names(figdata)[2]<-'province'
names(figdata)[3]<-'cases'
names(figdata)[4]<-'pop'


library(ggplot2)
library(dplyr)
library(scales)

head(figdata)
str(figdata)
figdata$day<-1
figdata$year1<-figdata$year-1
figdata$date <- as.Date(paste(figdata$year,figdata$month,figdata$day,sep = '-'))

province<-unique(figdata$province)

class(province[i])

##############分面图##############

figdata$incidence<-figdata$cases/figdata$pop*100000
figdata$m<-figdata$mean/figdata$pop*100000
figdata$l<-figdata$lci/figdata$pop*100000
figdata$u<-figdata$uci/figdata$pop*100000

pdf("Fig_S7.pdf", width=15, height=15)
ggplot(data = figdata) + geom_line(mapping = aes(x = date, y = m),lwd = 1.2) + 
  geom_ribbon(mapping = aes(x = date, y = m, ymin = l, ymax = u), alpha = 0.5, fill = 'gray') + 
  geom_point(mapping = aes(x = date, y = incidence), color = 'red') + 
  xlab('') + ylab('') + 
  facet_wrap(~figdata$province, ncol = 5,
             scales = 'free')+
  theme(
    axis.line = element_line(size = 0.5), 
    panel.grid.minor  = element_blank(), panel.grid.major = element_blank(), 
    axis.text.y = element_text(size = 13,family = 'sans'), 
    axis.title.y.right = element_text(angle=90), axis.ticks = element_line(size = 0.5), 
    panel.background = element_blank(), panel.border = element_rect(fill = NA,size = 0.5),legend.position = '',
    axis.title = element_text(family = 'sans',size = 15),plot.title = element_text(hjust = 0.5 , size = 16),
    strip.text = element_text(size=16),
    axis.text.x = element_text(size = 13,family = 'sans',angle = 30))+
  scale_x_date(breaks  =seq.Date(from=as.Date("2009-01-01"),
                                 to=as.Date("2018-01-01"),
                                 by="2 year"),date_labels = '%Y')+
  labs(y='HFMD incidence')

dev.off()

xx<-read_excel("figdata1_1106.xlsx")
xx$ilci<-(xx$lci/xx$pop)*100000
xx$iuci<-(xx$uci/xx$pop)*100000
ggplot(data = xx) +geom_ribbon(aes(ymin = ilci,ymax=iuci,x = date),fill="red",alpha = 0.1)+ 
  geom_line(mapping = aes(x = date, y = incidence), color = '#FF0000',alpha = 0.7)+
  geom_line(mapping = aes(x = date, y = m,), color = "#006400",alpha = 0.7,linetype=4,size=0.5) + 
  geom_line(mapping = aes(x = date, y = outbreak), color = '#191970',alpha = 0.7,linetype=6,size=0.5) + 
  xlab('') + ylab('') + 
  facet_wrap(~xx$province, ncol = 6,
             scales = 'free')+theme(panel.background = element_blank(), axis.line = element_line(size = 0.5), 
                                    panel.grid.minor  = element_blank(), panel.grid.major = element_blank(), 
                                    axis.text.y = element_text(size = 13,family = 'sans'), 
                                    axis.title.y.right = element_text(angle=90), axis.ticks = element_line(size = 0.5))
  theme(
    axis.line = element_line(size = 0.5), 
    panel.grid.minor  = element_blank(), panel.grid.major = element_blank(), 
    axis.text.y = element_text(size = 13,family = 'sans'), 
    axis.title.y.right = element_text(angle=90), axis.ticks = element_line(size = 0.5), 
    panel.background = element_blank(), panel.border = element_rect(fill = NA,size = 0.5),legend.position = '',
    axis.title = element_text(family = 'sans',size = 15),plot.title = element_text(hjust = 0.5 , size = 16),
    strip.text = element_text(size=22),
    axis.text.x = element_text(size = 13,family = 'sans',angle = 30))+
  scale_x_date(breaks  =seq.Date(from=as.Date("2009-01-01"),
                                 to=as.Date("2018-01-01"),
                                 by="2 year"),date_labels = '%Y')+
  labs(y='HFMD incidence')+
  theme(strip.background =element_blank())

data<-xx %>% group_by(month,province) %>% mutate(outbreak=quantile(incidence,0.75))
data$outbreak<-data$outbreak/data$pop*100000