DLNM在时间序列中的应用：多个城市

###1.读取数据####
library(dlnm)
library(mvmeta)
library(splines)
alldata<-read.csv("dataset.csv" )

###2.数据库转换####
#为方便对每个地区进行循环建模，每个地区作为一个dataframe，并存在一个list内
cities <- unique(alldata$area) #获得每个地区的唯一值

##对每个地区进行切割
data <- lapply(cities, function(x){ 
  alldata[alldata$area==x,]
})
names(data) <- cities


###3.交叉基参数定义####
##暴露参数
varfun = "bs"
vardegree = 2
varper <- c(25,50,75) 
cenper <- 75

##滞后参数
lag <- 21
arglag <- list(knots=logknots(21,3))

###4.模型定义####
##长期趋势的自由度
dft <- 8
##模型公式
formula <-deathall~cb.tm + ns(rh,3) + factor(dow)+ns(time,df=dft*7)  

###5.模型结果保存设置####
#需要将各地区建立模型的系数和协方差提取出来，为加速运行速度，我们预先建立系数和协方差的储存框架
#用NA填充矩阵，因为用bs作为暴露参数，设置了2个degree和3个节点，因此有5个参数（即b1-b5）
coef<- matrix(NA,length(data),length(varper)+vardegree),
dimnames=list(cities,paste("b",seq(5),sep=""))) 

#协方差矩阵保存在list中
vcov <- vector("list",length(data))
names(vcova) <- cities

###6.建立循环，拟合各个地区的模型####

for(i in seq(length(data))) {
  # 输出当前地区
  print(paste(i,cities[i],sep = "-"))

  # 提取单个地区数据
  sub <- data[[i]]
  sub <- sub[order(sub$date),] #根据时间进行排序
  #定义长期趋势变量time
  sub$time<-c(1:length(sub$tm))
  # 设置交叉基 
  argvar <- list(fun=varfun,degree=vardegree,
                 knots=quantile(sub$tm,varper/100,na.rm=T))
  cb.tm <- crossbasis(sub$tm,lag=lag,argvar=argvar,arglag=arglag) #命名需与公式定义相对应

  # 运行模型
  model <- glm(formula.all,family=quasipoisson,sub,na.action="na.exclude")

  # 降维为累积维度（即仅保留滞后累积结果）
  red <- crossreduce(cb.tm,model,cen=quantile(sub$tm,cenper/100,na.rm=T))
  coef[i,] <- coef(red)
  vcov[[i]] <- vcov(red)
}

###1.meta合并####

###1.meta合并####
mv <- mvmeta(coef,vcov,method = "reml") #利用随机效应模型（reml）进行meta合并

###2.绘制曲线图####
##获得所有地区的温度范围
bound.all <-range(alldata$tm,na.rm=T)

perct<-unique(alldata$tm)
##将位点25%、50%和75%tm提出，设置为节点（与first stage相对应）
knots1<- rowMeans(sapply(data,function(x) {
  quantile(x$tm,c(25,50,75)/100,na.rm=TRUE)
}))
##根据整体tm范围，建立整体的tm序列##
tperc <- seq(bound.all[1],bound.all[2],length=50) #50为序列长度，可自行设置

##设置交叉基参数onebasis
#因为上面crossreduce已经降维，交叉基需要用onebasis，其余参数一致
argvar1 <- list(fun=varfun,degree=vardegree,
                knots= quantile(alldata$tm,varper/100,na.rm=T),
                Boundary.knots=bound.all)
cb.tm1 <- crossbasis(tperc,maxlag=lag,argvar=argvar1,arglag=arglag)
bvar <- do.call("onebasis",c(list(x=tperc),attr(cb.tm1,"argvar")))

##利用meta合并的系数对温度序列进行预测
all.cp <- crosspred(bvar,coef=coef(mv),vcov=vcov(mv),model.link="log",
                    at=c(perct,knots1) )

##寻找最小死亡风险温度
MMT <- round(as.numeric(names(which.min(all.cp$allfit))),1)

##设置MMT，重新预测
all.cp<- crosspred(bvar,coef=coef(mv),vcov=vcov(mv),model.link="log",
                       at=c(perct) ,cen = MMT)


##对结果进行展示
plot(all.cp,"overall",col=2,cex.axis=1.50,xlab="温度(℃)",ylab="RR",
     main="",cex.main=1.5,cex.lab=1.5) 

blup <- blup(mv,vcov =TRUE)

plot(all.cp,"overall",col=2,cex.axis=1.50,xlab="温度(℃)",ylab="RR",
     main="",cex.main=1.5,cex.lab=1.5)

for(i in seq(length(cities))) {
  prec <- crosspred(bvar,coef=blup[[i]]$blup,vcov=blup[[i]]$vcov,cen =MMT,at=c(perct,tperc),
                    model.link = "log")
  lines(prec,col="lightblue4",lty=5)
}
lines(all.cp,"overall",col=2,lwd=4)