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

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

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

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

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

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

1.1 DLNM模型

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

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

1.2 多元meta分析

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

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

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

2 實例與R軟件分析

2.1 資料概述

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

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

2.2 數(shù)據(jù)的導入與數(shù)據(jù)的預處理

library(readxl)#加載readxl包

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

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)建一個數(shù)據(jù)框，命名為ranges，該數(shù)據(jù)框是各地級市平均氣溫的最大值和最小值

2.3 對單個地級市利用“dlnm”與“splines”包構建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)

#建立交叉基。對自變量與因變量的關系、滯后效應分別選擇合適的基函數(shù)，2個基函數(shù)的張力積即得交叉基。lag設置最大滯后天數(shù)90 d；參數(shù)argvar控制平均氣溫與肺結核報告發(fā)病率的關系，自由度設置為3，基函數(shù)采用二次B樣條函數(shù)，參考值采用平均氣溫的中位數(shù)；參數(shù)arglag控制滯后效應，自由度設置為3，采用自然立方樣條函數(shù)。summary(cb)可以查看交叉基參數(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)

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

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

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圖展示平均氣溫對每日肺結核報告發(fā)病數(shù)的暴露—滯后—反應關系。crosspred()函數(shù)中，參數(shù)by=1是指計算所有整數(shù)平均氣溫的效應值。x軸為平均氣溫，y軸為滯后天數(shù)，z軸為RR值。

圖1 不同滯后期不同氣溫對日發(fā)病數(shù)效應分布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)

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

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

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)

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

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

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)

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

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

2.4 對各地級市進行模型建立

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

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

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

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

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

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

names(Sall)<-regions

#將模型中估計得到的各地級市的協(xié)方差結果儲存于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)

})#對總體累積匯總

suppressWarnings({

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

})#對熱效應匯總

suppressWarnings({

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

})#對冷效應匯總

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

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

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

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

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

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

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

}

})

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

2.5 進行多元meta分析

library(mvmeta)#加載mvmeta 包

method <-"reml" #選取隨機效應模型，若用固定效應模型，采用“fixed”

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

mvhot <-mvmeta(yhot～1,Shot,method=method)#將各地級市熱效應結果進行合并

mvcold <-mvmeta(ycold～1,Scold,method=method)#將各地級市冷效應結果進行合并

m <-length(regions)#共有9個地級市

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)

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

圖5 9個地級市平均氣溫對日發(fā)病數(shù)的累積效應分布圖

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

(qall <-qtest(mvall))#對整體合并進行Q檢驗

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個地級市熱效應合并圖，如圖6所示。9個地級市熱效應對日發(fā)病數(shù)的滯后合并效應無統(tǒng)計學意義。

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

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

(qhot <-qtest(mvhot))#對熱效應合并進行Q檢驗

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個地級市冷效應合并圖，如圖7所示。滯后0～72 d冷效應對日發(fā)病數(shù)的效應無統(tǒng)計學意義，自72 d后，RR值大于1，具有統(tǒng)計學意義，冷效應是日發(fā)病數(shù)的危險因素。

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

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

(qcold <-qtest(mvcold))#對冷效應合并進行Q檢驗

3 討論

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