以小兒鹽酸麻黃堿滴鼻液為例探索用R 語言預測兒童用藥需求

田 剛

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

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

1 研究工具

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

2 研究方法

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

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

表1 源數據

代碼和解釋如下：

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”的數據框；

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

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

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

圖1 y.csv 文件的內容

3 研究步驟

3.1 數據初步考察

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

圖2 每月出庫批數時序圖

代碼如下：

y1＜- as.numeric（unlist（yaopin））#轉換為數值型向量并賦值給y1；

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

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

圖3 歷年相同月份數據對比圖

代碼如下：

year2012=yaopin[1：12,] #將數據框yaopin 的前12個數據作為2012 年的出庫批次數據賦值給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 其余數據作為2013—2019 年數據賦值給向量year2013-year2019

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

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 年的數據用不同點繪圖；

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 年的數據用線段繪圖

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

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

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

代碼如下：

表2 時間序列的白噪聲檢驗

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

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

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

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

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

圖4 季節(jié)分解圖

代碼如下：

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

plot(jijie＄seasonal,xaxs="i",yaxs="i",xlab="年份",lwd=2,ylab=" 季節(jié)因子",xaxt="n",main="A 季節(jié)效應圖",xlim=c(2012,2020),ylim=c(0,2)) # 將向量jijie 的季節(jié)因子數據繪圖，見圖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")) #將橫坐標刻度分別用2012—2019 表示，見圖4-A；

grid(32,6,lty=3,col="black")#添加網格線。見圖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 的趨勢因子數據繪圖，見圖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")) # 將橫坐標刻度分別用2012—2019 表示，見圖4-B；

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

plot(jijie ＄random,xaxs="i",yaxs="i",xlab="年份",lwd=2,ylab=" 隨機因子",xaxt="n", main="C 隨機波動圖",xlim=c(2012,2020),ylim=c(0,2)) # 將列表jijie 的隨機因子數據繪圖，見圖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")) # 將橫坐標刻度分別用2012—2019 表示，見圖4-C；

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

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

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

3.4 自動預測 對上述數據進行5 種模型預測，預測時間為6 個月，置信區(qū)間為95%。方法：使用forecast 軟件包里的naive、snaive、hw 函數分別對向量shuju 進行自動建模同時并預測，預測后的數據分別賦值給yuce_naive、yuce_snaive、yuce_hw 的列表中。使用ets 和auto.arima 函數分別對向量shuju 進行自動建模，再用forecast 函數對建模后的數據進行預測，預測后的數據分別賦值給列表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)的模型是不能用來預測的。如果該模型是適用的，那么擬合殘差應該是一白噪聲序列［3］。為此，可以使用shapiro.test 函數分別對上述5 種函數建模結果的殘差值＄residuals 列元素進行W檢驗，從而判斷是否為白噪聲序列。通過W檢驗發(fā)現，snaive、auto.arima 這2 個函數建立的模型殘差P＜0.05，說明模型中有用的信息沒有被充分提取，予以舍棄。W檢驗結果見表4。

表4 五種函數建模結果的白噪聲檢查對照表

代碼如下：

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 函數對剩余3 個函數的預測平均值＄mean 列元素分別與向量shijizhi 里的實際值進行擬合優(yōu)度對比（MAPE 值越小，說明建模的誤差越?。７祷氐膟uce_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 三種函數預測誤差的考查 分別提取上述剩余3 個函數預測結果的＄mean 列元素，返回各個函數預測的6個月預測平均值。將該值與2019 年9—12 月實際數據值匯總并逐個進行比對。見表5。

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

代碼如下：

yuce_naive＄mean

yuce_hw＄mean

yuce_ets＄mean

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

提取列表yuce_hw 里的＄model 列元素，得到hw函數自動擬合模型的相關參數。

代碼如下：

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

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

3.7 對篩選后的模型數據與實際數據進行繪圖 使用plot 函數對該hw 函數建擬合與預測結果繪圖，觀察建模擬合值、預測平均值與實際值的直觀效果。結果顯示模型的擬合值、預測平均值圖形與實際值圖形比較接近。模型預測平均值與實際值雖然有一定誤差，但也在預測區(qū)間的上限和下限之內。見圖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 年歷史實際值用實線和實心圓點表示

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

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

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

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

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

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

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

圖5 hw 函數擬合預測結果與實際結果對照圖

4 討論

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

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

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

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

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

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