陳東真 尹 嘉 丁國永 劉志東 劉雪娜 陸 華 李曉梅

1 山東第一醫(yī)科大學(xué)(山東省醫(yī)學(xué)科學(xué)院)，271016 山東 泰安；2 聊城市疾病預(yù)防控制中心，252000 山東 聊城；3山東大學(xué)齊魯醫(yī)院，250012 山東 濟(jì)南；4 泰安市疾病預(yù)防控制中心，271001 山東 泰安

環(huán)境暴露對(duì)健康的效應(yīng)一般采用多個(gè)研究地點(diǎn)的數(shù)據(jù)進(jìn)行時(shí)間序列分析，這種情況通常采用2階段分析方法。首先采用回歸模型分析暴露反應(yīng)關(guān)系，然后采用meta分析將多個(gè)研究地點(diǎn)的效應(yīng)進(jìn)行合并。而環(huán)境暴露對(duì)健康產(chǎn)生的效應(yīng)一般會(huì)存在滯后，且暴露與結(jié)局的關(guān)系呈現(xiàn)非線性關(guān)系，因此在第一階段常采用分布滯后非線性模型(distributed lag non-linear model, DLNM)來模擬環(huán)境暴露與健康結(jié)局的關(guān)系。DLNM是在廣義線性模型、廣義相加模型的基礎(chǔ)上，加入交叉基函數(shù)，可以同時(shí)靈活地評(píng)估時(shí)間序列數(shù)據(jù)中非線性暴露—反應(yīng)關(guān)系以及滯后效應(yīng)[1]。針對(duì)非線性暴露—反應(yīng)關(guān)系，在第二階段可采用多元meta分析[2]，對(duì)各研究地區(qū)估計(jì)的參數(shù)進(jìn)行合并，獲取環(huán)境暴露對(duì)健康的總效應(yīng)。

楊金建等[2]針對(duì)DLNM以及多元meta模型的原理進(jìn)行了詳細(xì)闡述，但是在實(shí)例分析部分并未展現(xiàn)結(jié)果實(shí)現(xiàn)的具體步驟。因此，本研究使用R軟件，以2005—2016年山東省9個(gè)地級(jí)市平均氣溫與肺結(jié)核報(bào)告發(fā)病數(shù)據(jù)為實(shí)例，利用“dlnm”包實(shí)現(xiàn)DLNM模型和“mvmeta”包實(shí)現(xiàn)多元meta分析。

1 DLNM模型及多元meta分析方法簡介

為獲取環(huán)境暴露對(duì)健康的總效應(yīng)，可采用混合式分析方法進(jìn)行，包括非線性暴露—反應(yīng)關(guān)系的DLNM模型與提取，并合并各地區(qū)參數(shù)值的多元meta分析。

1.1 DLNM模型

廣義相加模型只考慮暴露反應(yīng)關(guān)系，分布滯后線性模型沒有解決暴露與結(jié)局的非線性關(guān)系，由此發(fā)展了同時(shí)在暴露和滯后維度引入非線性關(guān)系的DLNM[3]。DLNM是以交叉基為核心，暴露—反應(yīng)關(guān)系與滯后效應(yīng)分別選取相應(yīng)的基函數(shù)，計(jì)算2組基函數(shù)的張力積得到交叉基函數(shù)。基函數(shù)主要包括多項(xiàng)式函數(shù)、自然立方樣條函數(shù)、線性閾值函數(shù)、懲罰樣條函數(shù)，最常用的是自然立方樣條和懲罰樣條。DLNM基本表達(dá)式：

其中，μt=E(Yt)，Yt代表結(jié)局變量，該結(jié)局變量符合指數(shù)分布族，例如正態(tài)分布、Gamma分布和Poisson分布等，E(Yt)代表t時(shí)刻因變量Y的期望；g代表連接函數(shù)；sj代表xj與E(Yt)間的非線性函數(shù)；uk代表其他與E(Yt)間存在線性關(guān)系的變量；β、γ分別代表xj、uk的參數(shù)向量。sj代表自變量xj的基函數(shù)，通過選擇不同的基函數(shù)可將xj轉(zhuǎn)換成一個(gè)新的變量集，包含到模型的設(shè)計(jì)矩陣中，并進(jìn)行參數(shù)估計(jì)。

1.2 多元meta分析

meta分析是對(duì)多個(gè)研究目的相同的研究結(jié)果進(jìn)行綜合分析，以使用量化的平均效果來解決研究問題的一種研究方法[4-5]。meta分析通常假設(shè)各研究之間是獨(dú)立的，而許多研究之間并不是相互獨(dú)立的，因此在meta分析的基礎(chǔ)上擴(kuò)展出了多元meta分析。

多元meta分析可以利用研究之間和/或研究內(nèi)部的相關(guān)結(jié)構(gòu)產(chǎn)生更有效或更精確的估計(jì)，此方法考慮到了各結(jié)果之間的相關(guān)性，克服了單變量meta分析方法忽略相關(guān)結(jié)構(gòu)的局限性[6]。對(duì)于隨機(jī)或者信息缺失，需要匯總數(shù)據(jù)結(jié)局的研究，該方法是更好的選擇。使用多元meta分析進(jìn)行分析時(shí)，主要包括隨機(jī)效應(yīng)模型與固定效應(yīng)模型。異質(zhì)性檢驗(yàn)使用Q檢驗(yàn)中的P值，來選擇效應(yīng)模型。在第一階段得到的研究內(nèi)協(xié)方差矩陣，包括每個(gè)研究效應(yīng)的方差及其協(xié)方差，在多元meta分析階段時(shí)，將這些估計(jì)值進(jìn)行合并[7]。

本研究靈活利用R語言“dlnm”包與“mvmeta”包實(shí)現(xiàn)DLNM模型與多元meta分析。通過實(shí)例分析，展現(xiàn)采用本研究的R語言代碼，利用原數(shù)據(jù)快速、直接地實(shí)現(xiàn)DLNM模型與多元meta分析相結(jié)合的兩階段統(tǒng)計(jì)分析方法。

2 實(shí)例與R軟件分析

2.1 資料概述

現(xiàn)以山東省威海市、煙臺(tái)市、青島市、東營市、濰坊市、日照市、濱州市、淄博市和臨沂市的氣候變量與肺結(jié)核報(bào)告發(fā)病數(shù)據(jù)集為例，演示DLNM與多元meta分析，合并多個(gè)地級(jí)市的非線性暴露—反應(yīng)效應(yīng)。本研究數(shù)據(jù)集包括每日平均氣溫、相對(duì)濕度、日照時(shí)數(shù)、平均風(fēng)速、降水量和肺結(jié)核報(bào)告發(fā)病數(shù)，見表1。

表1 2005—2016年平均氣溫對(duì)肺結(jié)核發(fā)病影響的數(shù)據(jù)資料

2.2 數(shù)據(jù)的導(dǎo)入與數(shù)據(jù)的預(yù)處理

library(readxl)#加載readxl包

Shandong<-read_excel("D:/Shandong.xlsx")#數(shù)據(jù)集導(dǎo)入

regions <-as.character(unique(Shandong$cityname))

data <-lapply(regions,function(x)Shandong[Shandong$cityname==x,])

names(data)<-regions

#創(chuàng)建地區(qū)序列的列表

ranges <-t(sapply(data, function(x)range(x$Atemp,na.rm=T)))

#創(chuàng)建一個(gè)數(shù)據(jù)框，命名為ranges，該數(shù)據(jù)框是各地級(jí)市平均氣溫的最大值和最小值

2.3 對(duì)單個(gè)地級(jí)市利用“dlnm”與“splines”包構(gòu)建DLNM模型

library(dlnm)#加載dlnm包

city <-"Qingdao" #選取青島市，命名為city

cen<-quantile(Atemp,0.50)

cb<-crossbasis(data[[city]]$Atemp,lag=90, argvar =list(df=3,fun="bs",cen=cen), arglag =list(df=3,fun="ns"))

summary(cb)

#建立交叉基。對(duì)自變量與因變量的關(guān)系、滯后效應(yīng)分別選擇合適的基函數(shù)，2個(gè)基函數(shù)的張力積即得交叉基。lag設(shè)置最大滯后天數(shù)90 d；參數(shù)argvar控制平均氣溫與肺結(jié)核報(bào)告發(fā)病率的關(guān)系，自由度設(shè)置為3，基函數(shù)采用二次B樣條函數(shù)，參考值采用平均氣溫的中位數(shù)；參數(shù)arglag控制滯后效應(yīng)，自由度設(shè)置為3，采用自然立方樣條函數(shù)。summary(cb)可以查看交叉基參數(shù)設(shè)置情況。

library(splines)#加載splines 包

model<-glm(PTB～cb+ns(time,7*12)+dow+Season+ns(SD,3)+ns(ARH,3)

+ns(AWS,3)+ns(PRE,3),family = quasipoisson(),data[[city]])

summary(model)

#模型的建立。調(diào)用自然立方樣條函數(shù)，運(yùn)用quasipoisson控制Poisson回歸的過離散現(xiàn)象，時(shí)間的自由度選取最常用的7，氣象變量的自由度選取3。

bound <-colMeans(ranges)#設(shè)置邊界值

cp <-crosspred(cb,model,from=bound[1],to=bound[2],by=1)

d3 <-plot(cp,xlab="Tempeature/℃", ylab="Lag/d", zlab="RR",theta=200, phi=40, lphi=150, shade=0.4)

#建立3-D圖，如圖1所示。3-D圖展示平均氣溫對(duì)每日肺結(jié)核報(bào)告發(fā)病數(shù)的暴露—滯后—反應(yīng)關(guān)系。crosspred()函數(shù)中，參數(shù)by=1是指計(jì)算所有整數(shù)平均氣溫的效應(yīng)值。x軸為平均氣溫，y軸為滯后天數(shù)，z軸為RR值。

圖1 不同滯后期不同氣溫對(duì)日發(fā)病數(shù)效應(yīng)分布3-D圖

P90<-as.numeric(quantile(Atemp,0.90))

crvar<-crossreduce(cb,model,type="var",value=P90,from=bound[1], to=bound[2],bylag=0.2)

plot(crvar,xlab="Lag",ylab="RR",col=2,lwd=2,tcl=0.5)

mtext(text=paste("Predictor-specific association at temperature ",P90, "℃",sep=""),cex=1)

#繪制特定平均氣溫時(shí)滯后—反應(yīng)關(guān)系圖，如圖2所示。橫坐標(biāo)代表滯后天數(shù)，縱坐標(biāo)代表RR值。滯后0 d時(shí)，平均氣溫第90百分位數(shù)(25.9 ℃)對(duì)日發(fā)病數(shù)的影響RR值小于1，隨著滯后天數(shù)的增加RR值逐漸變大，滯后90 d時(shí)RR值大于1，但均無統(tǒng)計(jì)學(xué)意義。

圖2 平均氣溫第90百分位數(shù)(25.9 ℃)時(shí)對(duì)日發(fā)病數(shù)的滯后效應(yīng)分布圖

crlag <-crossreduce(cb,model,type="lag",value=90,from=bound[1],

to=bound[2],bylag=0.2)

plot(crlag,xlab="Tempeature/℃",ylab="RR",col=2,ylim=c(0.9,1.15),lwd=2,tcl=0.5)

mtext(text="Lag-specific association at lag 90",cex=1)

#繪制特定滯后時(shí)間時(shí)暴露—反應(yīng)關(guān)系圖，如圖3所示。橫坐標(biāo)代表平均氣溫，縱坐標(biāo)代表RR值。滯后90 d時(shí)，平均氣溫對(duì)日發(fā)病數(shù)影響效應(yīng)無統(tǒng)計(jì)學(xué)意義。

圖3 滯后90 d平均氣溫與日發(fā)病數(shù)的暴露—反應(yīng)關(guān)系圖

crall <-crossreduce(cb,model,from=bound[1],to=bound[2],by=0.2)

plot(crall,xlab="Tempeature/℃",ylab="RR",ylim=c(0,40),col=2,lwd=2,tcl=0.5)

mtext(text="Overall cumulative association",cex=1)

#繪制累積效應(yīng)圖，如圖4所示。將每個(gè)滯后時(shí)間的滯后效應(yīng)相加便得到累積滯后效應(yīng)。不同平均氣溫對(duì)日發(fā)病數(shù)的累積效應(yīng)不同，但均無統(tǒng)計(jì)學(xué)意義。

圖4 不同平均氣溫對(duì)日發(fā)病數(shù)的累積90 d效應(yīng)分布圖

2.4 對(duì)各地級(jí)市進(jìn)行模型建立

lag <-c(0,90)#滯后范圍0～90 d

argvar <-list(df=3,fun="bs",cen=cen)#自變量與因變量的關(guān)系，自由度設(shè)置為3，基函數(shù)采用二次B樣條函數(shù)，參考值采用平均氣溫的中位數(shù)；

arglag <-list(df=3,fun="ns")#參數(shù)arglag控制滯后效應(yīng)，自由度設(shè)置為3，采用自然立方樣條函數(shù)

yall<-matrix(NA,length(data),3,dimnames=list(regions,paste("b",seq(3),sep="")))

#將模型中估計(jì)得到的各地級(jí)市系數(shù)結(jié)果儲(chǔ)存于yall矩陣中

Sall <-vector("list",length(data))

names(Sall)<-regions

#將模型中估計(jì)得到的各地級(jí)市的協(xié)方差結(jié)果儲(chǔ)存于Sall矩陣中

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)}

#求Q-AIC值

system.time({

for(i in seq(data)){

sub <-data[[i]]

suppressWarnings({

cb <-crossbasis(sub$Atemp,lag=lag,argvar=argvar,arglag=arglag)

})#建立交叉基

mfirst<-glm(PTB～cb+ns(time,7*12)+DOW+Season+ns(SD,3)+ns(ARH,3)+ns(AWS,3)+ns(PRE,3),family = quasipoisson(),sub)#建立模型

suppressWarnings({

crall <-crossreduce(cb,mfirst)

})#對(duì)總體累積匯總

suppressWarnings({

crhot <-crossreduce(cb,mfirst,type="var",value=25.9)

})#對(duì)熱效應(yīng)匯總

suppressWarnings({

crcold <-crossreduce(cb,mfirst,type="var",value=-1)

})#對(duì)冷效應(yīng)匯總

yall[i,]<-coef(crall)#整體模型系數(shù)

Sall[[i]]<-vcov(crall)#整體協(xié)方差

yhot[i,]<-coef(crhot)#熱效應(yīng)系數(shù)

Shot[[i]]<-vcov(crhot)#熱效應(yīng)協(xié)方差

ycold[i,]<-coef(crcold)#冷效應(yīng)系數(shù)

Scold[[i]]<-vcov(crcold)#冷效應(yīng)協(xié)方差

qaic[i]<-fqaic(mfirst)#求模型Q-AIC

}

})

#此為循環(huán)語句，對(duì)各地級(jí)市進(jìn)行模型建立、提取模型系數(shù)以及協(xié)方差，為多元meta分析效應(yīng)合并做鋪墊。

2.5 進(jìn)行多元meta分析

library(mvmeta)#加載mvmeta 包

method <-"reml" #選取隨機(jī)效應(yīng)模型，若用固定效應(yīng)模型，采用“fixed”

mvall <-mvmeta(yall～1,Sall,method=method)#將各地級(jí)市整體結(jié)果進(jìn)行合并

mvhot <-mvmeta(yhot～1,Shot,method=method)#將各地級(jí)市熱效應(yīng)結(jié)果進(jìn)行合并

mvcold <-mvmeta(ycold～1,Scold,method=method)#將各地級(jí)市冷效應(yīng)結(jié)果進(jìn)行合并

m <-length(regions)#共有9個(gè)地級(jí)市

xvar <-seq(bound[1],bound[2],by=0.1)

bvar <-do.call("onebasis",c(list(x=xvar),attr(cb,"argvar")))

cpall <-crosspred(bvar,coef=coef(mvall),vcov=vcov(mvall),

model.link="log",by=0.1,from=bound[1],to=bound[2])

tperc <-seq(bound[1],bound[2],by=0.1)

bperc <-do.call("onebasis",c(list(x=tperc),attr(cb,"argvar")))

plot(cpall,type="n",ylab="RR",ylim=c(0,30),xlab="Temperature/℃",tcl=0.5)

for(i in seq(m)){

suppressWarnings(

lines(crosspred(bperc,coef=yall[i,],vcov=Sall[[i]],model.link="log"),

col=grey(0.8),lty=5)

)

}

lines(cpall,"overall",col=1,lwd=2)

abline(h=1)

lines(cpall,col=2,lwd=2)

mtext("Main model: first-stage and pooled estimates",cex=1)

legend("top",c("Pooled(with 95%CI)","First-stage region-specific"),

lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)

#將所有地級(jí)市在各具體氣溫情況下的累計(jì)效應(yīng)進(jìn)行合并，繪制合并累積效應(yīng)圖，如圖5所示。9個(gè)地級(jí)市平均氣溫對(duì)日發(fā)病數(shù)的合并累計(jì)效應(yīng)先下降再上升，但無統(tǒng)計(jì)學(xué)意義。

圖5 9個(gè)地級(jí)市平均氣溫對(duì)日發(fā)病數(shù)的累積效應(yīng)分布圖

round(with(cpall,cbind(allRRfit,allRRlow,allRRhigh)["10",]),3)#不同氣溫的累積效應(yīng)

(qall <-qtest(mvall))#對(duì)整體合并進(jìn)行Q檢驗(yàn)

plot(cphot,type="n",ylab="RR",ylim=c(0.95,1.1),xlab="Lag",tcl=0.5)

xlag <-0:810/9

blag <-do.call("onebasis",c(list(x=xlag),attr(cb,"arglag")))

reghot <-apply(yhot,1,function(x)exp(blag%*%x))

matplot(xlag,reghot,type="l",col=grey(0.5),lty=2,add=T)

abline(h=1)

lines(cphot,col=2,lwd=2)

legend("top",c("Pooled(with 95%CI)","First-stage region-specific"),

lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)

mtext(text=paste("Predictor-specific summary for temperature = ",25.9,

"℃",sep=""),cex=1)

#繪制9個(gè)地級(jí)市熱效應(yīng)合并圖，如圖6所示。9個(gè)地級(jí)市熱效應(yīng)對(duì)日發(fā)病數(shù)的滯后合并效應(yīng)無統(tǒng)計(jì)學(xué)意義。

圖6 熱效應(yīng)(平均氣溫25.9 ℃)時(shí)9個(gè)城市日發(fā)病數(shù)的滯后效應(yīng)分布圖

round(with(cphot,cbind(matRRfit,matRRlow,matRRhigh)["90",]),3)#在熱效應(yīng)時(shí)，不同滯后時(shí)間的效應(yīng)

(qhot <-qtest(mvhot))#對(duì)熱效應(yīng)合并進(jìn)行Q檢驗(yàn)

plot(cpcold,type="n",ylab="RR",ylim=c(0.95,1.1),xlab="Lag",tcl=0.5)

xlag <-0:810/9

blag <-do.call("onebasis",c(list(x=xlag),attr(cb,"arglag")))

regcold <-apply(ycold,1,function(x)exp(blag%*%x))

matplot(xlag,regcold,type="l",col=grey(0.5),lty=2,add=T)

abline(h=1)

lines(cpcold,col=2,lwd=2)

legend("top",c("Pooled(with 95%CI)","First-stage region-specific"),

lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)

mtext(text=paste("Predictor-specific summary for temperature = ",-1, "℃",sep=""),cex=1)

#繪制9個(gè)地級(jí)市冷效應(yīng)合并圖，如圖7所示。滯后0～72 d冷效應(yīng)對(duì)日發(fā)病數(shù)的效應(yīng)無統(tǒng)計(jì)學(xué)意義，自72 d后，RR值大于1，具有統(tǒng)計(jì)學(xué)意義，冷效應(yīng)是日發(fā)病數(shù)的危險(xiǎn)因素。

圖7 冷效應(yīng)(平均氣溫-1℃)時(shí)9個(gè)城市日發(fā)病數(shù)的滯后效應(yīng)分布圖

round(with(cpcold,cbind(matRRfit,matRRlow,matRRhigh)["90",]),3)#冷效應(yīng)時(shí)，不同滯后時(shí)間的效應(yīng)。

(qcold <-qtest(mvcold))#對(duì)冷效應(yīng)合并進(jìn)行Q檢驗(yàn)

3 討論

本研究結(jié)合實(shí)例，介紹了DLNM與多元meta分析在R軟件中的實(shí)現(xiàn)過程。2階段時(shí)間序列分析是環(huán)境流行病學(xué)中一個(gè)強(qiáng)有力的分析方法，在復(fù)雜關(guān)聯(lián)的研究中提供了更大的靈活性且有更高的效率，可以減少選擇性結(jié)果報(bào)告偏倚，并允許對(duì)各結(jié)果之間的關(guān)聯(lián)進(jìn)行建模，實(shí)現(xiàn)對(duì)結(jié)果的聯(lián)合推斷以及從一個(gè)結(jié)果預(yù)測另一個(gè)結(jié)果[8]。DLNM結(jié)合多元meta分析的2階段分析，被國內(nèi)外多項(xiàng)研究運(yùn)用于氣象因素、污染因素與疾病健康之間關(guān)系的研究[9-13]，與單一地區(qū)研究相比，將多個(gè)研究地區(qū)結(jié)果進(jìn)行合并，可以減少研究的異質(zhì)性，而且多元meta分析提供了更好統(tǒng)計(jì)特性的參數(shù)估計(jì)，特別是通過調(diào)整研究間協(xié)方差結(jié)構(gòu)的估計(jì)，潛在地提高了研究結(jié)果的精度[7]。此外，利用R軟件繪制各研究地點(diǎn)氣象因素與疾病的暴露—反應(yīng)曲線，以及合并后的曲線圖，更加直觀地探討氣象因素對(duì)疾病發(fā)病影響的非線性關(guān)系以及滯后效應(yīng)。