使用R語言對兩種結核病發病預測模型的對比分析

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

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

1 資料與方法

1.1 一般資料 資料來源于重慶市結核病防治所1993—2009年的結核病發病登記數據。

1.2 研究方法

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

1.2.2 BP神經網絡的建模過程 BP神經網絡通常有一個或多個sigmoid隱層和線性輸出層，能夠對具有有限個不連續點的函數進行逼近[4]。其學習過程由兩部分組成：正向傳播與反向傳播。正向傳播讓輸入信息在相應權值、閾值和激活函數的作用下傳遞到輸出層，當輸出的結果和期望值的誤差大于給定精度時，則將誤差反向傳播。在誤差返回過程中，網絡修正各層的權值和閾值。如此反復迭代，最后使傳遞信號的誤差達到允許精度。模型的設計包括：網絡類型及層數的確定；輸入及輸出變量的選擇；隱含層神經元數目的確定；激活函數的選擇。常用的三層BP神經網絡模型結構見圖1。

圖1 三層BP網絡結構

2 結果

2.1 結核病發病登記數據 重慶市結核病防治所1993—2009年的結核病發病登記數據見表1。

表1 重慶市1993—2009年的結核病發病登記數據(含涂陰、涂陽病例)

2.2 神經網絡模型

2.2.1 數據準備 使用重慶市結核病防治所1993—2009年的結核病發病登記數據，將1993—2007年的登記數據作為建模數據，2008—2009年的登記數據用于驗證預測的準確性。按照此次神經網絡訓練的模式設計，輸入層有3個結點，隱含層結點數為3，隱含層的激活函數為tansig；輸出層結點數為1個，輸出層的激活函數為logsig，設置學習速率為0.05，收斂誤差界值為0.005。

依據下列公式對原始數據進行初始化。

Pmax=max〔P〕

Tmax=max〔T〕

Pmn=P/Pmax

Tmn=T/Tmax

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

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

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

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

for(i in 1：len){

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

〕 #循環處理每一個歸一化元素

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

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年結核病登記病例數歸一化值。

Tmax=max(T)；

Tdata=c()；

len=length (T)；

for(i in 1：len) {

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

}

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

利用R語言中的神經網絡模型(AMORE)模塊建立神經網絡模型的代碼如下：

library(AMORE)； #加載神經網絡算法包

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”)；#新建一個神經網絡

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

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

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

預測結果如下：

[，1]

[1，] 0.8707301

[2，] 0.8506458

[3，] 0.8381221

將2008—2009年的歸一化預測結果即[1，1]和[2，1]的值與Pmax相乘得到相應的預測結，見表2。

表2 神經網絡模型對重慶市2008年和2009年的結核病預測結果

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

年份登記病例數Pmax歸一化預測結果預測結果誤差200827098296440.87073025811.920124.75%200925010296440.85064625216.550020.83%

注：2010年的實際登記病例數未獲得，不能評價本模型的誤差

2.3 灰色模型(GM1，1)

2.3.1 模型函數準備

R語言中暫無灰色模型，可手工實現R語言的灰色模型函數，其代碼如下：

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 預測結果 調用上節創建的灰色模型函數進行結核病發病情況預測，代碼如下：gmData=c(8135，10329，11174，12264，15781，14387，15674，16402，18018，17187，23505，28836，29644，28349，29119)； #定義變量名：gmData，用于存儲1993—2007年結核病發病數，共15個數據。gm(gmData，17)；#調用函數，使用1993—2007年的作為訓練樣本，預測未來兩年的結核病發病數，共17個數據。

預測出2008—2009年發病例數分別為34 437.39、37 475.49，結果見表3。

表3 灰色模型對重慶市2008年和2009年的結核病預測結果

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

年份登記病例數預測結果誤差20082709834437.3927.08%20092501037475.4949.84%

3 討論

通過R語言的神經網絡模型和灰色模型對結核病發病情況預測，結果發現神經網絡模型對結核病發病人數的預測誤差不超過5%，而灰色模型對結核病發病數的預測誤差較大，分別為27.08%和49.84%。說明神經網絡模型對于結核病發病情況預測的準確性遠遠高于灰色模型，是預測結核病發病情況的較優模型，這與易靜等[7]研究結果一致。

從預測結果來看，R語言的預測結果并不亞于其他價格高昂的商業軟件，商業軟件具備強大的數據統計、分析、預測和周邊數據處理等功能，若要完全或深入掌握商業軟件各項功能的使用，需要耗費大量的時間和費用進行專門的培訓。在實際分析預測的工作中，不會用到商業軟件的所有功能。

本文采用免費開源的R語言中的神經網絡分析模型和灰色模型完成分析與預測，并得到了滿意的預測結果，相對于使用商業軟件而言，可大大降低辦公或科研成本。

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

1 張忠海.關于結核病防治現狀的思考[J].齊齊哈爾醫學院學報，2011，32(7)：1113-1114.

2 石發恩，任如山，蔣達華.R語言在消聲器設計試驗中的應用[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 徐國祥.統計預測和決策[M].上海：上海財經大學出版社，2001：113-131.

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

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

7 易靜，杜昌廷，王潤華，等.彈性BP神經網絡在結核病發病率預報中的應用[J].現代預防醫學，2007，34(19)：3699-3701.

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