田 剛

（天津市兒童醫(yī)院，天津 300134）

藥品生產(chǎn)管理者在制定生產(chǎn)計(jì)劃時要以銷定產(chǎn)，為此，要根據(jù)以往市場需求量來預(yù)測未來市場需求量，但單純依靠經(jīng)驗(yàn)預(yù)測未來需求量的做法存在一定局限性。這是醫(yī)院制劑室和藥廠普遍遇到的問題。為獲得精準(zhǔn)的預(yù)測，傳統(tǒng)上一般根據(jù)管理經(jīng)濟(jì)學(xué)提供的ARIMA模型和Holt-Winters 模型理論借助于SPSS 軟件計(jì)算出來。但該軟件提供的自動建模功能較少。而R 語言則能提供多種強(qiáng)大的函數(shù)來實(shí)現(xiàn)SPSS 軟件實(shí)現(xiàn)不了的功能，但缺點(diǎn)在于技術(shù)資料少，目前國內(nèi)很少見到將R 語言應(yīng)用到藥品管理領(lǐng)域的文獻(xiàn)。本次研究嘗試使用R 語言編程進(jìn)行歷史數(shù)據(jù)分析并建模和預(yù)測，從而摸索出一套數(shù)據(jù)分析模式，再將這種方法推廣到本院其他制劑生產(chǎn)預(yù)測中，改進(jìn)藥品管理決策。

1 研究工具

本次研究使用R 語言軟件（版本3.6.0）提供的auto.arima、naive、snaive、ets、hw 這5 種自動建模函數(shù)。其中，auto.arima 函數(shù)用于自動尋找最佳ARIMA 模型，naive 函數(shù)用于建立不帶漂移項(xiàng)的隨機(jī)游走模型并返回預(yù)測數(shù)據(jù)，snaive 函數(shù)用于建立季節(jié)性隨機(jī)游走模型并返回預(yù)測數(shù)據(jù)。此外本次研究還將使用2 種Holt-Winters 模型的建模函數(shù)，一個是ets 函數(shù)，另一個是hw函數(shù)。hw 函數(shù)是ets 函數(shù)的便捷包裝，函數(shù)中有事先默認(rèn)設(shè)定的參數(shù)值［1］。

2 研究方法

2.1 研究模型的建立 本次研究選取本院制劑室生產(chǎn)的小兒鹽酸麻黃堿滴鼻液作為研究兒童用藥銷售規(guī)律的切入點(diǎn)，采集在2012 年1 月—2019 年12 月每月出庫數(shù)量作為時間序列的數(shù)據(jù)，以2012 年1 月—2019 年6 月數(shù)據(jù)作為歷史數(shù)據(jù)建模并預(yù)測未來6 個月的需求量，以2019 年7—12 月數(shù)據(jù)用于檢驗(yàn)預(yù)測結(jié)果的對照。本次研究將使用上述5 種函數(shù)分別對歷史數(shù)據(jù)自動建立ARIMA 模型和Holt-Winters 模型來預(yù)測6 個月的市場趨勢，用W檢驗(yàn)對這些模型的平穩(wěn)性進(jìn)行檢查，使用MAPE 函數(shù)比較預(yù)測結(jié)果與對照數(shù)據(jù)的誤差，選取誤差最小的模型為最佳模型，從而建立起時間序列預(yù)測的方法學(xué)。

2.2 數(shù)據(jù)來源和預(yù)處理 本次研究將小兒鹽酸麻黃堿滴鼻液在2012 年1 月—2019 年12 月期間每月從制劑室出庫到藥房的銷量折算成每月出庫的生產(chǎn)批次數(shù)量，見表1。將這些出庫批數(shù)部分保存在y.csv 文件中，示意圖見圖1（限于篇幅，圖1 僅列舉前6 個月的數(shù)據(jù)）。

表1 源數(shù)據(jù)

代碼和解釋如下：

install.packages（"Metrics"）#下載并安裝Metrics軟件包；

install.packages（"forecast"） #下載并安裝forecast軟件包；

library（Metrics）#加載Metrics 軟件包；

library（forecast）#加載forecast 軟件包；

yaopin = read.csv ('s：\y.csv',header=T) #讀取S 盤的y.csv 文件并賦值給“yaopin”的數(shù)據(jù)框；

shixubiao ＜-ts（yaopin，frequency=12，start=2012）#將yaopin 數(shù)據(jù)框的數(shù)據(jù)建立時間序列并將結(jié)果賦值給“shixubiao”的向量；

shuju＜-window（shixubiao，start=2012，end=2019+5/12）#將2012 年1 月—2019 年6 月數(shù)據(jù)用于建模，并用window 函數(shù)將數(shù)據(jù)賦值給“shuju”的向量；

shijizhi＜-window (shixubiao,start=2019+6/12) # 將2019 年7—12 月數(shù)據(jù)賦予“shijizhi”的向量，用于與預(yù)測值的比對。

圖1 y.csv 文件的內(nèi)容

3 研究步驟

3.1 數(shù)據(jù)初步考察

3.1.1 繪制時序圖 求得數(shù)據(jù)框yaopin 的8 年所有數(shù)據(jù)平均值，繪制向量shixubiao 的時序圖，從而獲得對整個數(shù)據(jù)的初步直觀印象。即：每月該藥品的出庫批數(shù)起伏很大，但在2012—2017 年基本保持在一個范圍之中。從2018 年下半年開始到2019 年底，出庫批次有了明顯的增加。具體見圖2。

圖2 每月出庫批數(shù)時序圖

代碼如下：

y1＜- as.numeric（unlist（yaopin））#轉(zhuǎn)換為數(shù)值型向量并賦值給y1；

plot（shixubiao，xaxs="i",yaxs="i",xlab="年份",lty=1,lwd=2,type="o", ylab="出庫批次", xlim=c(2012,2020),ylim=c(0,8)) # 對8 年的時間序列數(shù)據(jù)作圖，見圖2；axis(1,c(2012,2013,2014,2015,2016,2017,2018,2019),labels=c("2012","2013","2014","2015","2016","2017","2018", "2019")) # 圖2 橫坐標(biāo)分設(shè)別為2012—2019 年；abline(h=mean(y1),lty=2,lwd=2) #求所有出庫批次平均值并將數(shù)值加載圖上，見圖2；grid（32,8,lty=3,col="black"）#添加網(wǎng)格線；legend（2012,8,c（"出庫批次","出庫批次平均值"），lwd=2，lty=c（1，2））#添加圖例。

3.1.2 繪制歷年相同月份數(shù)據(jù)對比圖 以12 個月份為橫坐標(biāo)，以出庫批次為縱坐標(biāo)，將2012—2019 年的數(shù)據(jù)繪制于圖3，每條線代表某一年的出庫批次變化。該藥品呈現(xiàn)明顯的季節(jié)性變化，每年藥品出庫幾乎都有兩個高峰，說明該藥品需求受到季節(jié)的影響?？梢詫υ摃r間序列進(jìn)行進(jìn)一步的季節(jié)分解。具體見圖3。

圖3 歷年相同月份數(shù)據(jù)對比圖

代碼如下：

year2012=yaopin[1：12,] #將數(shù)據(jù)框yaopin 的前12個數(shù)據(jù)作為2012 年的出庫批次數(shù)據(jù)賦值給year2012；

year2013 = yaopin [13：24,]；year2014 = yaopin[25：36,];year2015=yaopin[37：48,];year2016=yaopin[49：60,];year2017 = yaopin [61：72,]; year2018 = yaopin[73：84,];year2019=yaopin [85：96,]#分別將yaopin 其余數(shù)據(jù)作為2013—2019 年數(shù)據(jù)賦值給向量year2013-year2019

plot(year2012,pch=5,ylim=c(0,7),xlim=c(0,13),xaxs="i",yaxs="i",xaxt="n", xlab="月份", ylab="出庫批次") #將2012 年的12 個月數(shù)據(jù)描點(diǎn)繪圖；

points (year2013,pch=1);points (year2014,pch=4);points(year2015,pch=25);points(year2016,pch=16);points(year2017,pch =15);points (year2018,pch =11);points(year2019,pch=12) #分別將2013—2019 年的數(shù)據(jù)用不同點(diǎn)繪圖；

lines(year2012,lty=1); lines(year2013,lty=2); lines(year2014,lty=3);lines(year2015,lty=4);lines(year2016,lty=1);lines(year2017,lty=6);lines(year2018,lty=7);lines(year2019,lty=8)#分別將2012—2019 年的數(shù)據(jù)用線段繪圖

axis(1, c(1：12), labels=c("1 月","2 月","3 月","4 月","5 月","6 月","7 月","8 月","9 月","10 月","11 月","12月")) #分別用1—12 月標(biāo)記X軸刻度

grid(13,7,col="black",lty=3)#添加網(wǎng)格線

legend(0,7,c("2012","2013","2014","2015","2016","2017","2018","2019"),ncol=3,pch=c(5,1,4,25,16,15,11,12),lty=c(1,2,3,4,1,6,7,8)) #添加圖例。

3.2 檢查白噪聲 在進(jìn)行時間序列分析之前，必須先計(jì)算相關(guān)Ljung-Box 統(tǒng)計(jì)量的P值來證明該序列不是白噪聲。如果P值顯著大于顯著性水平，那么就說明該序列沒有任何統(tǒng)計(jì)規(guī)律可循，就得停止對該序列進(jìn)行統(tǒng)計(jì)分析［2］。檢查結(jié)果見表2。使用Box.test 函數(shù)對shixubiao 數(shù)據(jù)進(jìn)行上述檢查，發(fā)現(xiàn)該序列不是白噪聲，可被預(yù)測和建模。

代碼如下：

表2 時間序列的白噪聲檢驗(yàn)

Box.test (shixubiao, type="Ljung-Box",lag=6) # 檢查Ljung-Box 統(tǒng)計(jì)量，拖尾因子為6;

Box.test(shixubiao, type="Ljung-Box",lag=12) # 檢查Ljung-Box 統(tǒng)計(jì)量，拖尾因子為12。

3.3 季節(jié)性周期考查

3.3.1 季節(jié)性周期分解與繪圖 使用decompose 函數(shù)對該時間序列的季節(jié)性周期進(jìn)行分解，可以得到季節(jié)因子、趨勢因子和隨機(jī)因子，對這些因子數(shù)據(jù)進(jìn)行作圖，見圖4。匯總季節(jié)因子數(shù)據(jù)見表3。

表3 季節(jié)因子數(shù)據(jù)表

圖4 季節(jié)分解圖

代碼如下：

jijie＜-decompose(shixubiao, type="multiplicative")# 將向量shixubiao 進(jìn)行季節(jié)性分解并將結(jié)果賦值給jijie 的列表；

plot(jijie＄seasonal,xaxs="i",yaxs="i",xlab="年份",lwd=2,ylab=" 季節(jié)因子",xaxt="n",main="A 季節(jié)效應(yīng)圖",xlim=c(2012,2020),ylim=c(0,2)) # 將向量jijie 的季節(jié)因子數(shù)據(jù)繪圖，見圖4-A；

axis(1,c(2012,2013,2014,2015,2016,2017,2018,2019,2020),labels=c("2012","2013","2014","2015","2016","2017","2018","2019","2020")) #將橫坐標(biāo)刻度分別用2012—2019 表示，見圖4-A；

grid(32,6,lty=3,col="black")#添加網(wǎng)格線。見圖4-A。

plot(jijie＄trend,xaxs="i", yaxs="i",xlab="年份",lwd=2,ylab="趨勢因子", xaxt="n", main="B 趨勢圖",xlim=c(2012,2020),ylim=c(0,3)) #將向量jijie 的趨勢因子數(shù)據(jù)繪圖，見圖4-B

axis(1,c(2012,2013,2014,2015,2016,2017,2018,2019,2020),labels=c("2012","2013","2014","2015","2016","2017","2018","2019","2020")) # 將橫坐標(biāo)刻度分別用2012—2019 表示，見圖4-B；

grid(32,6,lty=3,col="black")#添加網(wǎng)格線，見圖4-B。

plot(jijie ＄random,xaxs="i",yaxs="i",xlab="年份",lwd=2,ylab=" 隨機(jī)因子",xaxt="n", main="C 隨機(jī)波動圖",xlim=c(2012,2020),ylim=c(0,2)) # 將列表jijie 的隨機(jī)因子數(shù)據(jù)繪圖，見圖4-C；

axis(1,c(2012,2013,2014,2015,2016,2017,2018,2019,2020),labels =c ("2012","2013","2014","2015","2016","2017","2018","2019","2020")) # 將橫坐標(biāo)刻度分別用2012—2019 表示，見圖4-C；

grid(32,6,lty=3,col="black")#添加網(wǎng)格線，見圖4-C。

jijie＄figure #提取列表jijie 的季節(jié)因子數(shù)據(jù)

3.3.2 季節(jié)性考查結(jié)果 圖4-A 可以看出，該時間序列呈現(xiàn)每年兩次的季節(jié)性波動。圖4-B 圖可以看出，出庫批次的發(fā)展趨勢是2012—2018 年上半年相對穩(wěn)定，從2018 年下半年開始有明顯的升高。圖4-C 圖可以看出，該藥品隨機(jī)因子基本上保持在0.5～1.5 之內(nèi)。表3 可以看出，11 月的數(shù)值最高，其次是12 月和4 月，而8 月是全年最低值，這提示該藥品在每年的用藥高峰出現(xiàn)在11—12 月，其次是4 月。8 月份為全年用藥低谷。

3.4 自動預(yù)測 對上述數(shù)據(jù)進(jìn)行5 種模型預(yù)測，預(yù)測時間為6 個月，置信區(qū)間為95%。方法：使用forecast 軟件包里的naive、snaive、hw 函數(shù)分別對向量shuju 進(jìn)行自動建模同時并預(yù)測，預(yù)測后的數(shù)據(jù)分別賦值給yuce_naive、yuce_snaive、yuce_hw 的列表中。使用ets 和auto.arima 函數(shù)分別對向量shuju 進(jìn)行自動建模，再用forecast 函數(shù)對建模后的數(shù)據(jù)進(jìn)行預(yù)測，預(yù)測后的數(shù)據(jù)分別賦值給列表yuce_ets 和yuce_arima。

代碼如下：

yuce_naive＜-naive(shuju,h=6 ,level = c(95))

yuce_snaive＜-snaive(shuju,h=6,level = c(95))

yuce_hw＜-hw(shuju,h=6,seasonal= "multiplicative",level = c(95))

yuce_ets＜- forecast(ets(shuju),h=6,level = c(95))

yuce_arima＜- forecast(auto.arima(shuju),h=6,level =c(95))

3.5 模型檢查 時序擬合某一模型之后，需要檢查平穩(wěn)性，不平穩(wěn)的模型是不能用來預(yù)測的。如果該模型是適用的，那么擬合殘差應(yīng)該是一白噪聲序列［3］。為此，可以使用shapiro.test 函數(shù)分別對上述5 種函數(shù)建模結(jié)果的殘差值＄residuals 列元素進(jìn)行W檢驗(yàn)，從而判斷是否為白噪聲序列。通過W檢驗(yàn)發(fā)現(xiàn)，snaive、auto.arima 這2 個函數(shù)建立的模型殘差P＜0.05，說明模型中有用的信息沒有被充分提取，予以舍棄。W檢驗(yàn)結(jié)果見表4。

表4 五種函數(shù)建模結(jié)果的白噪聲檢查對照表

代碼如下：

shapiro.test (yuce_naive ＄residuals)

shapiro.test (yuce_snaive ＄residuals)

shapiro.test (yuce_hw ＄residuals)

shapiro.test (yuce_ets ＄residuals)

shapiro.test (yuce_arima ＄residuals)

3.6 模型篩選

3.6.1 擬合優(yōu)度對比 使用MAPE 函數(shù)對剩余3 個函數(shù)的預(yù)測平均值＄mean 列元素分別與向量shijizhi 里的實(shí)際值進(jìn)行擬合優(yōu)度對比（MAPE 值越小，說明建模的誤差越小）。返回的yuce_naive、yuce_hw、yuce_ets 的MAPE 值分別為0.549 6、0.256 5、0.376 8。

代碼如下：

mape(yuce_naive ＄mean,shijizhi)

mape(yuce_hw ＄mean,shijizhi)

mape(yuce_ets ＄mean,shijizhi)

3.6.2 三種函數(shù)預(yù)測誤差的考查 分別提取上述剩余3 個函數(shù)預(yù)測結(jié)果的＄mean 列元素，返回各個函數(shù)預(yù)測的6個月預(yù)測平均值。將該值與2019 年9—12 月實(shí)際數(shù)據(jù)值匯總并逐個進(jìn)行比對。見表5。

表5 三種函數(shù)建模所得到的近6 個月預(yù)測平均值和實(shí)際值與誤差對照表

代碼如下：

yuce_naive＄mean

yuce_hw＄mean

yuce_ets＄mean

3.6.3 建模函數(shù)篩選結(jié)果 根據(jù)“3.6.1”項(xiàng)下結(jié)果，hw函數(shù)的預(yù)測平均值與實(shí)際值誤差的MAPE 值最小。而表5 中的3 種函數(shù)預(yù)測平均值與實(shí)際值的詳細(xì)對比也證實(shí)了hw 函數(shù)建模的誤差最小。2019 年7—12 月hw函數(shù)預(yù)測平均值為2.4 個批次，實(shí)際平均值3.2 個批次，兩者平均誤差-12.9%。因此，本次研究選擇hw 函數(shù)模型作為擬合模型。

提取列表yuce_hw 里的＄model 列元素，得到hw函數(shù)自動擬合模型的相關(guān)參數(shù)。

代碼如下：

yuce_hw＄model#提取列表yuce_hw 里的＄model 列元素；

該平滑參數(shù)為：α= 0.584 8，β=3×10-4，γ=1×10-4。

3.7 對篩選后的模型數(shù)據(jù)與實(shí)際數(shù)據(jù)進(jìn)行繪圖 使用plot 函數(shù)對該hw 函數(shù)建擬合與預(yù)測結(jié)果繪圖，觀察建模擬合值、預(yù)測平均值與實(shí)際值的直觀效果。結(jié)果顯示模型的擬合值、預(yù)測平均值圖形與實(shí)際值圖形比較接近。模型預(yù)測平均值與實(shí)際值雖然有一定誤差，但也在預(yù)測區(qū)間的上限和下限之內(nèi)。見圖5。

代碼如下：

plot(shixubiao,lwd=2,pch=19,type="b",lty=1，xaxs="i",yaxs="i",ylab="出庫批次",xlab="年份",main="",xlim=c(2012,2020),ylim = c(0,9)) #將8 年歷史實(shí)際值用實(shí)線和實(shí)心圓點(diǎn)表示

lines(yuce_hw ＄fitted,lwd=2,lty=2,pch=1,type="b") #將hw 函數(shù)擬合的7 年半模型數(shù)據(jù)用虛線和空心圓圈繪圖

lines (yuce_hw ＄mean,lwd=2,lty=2,pch=6,type="b")#將hw 函數(shù)預(yù)測的6 個月平均值用虛線和空心倒三角形點(diǎn)繪圖

U = yuce_hw＄upper ; L = yuce_hw ＄lower # 將hw 函數(shù)中的預(yù)測區(qū)間上限和下限數(shù)據(jù)分別賦值給U和L

points(U,pch=0); points(L,pch=0) #將hw 函數(shù)預(yù)測區(qū)間的上限和下限用空心方形點(diǎn)繪圖，見圖5；

lines(U, lwd=2,lty=4); lines(L, lwd=2,lty=4) # 將預(yù)測區(qū)間的上限和下限的數(shù)據(jù)用點(diǎn)劃線繪圖；axis(1,c(2013,2014,2015,2016,2017,2018,2019),labels=c("2013","2014","2015","2016","2017","2018","2019")) # 將橫坐標(biāo)中間的刻度分別用2013-2019 表示

grid(32,9,col="black",lty=3) #添加網(wǎng)格線

legend(2012,8.5,c("實(shí)際值","模型擬合值","模型預(yù)測平均值"," 預(yù)測區(qū)間的上限和下限"),pch=c(19,1,6,0),lty=c(1,2,2,4)) #添加圖例。

圖5 hw 函數(shù)擬合預(yù)測結(jié)果與實(shí)際結(jié)果對照圖

4 討論

4.1 數(shù)據(jù)預(yù)處理 由于藥品生產(chǎn)總是以批次為單位的，生產(chǎn)管理者更關(guān)注未來應(yīng)當(dāng)生產(chǎn)多少個批次的藥品而不是多少支藥，因此本次研究在對藥品出庫數(shù)據(jù)建模時使用了出庫量與每批藥品理論產(chǎn)量之比。例如，假設(shè)每批滴鼻液生產(chǎn)的理論數(shù)量是2 000 支，而當(dāng)月制劑室出庫了400 支給藥房，那么當(dāng)月的出庫批次就是400/2 000=0.25 個批次。

4.2 建模函數(shù)的選擇 雖然auto.arima 函數(shù)和ets 函數(shù)分別包含了naive、snaive、hw 函數(shù)自動建模的功能，但是其自動篩選的模型都有局限性。自動建模得到的結(jié)果并不一定是“最優(yōu)”，理想的模型仍需研究者自行比較和判斷［4］。

4.3 建模函數(shù)的預(yù)測誤差 本次研究結(jié)果表明藥品出庫批次的實(shí)際值與預(yù)測平均值存在一定的誤差，這些誤差從-28.9%到20.8%之間。雖然這個誤差范圍較大，但是對于生產(chǎn)管理和決策者來說，已經(jīng)足可以避免斷藥或藥品生產(chǎn)過剩。

4.4 本藥品的銷售特征 小兒鹽酸麻黃堿滴鼻液主要用于過敏性鼻炎和上呼吸道感染引起的鼻塞。在冬季的11 ～12 月，上呼吸道感染患者大量增加。而在每年春季的4 月，過敏性鼻炎患者增多，因而該藥品銷量也相應(yīng)增多，在季節(jié)效應(yīng)圖上反映出了每年的這兩個季節(jié)性波動。這提示管理者可以根據(jù)研究結(jié)果在高峰或低谷到來之前提前做好藥品生產(chǎn)計(jì)劃。

4.5 局限性 本次研究所采集的數(shù)據(jù)存在一定的局限性。首先，時間序列分析的研究前提是事件的發(fā)展通常具有一定的慣性［5］。但是就藥品銷售而言，這種需求慣性受到某種不確定性因素的限制，即一旦某種高致病性傳染病突然在社會上大面積流行傳播并且難以控制，就會影響相關(guān)藥品銷售數(shù)據(jù)預(yù)測的準(zhǔn)確性。所以，雖然本次研究選擇了hw 函數(shù)用于建模，但是一旦將來出庫數(shù)據(jù)發(fā)生了變化，那么建模結(jié)果未必還是hw 函數(shù)。其次，本次研究受采集條件所限，將院內(nèi)藥房從制劑室領(lǐng)取的藥品出庫數(shù)據(jù)作為分析基礎(chǔ)，但這只能間接反映藥品銷售，只有藥房銷售給患者的數(shù)據(jù)才是最直接的反映。因?yàn)樗幏款I(lǐng)取的院內(nèi)制劑總要積壓一部分以備不時之需，這種人為占壓庫存的因素可能導(dǎo)致采集到的數(shù)據(jù)未必反映市場實(shí)時需求。下一步的研究是要采集藥房的每月銷售數(shù)據(jù)進(jìn)行預(yù)測，探索更佳的擬合度。

綜上所述，本次研究的重點(diǎn)不在于能否摸索到某個具體的建模函數(shù)，而是在于試圖建立起預(yù)測藥品市場需求量的方法學(xué)，從而為科學(xué)管理提供依據(jù)。雖然R語言初學(xué)較難，但是一旦建立起這種方法學(xué)就可以快速自動預(yù)測，便于日常工作。因此，本次研究的模式可以推廣到本院制劑室的其他藥品預(yù)測中去，也可為藥廠等領(lǐng)域提供參考。