使用R語言對兩種結(jié)核病發(fā)病預(yù)測模型的對比分析

汪 洋，胡代玉，楊德香，劉 瑛，王潤華，易 靜

2000年全國結(jié)核病流行病學(xué)調(diào)查顯示，我國目前約有5.5億人感染了結(jié)核分枝桿菌，我國現(xiàn)有活動性肺結(jié)核病患者約450萬人，每年新增結(jié)核病患者約150萬人，每年死于結(jié)核病的患者約有25萬，仍為各類傳染病之首[1]。根據(jù)結(jié)核病變化規(guī)律建立結(jié)核病疫情發(fā)展預(yù)測模型，逼近準(zhǔn)確的結(jié)核發(fā)病率，可以量化對結(jié)核的控制和管理的決策，是有效的結(jié)核病防治工作的前提。本研究采用免費(fèi)開源的R語言中的神經(jīng)網(wǎng)絡(luò)分析模型和灰色模型進(jìn)行結(jié)核病發(fā)病情況的分析與預(yù)測，現(xiàn)報(bào)道如下。

1 資料與方法

1.1 一般資料 資料來源于重慶市結(jié)核病防治所1993—2009年的結(jié)核病發(fā)病登記數(shù)據(jù)。

1.2 研究方法

1.2.1 R語言程序 R是由AT&T貝爾實(shí)驗(yàn)室所創(chuàng)的s語言發(fā)展出的有著統(tǒng)計(jì)分析及強(qiáng)大作圖功能的軟件系統(tǒng)。R語言內(nèi)含的作圖函數(shù)能將產(chǎn)生的圖片展示在一個(gè)獨(dú)立的窗口中，并能將之保存為各種形式的文件(jpg、png、bmp、ps、pdf、emf、pictex、xfig，具體形式取決于操作系統(tǒng))。統(tǒng)計(jì)分析的結(jié)果也能直接顯示出來，一些中間結(jié)果(如P值、回歸系數(shù)、殘差等)既可保存到專門的文件中，也可以直接用作進(jìn)一步的分析[2]。與SPSS、MATLAB等相比，R作為一個(gè)開源的免費(fèi)軟件，其價(jià)值得到不斷的延伸，主要表現(xiàn)在：可以跨平臺運(yùn)行，對矩陣的操作強(qiáng)大而高速，擁有許多可用的附加包，靈活的編程環(huán)境，對于已經(jīng)廣泛使用的統(tǒng)計(jì)軟件是可兼容的[3]。

1.2.2 BP神經(jīng)網(wǎng)絡(luò)的建模過程 BP神經(jīng)網(wǎng)絡(luò)通常有一個(gè)或多個(gè)sigmoid隱層和線性輸出層，能夠?qū)哂杏邢迋€(gè)不連續(xù)點(diǎn)的函數(shù)進(jìn)行逼近[4]。其學(xué)習(xí)過程由兩部分組成：正向傳播與反向傳播。正向傳播讓輸入信息在相應(yīng)權(quán)值、閾值和激活函數(shù)的作用下傳遞到輸出層，當(dāng)輸出的結(jié)果和期望值的誤差大于給定精度時(shí)，則將誤差反向傳播。在誤差返回過程中，網(wǎng)絡(luò)修正各層的權(quán)值和閾值。如此反復(fù)迭代，最后使傳遞信號的誤差達(dá)到允許精度。模型的設(shè)計(jì)包括：網(wǎng)絡(luò)類型及層數(shù)的確定；輸入及輸出變量的選擇；隱含層神經(jīng)元數(shù)目的確定；激活函數(shù)的選擇。常用的三層BP神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)見圖1。

圖1 三層BP網(wǎng)絡(luò)結(jié)構(gòu)

2 結(jié)果

2.1 結(jié)核病發(fā)病登記數(shù)據(jù) 重慶市結(jié)核病防治所1993—2009年的結(jié)核病發(fā)病登記數(shù)據(jù)見表1。

表1 重慶市1993—2009年的結(jié)核病發(fā)病登記數(shù)據(jù)(含涂陰、涂陽病例)

2.2 神經(jīng)網(wǎng)絡(luò)模型

2.2.1 數(shù)據(jù)準(zhǔn)備 使用重慶市結(jié)核病防治所1993—2009年的結(jié)核病發(fā)病登記數(shù)據(jù)，將1993—2007年的登記數(shù)據(jù)作為建模數(shù)據(jù)，2008—2009年的登記數(shù)據(jù)用于驗(yàn)證預(yù)測的準(zhǔn)確性。按照此次神經(jīng)網(wǎng)絡(luò)訓(xùn)練的模式設(shè)計(jì)，輸入層有3個(gè)結(jié)點(diǎn)，隱含層結(jié)點(diǎn)數(shù)為3，隱含層的激活函數(shù)為tansig；輸出層結(jié)點(diǎn)數(shù)為1個(gè)，輸出層的激活函數(shù)為logsig，設(shè)置學(xué)習(xí)速率為0.05，收斂誤差界值為0.005。

依據(jù)下列公式對原始數(shù)據(jù)進(jìn)行初始化。

Pmax=max〔P〕

Tmax=max〔T〕

Pmn=P/Pmax

Tmn=T/Tmax

集合變量P為1993—2007年結(jié)核病登記發(fā)病數(shù)，Pmax為集合中的最大值，同樣，T為輸出量集合，Tmax為輸出量集合中的最大值，Pmn和Tmn為集合P、T經(jīng)過歸一化處理后的實(shí)驗(yàn)數(shù)據(jù)，代碼如下：P=c(8135，10329，11174，12264，15781，14387，15674，16402，18018，17187，23505，28836，29644，28349，29119)；#1993—2007年發(fā)病數(shù)，賦值給P

Pmn=c()；#歸一化數(shù)據(jù)存儲變量

Pmax=max(P)；#取最大值

len=length(P)；#計(jì)算長度

for(i in 1：len){

Pmn[i]= P[i]/Pmax；

〕 #循環(huán)處理每一個(gè)歸一化元素

Pdata=matrix(nrow=12，ncol=3) #創(chuàng)建12*3的矩陣用于存儲訓(xùn)練數(shù)據(jù)的輸入集，以1993、1994、1995年的歸一化數(shù)據(jù)為一組，1994、1995、1996年的歸一化數(shù)據(jù)為一組，以此類推共12組數(shù)據(jù)。

k=1；

for(i in 1：12){

for(j in 1：3) {

Pdata[i，j]=Pmn[k]；

k= k+1；

}

k=k-2；

}

T=c(12264，15781，14387，15674，16402，18018，17187，23505，28836，29644，28349，29119)；#輸出集，為1996—2007年結(jié)核病登記病例數(shù)歸一化值。

Tmax=max(T)；

Tdata=c()；

len=length (T)；

for(i in 1：len) {

Tdata[i]=T[i] / Tmax；

}

2.2.2 建模與預(yù)測結(jié)果 有超過1 800個(gè)免費(fèi)的預(yù)測分析功能包供下載，地址：http：//cran.r-project.org[6]。此處需要：AMORE、Rserve，下載地址：http：//cran.csdb.cn/web/packages/available_packages_by_name.html。

利用R語言中的神經(jīng)網(wǎng)絡(luò)模型(AMORE)模塊建立神經(jīng)網(wǎng)絡(luò)模型的代碼如下：

library(AMORE)； #加載神經(jīng)網(wǎng)絡(luò)算法包

net=newff(n.neurons=c(3，3，1，1)，learning.rate.global=1e-2，momentum.global=0.5，

error.criterium=“LMS”，Stao=NA，hidden.layer=“tansig”，

output.layer=“purelin”，method=“ADAPTgd”)；#新建一個(gè)神經(jīng)網(wǎng)絡(luò)

result=train(net，Pdata，Tdata，error.criterium="LMS"，report=TRUE，show.step=100，n.shows=5)；#訓(xùn)練結(jié)果。

P2=c(29644，28349，29119，27098，25010)；#預(yù)測2008、2009年的結(jié)核病發(fā)病例數(shù)。為滿足輸入層有3個(gè)結(jié)點(diǎn)，隱含層結(jié)點(diǎn)數(shù)為3的條件，使用2005—2009年的結(jié)核病登記病例數(shù)作為輸入集。

P2max=max(P2)；

len=length(P2)；

P2mn=c()；

for(i in 1：len) {

P2mn[i]=P2[i]/P2max；

}

Pdata2=matrix(ncol=3，nrow=3)；

k=1；

for(i in 1：3) {

for(j in 1：3) {

Pdata2[i，j]=P2mn[k]；

k= k+1；

}

k= k-2；

}

y=sim(result$net，Pdata2)；

預(yù)測結(jié)果如下：

[，1]

[1，] 0.8707301

[2，] 0.8506458

[3，] 0.8381221

將2008—2009年的歸一化預(yù)測結(jié)果即[1，1]和[2，1]的值與Pmax相乘得到相應(yīng)的預(yù)測結(jié)，見表2。

表2 神經(jīng)網(wǎng)絡(luò)模型對重慶市2008年和2009年的結(jié)核病預(yù)測結(jié)果

Table2 Forecast results of tuberculosis by neural network model in Chongqing city in 2008 and 2009

年份登記病例數(shù)Pmax歸一化預(yù)測結(jié)果預(yù)測結(jié)果誤差200827098296440.87073025811.920124.75%200925010296440.85064625216.550020.83%

注：2010年的實(shí)際登記病例數(shù)未獲得，不能評價(jià)本模型的誤差

2.3 灰色模型(GM1，1)

2.3.1 模型函數(shù)準(zhǔn)備

R語言中暫無灰色模型，可手工實(shí)現(xiàn)R語言的灰色模型函數(shù)，其代碼如下：

gm<-function(x0，t){

x1<-cumsum(x0)

b<-numeric(length(x0)-1)

n<-length(x0)-1

for(i in 1：n){

b[i]<--(x1[i]+x1[i+1])/2

b}

D<-numeric(length(b))

D[]<-1

B<-cbind(b，D) BT<-t(B)

M<-solve(BT%*%B)

YN<-numeric(length(x1)-1)

YN<-x0[2：length(x0)]

alpha<-M%*%BT%*%YN

alpha2<-matrix(alpha，ncol=1)

a<-alpha2[1]

u<-alpha2[2]

cat("Model parameters："，′/n′，"a="，a，"u="，u，′/n′，′/n′)

y<-numeric(length(c(1：t)))

y[1]<-x1[1]

for(w in 1：(t-1)){

y[w+1]<-(x1[1]-u/a)*exp(-a*w)+u/a

}

cat("1AGO predictions："，′/n′，y，′/n′，′/n′)

xy<-numeric(length(y))

xy[1]<-y[1]

for(o in 2：t){

xy[o]<-y[o]-y[o-1]

}

cat("The predict values："，′/n′，xy，′/n′，′/n′)

y0<-numeric(length(x0))

y0[1]<-x1[1]

for (q in 1：n){

cc<-exp(-a*q)

y0[q+1]<-(x1[1]-u/a)*cc+u/a

}

xp<-numeric(length(x0))

m<-length(x0)

xp[1]<-y0[1]

for(j in 2：m){

xp[J]<-y0[J]-y0[j-1]

}

e<-numeric(length(x0))

for (l in 1：m){

e[l]<-xp[l]-x0[l]

}

cat("Residuals："，′/n′，e，′/n′，′/n′)

se<-sd(e)

sx<-sd(x0)

cv<-se/sx

cat("Test for raw data："，′/n′，cv，′/n′，′/n′)

if(cv<0.35) cat("prediction is good"，′/n′)

else cat("prediction is not good"，′/n′)

plot(xy，col=′blue′，type=′b′，pch=16，xlab=′Time series′，ylab=′Values′)

points(x0，col=′red′，type=′b′，pch=18)

legend(′topright′，c(′Predictions′，′Raw data′)，pch=c(16，18)，lty=l，col=c(′blue′，′red′))

}

2.3.2 預(yù)測結(jié)果 調(diào)用上節(jié)創(chuàng)建的灰色模型函數(shù)進(jìn)行結(jié)核病發(fā)病情況預(yù)測，代碼如下：gmData=c(8135，10329，11174，12264，15781，14387，15674，16402，18018，17187，23505，28836，29644，28349，29119)； #定義變量名：gmData，用于存儲1993—2007年結(jié)核病發(fā)病數(shù)，共15個(gè)數(shù)據(jù)。gm(gmData，17)；#調(diào)用函數(shù)，使用1993—2007年的作為訓(xùn)練樣本，預(yù)測未來兩年的結(jié)核病發(fā)病數(shù)，共17個(gè)數(shù)據(jù)。

預(yù)測出2008—2009年發(fā)病例數(shù)分別為34 437.39、37 475.49，結(jié)果見表3。

表3 灰色模型對重慶市2008年和2009年的結(jié)核病預(yù)測結(jié)果

Table3 Forecast results of tuberculosis by grey model in Chongqing city in 2008 and 2009

年份登記病例數(shù)預(yù)測結(jié)果誤差20082709834437.3927.08%20092501037475.4949.84%

3 討論

通過R語言的神經(jīng)網(wǎng)絡(luò)模型和灰色模型對結(jié)核病發(fā)病情況預(yù)測，結(jié)果發(fā)現(xiàn)神經(jīng)網(wǎng)絡(luò)模型對結(jié)核病發(fā)病人數(shù)的預(yù)測誤差不超過5%，而灰色模型對結(jié)核病發(fā)病數(shù)的預(yù)測誤差較大，分別為27.08%和49.84%。說明神經(jīng)網(wǎng)絡(luò)模型對于結(jié)核病發(fā)病情況預(yù)測的準(zhǔn)確性遠(yuǎn)遠(yuǎn)高于灰色模型，是預(yù)測結(jié)核病發(fā)病情況的較優(yōu)模型，這與易靜等[7]研究結(jié)果一致。

從預(yù)測結(jié)果來看，R語言的預(yù)測結(jié)果并不亞于其他價(jià)格高昂的商業(yè)軟件，商業(yè)軟件具備強(qiáng)大的數(shù)據(jù)統(tǒng)計(jì)、分析、預(yù)測和周邊數(shù)據(jù)處理等功能，若要完全或深入掌握商業(yè)軟件各項(xiàng)功能的使用，需要耗費(fèi)大量的時(shí)間和費(fèi)用進(jìn)行專門的培訓(xùn)。在實(shí)際分析預(yù)測的工作中，不會用到商業(yè)軟件的所有功能。

本文采用免費(fèi)開源的R語言中的神經(jīng)網(wǎng)絡(luò)分析模型和灰色模型完成分析與預(yù)測，并得到了滿意的預(yù)測結(jié)果，相對于使用商業(yè)軟件而言，可大大降低辦公或科研成本。

研究表明，R語言能滿足結(jié)核病發(fā)病數(shù)的預(yù)測需求，可推廣到其他傳染病發(fā)病數(shù)的預(yù)測或其他行業(yè)的數(shù)據(jù)分析與預(yù)測。下一步研究方向是根據(jù)決策需求和數(shù)據(jù)來源(如：《中國疾病預(yù)防控制信息系統(tǒng)(網(wǎng)絡(luò)直報(bào)信息系統(tǒng))》)將R語言的各個(gè)功能抽象化，通過高級編程語言(如：.net、Java等)將分析模型、數(shù)據(jù)源和需求有機(jī)地整合起來，提供一個(gè)圖形化操作界面，形成一套易用的決策支撐與預(yù)測系統(tǒng)，為衛(wèi)生資源調(diào)配、科研工作等提供可靠的決策支撐。隨著計(jì)算機(jī)硬件技術(shù)和軟件工具的持續(xù)發(fā)展，R語言也會在保持目前成就的情況下，通過不斷改善計(jì)算環(huán)境，成為有力的科研工具[8]。

1 張忠海.關(guān)于結(jié)核病防治現(xiàn)狀的思考[J].齊齊哈爾醫(yī)學(xué)院學(xué)報(bào)，2011，32(7)：1113-1114.

2 石發(fā)恩，任如山，蔣達(dá)華.R語言在消聲器設(shè)計(jì)試驗(yàn)中的應(yīng)用[J].金屬礦山，2008，20(11)：97.

3 Hiroyoshi Arai.A function for the R programming language to recast garnet analyses into end-members：Revision and porting of Muhling and Griffin S method[J].Computers& Geosciences，2010，36：406-409.

4 徐國祥.統(tǒng)計(jì)預(yù)測和決策[M].上海：上海財(cái)經(jīng)大學(xué)出版社，2001：113-131.

5 梁智勇.用Matlab實(shí)現(xiàn)GM(1，1)灰色模型的供電量預(yù)測[J].電腦編程技術(shù)與維護(hù)，2009，16(24)：93.

6 Robert I.Kabacoff.R in Action[M].New York：Manning Publications Co，2009：8.

7 易靜，杜昌廷，王潤華，等.彈性BP神經(jīng)網(wǎng)絡(luò)在結(jié)核病發(fā)病率預(yù)報(bào)中的應(yīng)用[J].現(xiàn)代預(yù)防醫(yī)學(xué)，2007，34(19)：3699-3701.

8 羅玫，趙嵩正，蔣建洪.模糊綜合評價(jià)模型的R語言實(shí)現(xiàn)[J].航空計(jì)算技術(shù)，2011，41(4)：56.