library(metafor)
library(timeDate)
library(lubridate)
library(plyr)
library(ggplot2)
library(lme4)
library(dplyr)
library(Hmisc)
library(tidyr)
library(stringr)
library(splines)
library(nlme)


#读取文件
mydata <- read.csv("F:/280区县死亡数据/new_280区县死因_匹配/2280死因数据_匹配_站点监测值_提取_2013-2017_19种pm25相关死亡因素_9个污染物.csv",header = T)

#查看字段名
names(mydata)


#创建新的数据框用来存放估计值、标准误差、ER值和95%CI
out<- data.frame(mcode=unique(mydata$newmcode),
                 B=NA,SE=NA,pval=NA,pi=NA,lpi=NA,upi=NA,model="glm",people="hbp")  #########修改1
#查看列名，设置列名，显示前几列
names(out)
rownames(out)<- NULL  #行名设置为空，即没有
head(out)



#####ER值和95%CI相关#######
intoPI <- function(B, unit){  
  #计算ER值和95%CI
  PI <- (exp(B * unit)-1)*100 
  return(PI)
}

intoCIs<- function(est,se){
  #计算最小值和最大值
  CIs<- est+c(-1,1)*1.96*se  #返回最小和最大值
  return(CIs)
}

intoeff<- function(B,unit,B_se){
  B_ci <- intoCIs(est = B, se = B_se) 
  out <- sapply(c(B, B_ci), intoPI, unit = unit)  #intopi 函数作用于ER值和95%CI，然后返回一个向量结果
  #names(out) <- c("point_pi", "lower_pi", "upper_pi")
  return(out)
}

names(mydata)
names(out)


####glm 模型运行####
fitglmall<-function(mydata,name,out,df_tem,df_rh,df_time){
  
  subdata<-subset(mydata,mydata$newmcode==name) %>% mutate(dow=NA,time=NA) #mutate创建新列dow ,time  name = out$mcode[i],i = 1时候，对应2013-2017年的110101
  subdata$date<- as.Date(subdata$date,"%Y/%m/%d")
  subdata$dow<- as.POSIXlt(subdata$date)$wday  #0=sunday
  #subdata$dow<-format(subdata$date,'%w') #转化对应的星期
  #subdata$time <- scale(subdata$date, scale = FALSE, center = TRUE)  #scale进行数据的中心化和标准化  Scale(x,center,scale)  [x—即需要标准化的数据, center—表示是否进行中心化,scale—表示是否进行标准化]  中心化就是数据集中的各项数据减去数据集的均值,标准化则是中心化之后的数据在除以数据集的标准差
  subdata$time<-seq(length=1826, from=1, by=1) 
  
  mod <- glm(hbp~pm25.24h + as.factor(dow) + ns(tem, df_tem) + ns(rh, df_rh) + ns(time, df_time*5),  ############修改2
             family=quasipoisson(link = "log"),data=subdata) 
  
  a <- subset(out,out$mcode==name)
  
  out[rownames(a),2:4]<-c(summary(mod)$coef[2,1],summary(mod)$coef[2,2],summary(mod)$coef[2,4])  #β，se，p值
  out[rownames(a),5:7]<-intoeff(out[rownames(a),2],10,out[rownames(a),3])
  
  return(out[rownames(a),])
}


####计算280各区县的估计值，标准误差，RR值，95%CI值

####计算280各区县的估计值，标准误差，RR值，95%CI值
#levels(as.factor(out$mcode))
for( i in 1:280){
  #循环过程中跳过出错的区县 23个
  if(i== 23 || i== 29 ||i==35 ||i==40 ||i== 42 ||i== 50 ||i==55 ||i==65 ||i==158 ||i==159 ||i==160 ||i==163 ||i== 194 ||i==202 ||i==207 ||i== 229 ||i==263 ||i==264 ||i==265 ||i==276 ||i==277 ||i==279 ||i==280){
    next
  }
  out[i,] <- fitglmall(mydata,out$mcode[i],out,3,3,7)  #时间自由度7
} 




write.csv(out,'C:/Users/zang/Desktop/quasipoisson/hbp_站点监测值.csv',row.names = F)   ##########修改3
#write.csv(out,'F:/280区县死亡数据/test_meta/统一区县_257个/时间自由度7-准泊松回归/站点监测值/lc_站点监测值.csv',row.names = F)   ##########修改3
#write.csv(out,'F:/280区县死亡数据/test_meta/res_TAP.csv',row.names = F)


#####meta分析####
#test1 <- read.csv("F:/280区县死亡数据/test_meta/统一区县_257个/时间自由度7-准泊松回归/站点监测值/lc_站点监测值.csv",header = T)   ###########修改4
test1 <- read.csv("C:/Users/zang/Desktop/quasipoisson/hbp_站点监测值.csv",header = T)   ###########修改4

vars = 'hbp'                                                                            ####################修改 5

metaresults2 <- data.frame(cat=vars, fit = NA, lower = NA, upper = NA)

regions <- as.character(unique(test1$people))  
#print(regions)
dlist <- lapply(regions, function(x) test1[test1$people==x,]) #生成类似nc结构文件
names(dlist) <- regions         #给上一步没有名字的文件改成regions的内容
print(length(dlist))
print(seq(length(dlist)))



for(i in seq(length(dlist))) {
  
  # PRINT
  cat(i,"")
  
  # EXTRACT THE DATA
  datarun <- dlist[[i]]
  
  
  fit_est<-datarun$B*10
  fit_se<-datarun$SE*10
  
  #固定效应模型
  metamod<-rma(yi=fit_est,sei=fit_se,data=datarun,method="REML") #yi表示效应值大小或结果测量，sei表示标准误的大小，data	表示包含该函数中yi，vi或sei等的数据框
  
  metaresults2[i,2]=predict(metamod, transf = exp, digits = 4)$pred
  metaresults2[i,3]=predict(metamod, transf = exp, digits = 4)$ci.lb
  metaresults2[i,4]=predict(metamod, transf = exp, digits = 4)$ci.ub
  
  #metaresults2$Lag<-"Lag0"
}  
write.csv(metaresults2,"C:/Users/zang/Desktop/quasipoisson/metaresult_hbp_站点监测值.csv",row.names = F)   ##########修改6