周真　李翰斌　齊佳　馬德仲

摘要：藥品貨架壽命的有效預(yù)測是藥品安全管理的關(guān)鍵問題。針對傳統(tǒng)的退化過程建模不能考慮同批產(chǎn)品中個體差異的問題，提出一種有效融合先驗退化數(shù)據(jù)和現(xiàn)場退化數(shù)據(jù)的方法，可對單片藥品及新藥品進(jìn)行貨架壽命預(yù)測。依據(jù)先驗信息確定模型的參數(shù)，引入Bayes融合現(xiàn)場退化數(shù)據(jù)進(jìn)行參數(shù)更新，在此基礎(chǔ)上對單片藥品進(jìn)行貨架壽命預(yù)測。結(jié)果表明，瓶裝和片裝的五號藥片的貨架壽命預(yù)測值分別為43?95周和47?47周，與試驗值的相對誤差分別為0?043和0?051。驗證了利用融合現(xiàn)場退化數(shù)據(jù)對單片藥品及新藥品貨架壽命預(yù)測方法的可行性。

關(guān)鍵詞：

退化模型;貨架壽命;Bayes理論;藥品

DOI：10.15938/j.jhust.2019.01.009

中圖分類號： TB114?3

文獻(xiàn)標(biāo)志碼： A

文章編號： 1007-2683（2019）01-0055-05

The Degradation Process Modeling and Shelf Life Prediction of Drug

ZHOU Zhen?1，2?，LI Han?bin?1，QI Jia?1，2?，MA De?zhong?1，2

（1?School of Measurement?control Technology and Communications Engineering， Harbin University of Science and Technology， Harbin 150080，China;

2?Higher Educational Key Laboratory for Measuring and Control Technology and Instrumentations of Heilongjiang Province， Harbin 150080， China）

Abstract：The shelf life prediction of drug is the key problem of drug safety management. To solve the problem of individual difference in the same batch cannot be considered during the traditional degradation process modeling， an effective method is proposed to fuse priori degenerate data with site degraded data which can predict shelf life for the monolithic drug and new drugs?Determining the distribution of model parameters based on prior information， fusion of field data to update parameters by Bayes?On this basis， predict the shelf life of the monolithic drug?The results show that the No?5 loaded and bottled drugs shelf life predicted values are 43?95 and 47?47 weeks?The relative errors are 0?043 and 0?051?And the feasibility of fusing degradation data to predict the shelf life of the monolithic drug and new drugs is verified

Keywords：degradation model; shelf life; bayesian theory; drug

0引言

藥品質(zhì)量和效期問題是比較重大的公共衛(wèi)生安全問題，它關(guān)乎著病人的生命安全。經(jīng)統(tǒng)計，我國目前藥品不良反映案例中，有1/3是由于過期藥品引起的。藥品過期后將發(fā)生內(nèi)在質(zhì)量變化，繼續(xù)使用會引起中毒。比如曾經(jīng)轟動全國的“梅花K”事件，就是過期變質(zhì)的四環(huán)素導(dǎo)致了患者腎小管功能損害、肝損害等，嚴(yán)重危害用藥者健康。用藥的最高原則應(yīng)是“安全、有效”。藥品貨架壽命是衡量藥品安全性和有效性的重要指標(biāo)，過期的藥品對正常藥品市場、人民群眾身體健康和生命安全都有著嚴(yán)重的影響[1-2]。因此，開展對藥品貨架壽命的有效預(yù)測，為控制和防范藥品安全危害，加強藥品安全監(jiān)管提供理論上的支持。

傳統(tǒng)的壽命預(yù)測方法是基于統(tǒng)計方法得出的共性的平均壽命，難以有效處理小樣本少量藥劑的失效數(shù)據(jù);且缺乏對單個藥片壽命預(yù)測的準(zhǔn)確性[3-4]。實際上，在藥品公司研制過程中積累了大量的藥品藥效退化的信息，利用這些藥效退化信息，結(jié)合性能可靠性理論與方法，可為藥品的貨架壽命預(yù)測提供一條新的途徑[5]。

基于退化軌跡建模及分析理論是解決產(chǎn)品壽命預(yù)測等問題的重要手段。近年來，基于退化模型的壽命預(yù)測方法研究得到國內(nèi)外學(xué)者廣泛關(guān)注。Freitas等[6]基于線性退化模型對產(chǎn)品的退化過程進(jìn)行研究。相比壽命數(shù)據(jù)，退化數(shù)據(jù)能夠提供更多的可靠性信息。周月閣等[7]基于性能退化和Monte?Carlo仿真的可靠性分析方法，闡述了基本的退化建模思想及評估方法，并結(jié)合產(chǎn)品的失效數(shù)據(jù)驗證了評估的有效性及正確性。Xu等[8]針對線性性能退化模型，確定參數(shù)的先驗分布，并對產(chǎn)品中位壽命運用Bayes方法進(jìn)行統(tǒng)計推斷。另外，還有學(xué)者對可靠性模型參數(shù)的更新進(jìn)行了研究。Gebraeel等[9-10]與Elwany等[11]研究軸承第二階段的性能退化軌跡，通過指數(shù)模型與不同軸承個體振動幅值的退化趨勢相擬合，并利用現(xiàn)場數(shù)據(jù)對性能退化模型在線更新，預(yù)測了軸承的剩余壽命。Kaiser & Gebraeel[12]根據(jù)更新后的產(chǎn)品剩余壽命分布對產(chǎn)品進(jìn)行壽命預(yù)測，在此基礎(chǔ)上對產(chǎn)品的預(yù)防性維修問題進(jìn)行研究。Chen等[13]提出了一種兩階段的性能退化模型，假定性能退化路徑分段線性且分段時間點隨機(jī)，模型參數(shù)的先驗分布運用Bayes原理進(jìn)行更新。盡管已有大量文獻(xiàn)基于性能退化數(shù)據(jù)壽命預(yù)測課題進(jìn)行研究，但仍有問題需進(jìn)一步解決：①傳統(tǒng)方法是對產(chǎn)品總體的壽命預(yù)測，而無法體現(xiàn)同批產(chǎn)品個體的差異性。對于高危和昂貴藥品，利用已有的退化數(shù)據(jù)對單片藥品進(jìn)行貨架壽命預(yù)測更符合藥品安全問題的實際需求。②在新產(chǎn)品的研制過程中存在小子樣問題，如何對新產(chǎn)品少量藥劑進(jìn)行貨架壽命預(yù)測，也是所需解決的難題。

針對以上問題，本文在退化數(shù)據(jù)的基礎(chǔ)上對藥品藥效退化過程進(jìn)行建模，融合現(xiàn)場退化數(shù)據(jù)，結(jié)合Bayes理論，實現(xiàn)對單片藥品及新藥品貨架壽命預(yù)測。最后，通過某種藥品貨架壽命預(yù)測實例驗證了本文提出的利用Bayes方法融合現(xiàn)場退化數(shù)據(jù)對單片藥品進(jìn)行貨架壽命預(yù)測的合理性。

1退化量統(tǒng)計模型

1?1常用退化軌跡模型

性能退化量是指試驗樣品在某一工作條件下性能的變化量，它與測量間隔時間之間建立的函數(shù)關(guān)系為退化軌跡模型。大部分性能退化產(chǎn)品的退化量，可以用以下5種常用的退化軌跡模型進(jìn)行擬合。在誤差允許的范圍內(nèi)，可以將實際的退化軌跡被看作擬合后的退化軌跡[14-15]。

線性函數(shù)模型：

y?j=ω?jt?j+λ?j（1）

指數(shù)函數(shù)模型：

y?j=λ?j?exp?（ω?jt?j）（2）

冪函數(shù)模型：

y?j=λ?jt?ωj?j（3）

對數(shù)函數(shù)模型：

y?j=ω?j?ln?t?j+λ?j（4）

Lloyd?Lipow?模型：

y?j=λ?j-ω?j/t?j（5）

式中：j=1， 2，…，n;n為性能退化量的個數(shù);y?j為性能退化量;ω?j，λ?j為待估模型參數(shù);t?j為測量間隔時間。

1?2退化量分布法的基本原理

基于統(tǒng)計分析是退化數(shù)據(jù)分析的主要方法，步驟及實現(xiàn)方法具體如下：

1）性能退化數(shù)據(jù)（t?1，y?1j?），（t?2，y?1j?），…，（t?n，y?1j?） （j=1， 2，…，n）在t?1，t?2，…，t?n時刻對n個樣品進(jìn)行測量，失效閾值D?f根據(jù)相應(yīng)的國家標(biāo)準(zhǔn)及使用方要求來確定。

2）分析收集到的試驗數(shù)據(jù)，對各時刻所有樣本的退化量（y?i1?，y?i2?，…，y?im?，i=1，2，…，m）依次分別應(yīng)用對數(shù)正態(tài)分布、威布爾分布、正態(tài)分布進(jìn)行假設(shè)檢驗，選擇其中最優(yōu)的分布函數(shù)，并求出該分布函數(shù)中的待定參數(shù)。

3）為了確定最優(yōu)的退化軌跡模型，運用最小二乘法將步驟2）中所求的分布函數(shù)中的參數(shù)y?j分別與5種常用的退化軌跡模型進(jìn)行擬合[16]。

2單個產(chǎn)品壽命預(yù)測

2?1Bayes原理

根據(jù)?Bayes?原理（式（6）），融合未知參數(shù)的先驗信息和樣本信息，得到后驗信息，在此基礎(chǔ)上對未知參數(shù)進(jìn)行推斷。獲取單個產(chǎn)品及新品的退化數(shù)據(jù)x=[x?ij?]，運用?Bayes?原理對退化模型的參數(shù)更新，從而對單個產(chǎn)品及新品壽命預(yù)測進(jìn)行更新。

f（ω，λ|x）=g（x|μ，θ）f（μ）f??′（θ）∫D（μ）∫D（θ）g（x|μ，θ）f（μ）f??′（θ）?d?μ?d?θ（6）

其中：f（ω）、f（λ）分別表示ω和λ的先驗概率密度函數(shù);f（ω，λ|x）表示后驗概率密度函數(shù);D（ω）和D（λ）分別表示ω和λ取值范圍;g（x|μ，θ）為似然函數(shù)[17]。

2?2線性退化模型參數(shù)更新

線性退化建模中參數(shù)ω?j、λ?j的估計是關(guān)鍵，在退化模型分析中，用最小二乘法估計ω?j、λ?j的值。然而在小樣本的情況下這種估計有較大誤差，得到的估計結(jié)果也不準(zhǔn)確，難以用于產(chǎn)品的壽命預(yù)測?？梢约僭O(shè)ω?j服從ω?j～N（μ?0，σ?2?0）;λ?j服從λ?j～N（μ?1，σ?2?1）的正態(tài)分布且ω?j與λ?j相互獨立。由隨機(jī)變量的函數(shù)分布關(guān)系，可得y是服從正態(tài)分布的隨機(jī)變量，記作：

y～N（μ， σ?2）（7）

式中：μ=μ?0t?j+μ?1，σ?2=σ?2?0t?2?j+σ?2?1。

假設(shè)ω?j，λ?j的先驗分布為正態(tài)分布，選取單個試驗樣本試驗，不同試驗時間內(nèi)的試驗結(jié)果記作：

y=（y（t?1），y（t?2），y（t?3），…，y（t?j））。

y的密度函數(shù)為：

f（y）=12?π?σ?e??（y?j-μ）?22σ?2

則f（y）=12?π?（σ?2?0t?2?j+σ?2?1）?e??（y?j-μ?0t?j-μ?1）?22（σ?2?0t?2?j+σ?2?1）

似然函數(shù)為：

L（Y|ω?j，λ?j）=∏n?j=1?12?π?（σ?2?0t?2?j+σ?2?1）?e??（y?j-μ?0t?j-μ?1）?22（σ?2?0t?2?j+σ?2?1）

ω?j，λ?j的密度函數(shù)如下：

f?（ω?j）?=12?π?σ?0?e??-?（ω?j-μ?0）?22σ?2?0

f?（λ?j）?=12?π?σ?1?e??-?（λ?j-μ?1）?22σ?2?1

由?Bayes?公式可得的ω?j，λ?j后驗聯(lián)合密度函數(shù)：

f?（ω?j，λ?j|Y）?=f?（ω?j）?f?（λ?j）?L（Y|ω?j，λ?j）L（Y）（8）

即：

f?（ω?j，λ?j|Y）?=f?（ω?j）?f?（λ?j）?L（Y|ω?j，λ?j）?∫??ω?j?∫??λ?j?f（ω?j）f（λ?j）L（Y|ω?j，λ?j）?d?ω?j?d?λ?j

則?ω?j，λ?j后驗聯(lián)合密度函數(shù)：

f?（ω?j，λ?j|Y）?∝f?（ω?j）?f?（λ?j）?L（Y|ω?j，λ?j）

∝?e??-?（ω?j-μ?0）?22σ?2?0?e??-?（λ?j-μ?1）?22σ?2?1?e??-∑n?j=1?（y?j-ω?jt?j-λj）?22（σ?2?0t?2?j+σ?2?1）?（9）

又因為ω?j，λ?j的后驗分布f?（ωj，λj|Y）?是二維正態(tài)分布，假設(shè)ω?j，λ?j的相關(guān)系數(shù)為r，后驗分布均值和方差分別為（μ?0′，μ?1′）和（σ?0?′2?，σ?1?′2?），其概率密度函數(shù)如下：

f?（ω?j，λ?j|Y）?=12?π?σ?′?0σ?′?11-r?2?exp?-12（1-r?2）

（?（ω?j-μ?0）?2σ?′2?0-2r（ω?j-μ?′?0）（λ?j-μ?′?0）σ?′?0σ?′?1+?（λ?j-μ?′?0）?2σ?′2?1）

其中：

μ?′?0=（∑n?j=1?y?jt?jσ?2?0+μ?0σ?2）（σ?2+nσ?2?1）A-B

-（∑n?j=1?t?jσ?2?0）（∑n?j=1?y?jσ?2?1+μ?1σ?2）A-B（10）

μ?′?1=（∑n?j=1?y?jσ?2?1+μ?1σ?2）（σ?2+σ?2?0∑n?j=1?t?2?j）A-B

-（∑n?j=1?t?jσ?2?1）（∑n?j=1?y?jσ?2?0+μ?0σ?2）A-B（11）

σ?′2?0=σ?2σ?2?0（σ?2+nσ?2?1）A-σ?2?0σ?2?1?（∑n?j=1?t?j）?2（12）

σ?′2?1=σ?2σ?2?1（σ?2+σ?2?0∑n?j=1?t?2?j）A-σ?2?0σ?2?1?（∑n?j=1?t?j）?2（13）

r=-σ?0σ?1∑n?j=1?t?2σ?2+nσ?2?1σ?2+σ?2?0∑n?j=1?t?2?j（14）

A=（σ?2+nσ?2?1）（σ?2+σ?2?0∑n?j=1?t?2?j）（15）

B=（∑n?j=1?t?jσ?2?0）（∑n?j=1?t?jσ?2?1）（16）

2?3壽命預(yù)測

結(jié)合單個樣本信息，確定退化模型參數(shù)的后驗分布，從而預(yù)測單個產(chǎn)品壽命，即估計性能退化參數(shù)達(dá)到失效閾值D?f的時間[18]。

1）性能參數(shù)y?j≥失效閾值D?f

將壽命視為退化軌跡第一次超過閾值D?f時對應(yīng)的時間，通常稱為首達(dá)時。因此，根據(jù)首達(dá)的概念，此時對應(yīng)的壽命T可以定義為T=?inf?t|y?ij?（t）≥D?f

2）性能參數(shù)y?j≤失效閾值D?f

將壽命視為退化軌跡第一次低于閾值D?f時對應(yīng)的時間，通常稱為首達(dá)時。因此，根據(jù)首達(dá)的概念，此時對應(yīng)的壽命T可以定義為T=?inf?t|y?ij?（t）≤D?f

3實例應(yīng)用

藥品貨架壽命的定義是在生產(chǎn)和包裝后，在指定的貯藏條件下能保持其安全性并適用的期限。而藥品的保存方式可分為常溫保存和冷柜保存[19]。本文選取某種僅適用于常溫保存的藥品為研究對象，在瓶裝和片裝兩種包裝下，分別在不同批藥品中進(jìn)行抽樣，瓶裝藥品和片狀藥品均隨機(jī)抽取5片，樣品放在貨架上（溫度10℃～30℃，RH45%～75%，避光），以藥效指標(biāo)作為退化量，當(dāng)藥效指標(biāo)下降到初始值的10%時判定藥品失效，藥效指標(biāo)下降10%所經(jīng)歷的時間，即為該藥品的貨架壽命[20]。對每個的樣品重復(fù)測量6次，其中前5次每隔三周測量一次，第6次測量是在第5次測量的6周后測試，測量的退化數(shù)據(jù)為?y?ij??（1）?與y?ij??（2）?，如表1所示。

分別建立藥品藥效退化模型，將藥效退化數(shù)據(jù)與5種模型（式（1）～式（5））運用最小二乘法擬合，擬合參數(shù)如表2所示，且由表2可知瓶裝和片裝藥品藥效退化數(shù)據(jù)與線性退化模型擬合程度最高，均方根誤差分別為0?364和0?201。即?y?ij?=ω?j?t?j?+λ?j?是瓶裝和片裝條件下該藥品藥效最合適的退化模型。

將表1數(shù)據(jù)和一個新的樣品退化數(shù)據(jù)（見表3）融合。運用Bayes方法更新藥品退化軌跡模型參數(shù)，按照式（7）～式（16）分別算出更新后退化模型系數(shù)?ω?j?（1）?，λ?j?（1）?的平均值和方差，ω?j?（2）?，λ?j?（2）?的平均值和方差，如表4所示。

監(jiān)測五號藥品失效前的時間為貨架壽命實驗值。片裝和瓶裝五號藥品的實驗貨架壽命分別為45?92周和45?17周[20]。二階Taylor近似估計方法是預(yù)測藥品貨架壽命傳統(tǒng)的方法，預(yù)測值為同批藥品貨架壽命的均值，具有參考價值。表5給出本文方法與二階Taylor近似估計方法對于貨架壽命預(yù)測的結(jié)果。通過與實驗值比較可得出預(yù)測相對誤差，如表6所示。分析可知，二階Taylor近似估計方法與實驗值的差距明顯。通過文中方法對五號藥品藥效退化參數(shù)進(jìn)行更新，藥品貨架壽命的預(yù)測值接近實驗值。預(yù)測貨架壽命精度更高，相對誤差更小。

4結(jié)語

本文提出一種基于退化數(shù)據(jù)的藥品貨架壽命預(yù)測方法。通過對藥效退化數(shù)據(jù)的分析，利用Bayes理論融合先驗信息和現(xiàn)場退化數(shù)據(jù)，對藥效退化數(shù)據(jù)和退化軌跡模型進(jìn)行了更新。實現(xiàn)了對單片藥品貨架壽命的預(yù)測，也可對研發(fā)過程中的新藥品、昂貴藥品的貨架壽命預(yù)測提供理論指導(dǎo)。該方法擺脫了傳統(tǒng)統(tǒng)計方法對大樣本統(tǒng)計信息的依賴，實現(xiàn)了對單個產(chǎn)品壽命的預(yù)測，為小樣本條件下產(chǎn)品壽命的預(yù)測提供了有效手段，具有工程應(yīng)用價值。

參 考 文 獻(xiàn)：

[1]胡穎廉?監(jiān)管和市場：我國藥品安全的現(xiàn)狀、挑戰(zhàn)及對策[J]. 中國衛(wèi)生政策研究，2013（7）：38

[2]林麗花?強化病區(qū)高危藥品安全管理的效果觀察[J]. 中國藥物經(jīng)濟(jì)學(xué)，2016，11（8）：191

[3]郭建英，孫永全，于曉洋?可靠性增長技術(shù)發(fā)展動態(tài)詮釋[J]. 哈爾濱理工大學(xué)學(xué)報，2011，16（2）：1

[4]朱顯輝，師楠，田德文，等?PWM驅(qū)動電機(jī)可靠性研究發(fā)展及展望[J]. 哈爾濱理工大學(xué)學(xué)報， 2014，19（4）：22

[5]胡姚剛，李輝，廖興林，等?風(fēng)電軸承性能退化建模及其實時剩余壽命預(yù)測[J]. 中國電機(jī)工程學(xué)報，2016，20（6）：1643

[6]FREITAS M A，DE TOLEDO M LG，Colosimo E A?Using Degradation Data to Assess Reliability：A Case Study on Train Wheel Degradation[J]. Quality and Reliability Engineering International，2009，25（5）：607

[7]周月閣，葉雪榮，翟國富?基于性能退化和Monte?Carlo仿真的系統(tǒng)性能可靠性評估[J]. 儀器儀表學(xué)報，2014（5）：1185

[8]XU A?，TANG Y?Objective Bayesian Analysis for Linear Degradationmodels[J]. Communications in Statistics?Theory and Methods，2012，41（21）：4034

[9]GEBRAEEL N Z，LAWLEY M A，Rong L，Jennifer K?Residual?Life Distributions from Component Degradation Signals：A Bayesian Approach[J].IIE Transactions， 2005， 37（6）： 543

[10]GEBRAEEL N?Sesory?updated Residual Life Distributions for Components with Exponential Degradation Patterns[J]. IEEE Transactions on Automation Science and Engineering，2006，3（4）：382

[11]ELWANY A?，GEBRAEEL N?Real?time Estimation of Mean Remaining Life Using Sensor?based Degradation Models[J].Journal of Manufacturing Science and Engineering，2009，131（5）：51

[12]KAISER K A，GEBRAEEL N Z?Predictive Maintenance Management Using Sensor Degration Models[J].IEEE Transaction on Systems，Man and Cybernetics，Part A：Systems and Humans，2009，39（4）：840

[13]CHEN N?，TSUI K?L?Condition Monitoring and Remaining Useful Life Prediction Using Degration Signals：Revisited[J]. IIE Transactions，2013，45（9）：939

[14]邱榮華， 巨孔亮， 董友耕，等?基于性能退化的小子樣電主軸可靠性試驗研究[J]. 中國機(jī)械工程， 2016， 27（20）：2738

[15]晁代宏，馬靜，陳淑英?應(yīng)用多元性能退化量評估光纖陀螺貯存的可靠性[J]. 光學(xué)精密工程，2011（1）：35

[16]蔡忠義， 陳云翔， 張諍敏，等?基于比例失效率退化模型的可靠性評估方法[J]. 系統(tǒng)工程與電子技術(shù)， 2015， 37（8）：1943

[17]王華偉，高軍，吳海橋?基于Bayes模型平均的航空發(fā)動機(jī)可靠性分析[J]. 航空動力學(xué)報，2014（2）：305

[18]彭寶華，周經(jīng)倫，潘正強?Wiener過程性能退化產(chǎn)品可靠性評估的Bayes方法[J]. 系統(tǒng)工程理論與實踐， 2010（3）： 543

[19]王笑琴?常溫保存的藥品存放于冷藏處對藥效的影響[J]. 醫(yī)療裝備，2016（18）：112

[20]SHAO J?and CHOW S?C?Constructing Release Targets for Drug Products：A Bayesian Decision Theory Approach[J]. Journal of the Royal Statistical Society， 1991，40（3）：381