##############################################################Descriptive analysis
rm(list=ls())
library("dlnm")
library(dplyr)
library(lubridate)
library(splines)
library(ggplot2)
library("scales")
library(ggsci)
windowsFonts(A = windowsFont("Times New Roman"))
setwd('xxx\\path')

data <- read.csv('xxx\\Raw_data.csv')
summary(data)
str(data)

cor(data[,c(3,4,5,12,13,10,7,11)]) 


###########################################################DLNM analysis
fqaic <- function(model) {
  loglik <- sum(dpois(model$y,model$fitted.values,log=TRUE))
  phi <- summary(model)$dispersion
  qaic <- -2*loglik + 2*summary(model)$df[3]*phi
  return(qaic)
}

fqbic <- function(model) {
  loglik <- sum(dpois(model$y,model$fitted.values,log=TRUE))
  phi <- summary(model)$dispersion
  qbic <- -2*loglik + log(sum(model$data$visits))*summary(model)$df[3]*phi
  return(qbic)
}

#total
data$date <- as.Date(data$date)
basis.temp <- crossbasis(data$temperature, lag=7,
                         argvar=list(fun="ns",df = 4), arglag=list(fun='ns'),df = 4) 
model <- glm(visits ~ basis.temp +ns(sunshine_duration,3)+ ns(relative_humidity,3) + 
               ns(precipitation,3) +
               ns(air_presure,3) + ns(date,6*7) + as.factor(dow)+PM2.5,  #7=df，14=年
             family=quasipoisson(), data = data)
fqaic(model)
fqbic(model)


pred.temp <- crosspred(basis.temp, model, by = 0.01)
which.min(pred.temp$allRRfit)
pred.temp <- crosspred(basis.temp, model, by = 0.01,cen=22.76)


tiff("Figure2.tiff", units="in", width=6.2, height=5, res=300)
par(mar=c(5,5,4,4.5), 0.1)
par(family='A')
plot(pred.temp,"overall",col='red',ylab="RR",xlab="Daily Mean Temperature (℃)",xlim=c(-5,35), ylim=c(-0.5,6),axes=F,lwd = 1.5,font=2,font.lab=2)
axis(1,at=-1:7*5)
axis(2,at=c(1:6))
par(new=T)
hist(data$temperature,xlim=c(-5,35),ylim=c(0,1200),axes=F,ann=F,col=grey(0.95),breaks=30)
abline(v=quantile(data$temperature,c(0.01,0.99)),lty=2)
abline(v=22.76,lty=3)
axis(4,at=0:2*100)
mtext("Freq",4,line=2.5,at=200,cex=0.8,font=2,font.lab=2)
dev.off()


#################################################Attribution analysis
source("C:\\Users\\86130\\Desktop\\答复信\\code\\attrdl.R")

cen=22.76
attrdl(data$temperature,basis.temp,data$visits,model,type="an",cen=cen)   #计算总归因人数
attrdl(data$temperature,basis.temp,data$visits,model,cen=cen)*100             #计算总归因分值


##################################################model power estimation
####power.str from Armstrong BG, Gasparrini A, Tobias A, Sera F. Sample size issues in time series regressions of counts on environmental exposures. BMC Med Res Methodol. 2020 Jan 28;20(1):15.
power.tsr <- function(ncases,usablesdx=1,coef,overdispersion=1,alpha=0.05) {
  se<-sqrt(overdispersion)/(sqrt(ncases)*usablesdx)
  zetaalpha2 <- qnorm(1-alpha/2)
  cum1<-1-pnorm(zetaalpha2+(coef/se))
  cum2<-pnorm(coef/se-zetaalpha2)
  power<-cum1+cum2
  results <- cbind(ncases,se, 100*power)  
  colnames(results) <- c("ncases","SE","Power%")
  fixedparameters <- c(usablesdx,coef,overdispersion,alpha)
  names(fixedparameters) <- c("usablesdx","coef","overdispersion","alpha")
  output <- list(fixedparameters,results)
  names(output) <- c("fixedparamters","results")
  return(output)  # list of results (matrix) and fixed input parameters (vector)
}

power1 <- power.tsr(ncases=10747, usablesdx=0.7,coef=0.05) 
power1