時間序列模型對細菌耐藥率變化趨勢預測價值

錢鑫，熊世娟，楊廷秀，胡方芳，謝娟

貴州省人民醫院，貴州 貴陽 550002

隨著抗菌藥物的廣泛暴露，病原微生物的不斷進化，細菌耐藥已成為全球公共健康領域的重大挑戰。近年來，隨著碳青霉烯類耐藥肺炎克雷伯菌（carbapenem-resistant klebsiella pneumoniae，CRKP）檢出率呈現明顯上升趨勢，已嚴重威脅人類健康，給臨床治療帶來巨大困難與挑戰。在大數據時代，如何利用所收集的數據為決策提供服務已經逐漸成為這個時代的主題。美國學者 Box 和英國統計學專家 Jenkins 于1976 年提出的自回歸移動平均模型（autoregressive moving average model，ARMA）［1］，目前在臨床流行病學方面應用較為廣泛。本研究旨在探討ARMA模型對CRKP耐藥率趨勢的預測價值，為今后醫院感染的預防和控制，以及進一步加強抗菌藥物管理提供相應的數據支持和理論參考。

1 資料與方法

1.1 研究資料

本次研究的資料來自貴州省人民醫院在2013年1 月至2018 年12 月送檢的病原學檢測結果，去除了同一個患者分離的重復菌株。標本均由各病區護理人員按標準采集，采集流程按照《全國臨床檢驗操作規程》［2］執行。采集的標本均按時送至微生物室，并行標本質量的人工復檢，合格的標本在初步培養的基礎上，采用法國梅里埃公司 VITEK 2 Compact 全自動微生物鑒定系統進行細菌鑒定及藥敏試驗，所分離的CRKP 菌株根據美國臨床與實驗室標準協會（clinical and laboratory standards institute，CLSI）標準對藥敏結果進行判定。

1.2 模型介紹

ARMA 模型是時間序列預測分析方法之一［3］，一般常用經差分的ARIMA（p，d，q）模型。AR 是自回歸模型，p 為自回歸系數；MA 為移動平均模型，q為移動平均系數，d 為成為平穩序列所需差分的階數。

1.3 模型建立

建模過程主要分3 個步驟進行：（1）確定模型的參數，時間序列模型分析需要保證時間序列的平穩性。如果時間序列表現得不夠平穩，則通過差分法可以使其均數與方差趨于平穩化，之后觀察所作自相關函數圖（autocorrelation function，ACF）和偏自相關函數圖（partial correlation function，PACF）的截尾、拖尾情況，初步進行定階，按低階到高階的順序，提供相關的粗放性模型進行篩選。（2）檢驗模型的參數，在模型擬合的過程中，如果存在多個模型均通過Box-Ljung 檢驗，則提示其殘差序列均為白噪聲序列［4］，根據赤池信息準則（akaike information criterion，AIC）［5］、貝葉斯信息準則（bayesian information criterion，BIC）［6］及許瓦茲貝葉斯準則（schwarz bayesianc riterion，SBC）［7］來選擇最優模型。若模型參數無統計學意義，那么分析流程則中止，一般需要返回模型識別階段重新定階，再進行模型參數的估計和診斷。（3）評估模型預測價值：比較擬選模型的預測值與實際值，對擬選模型間預測的準確性進行評價。

1.4 數據分析

以月為單位，將耐藥率數據錄入Excel 2016中，分別導入IBM SPSS Statistics 19.0 和SAS Statistical Discovery 13，將需要分析的原始數據進行平穩化處理后，采用時間序列模型進行擬合分析，以最終確定的模型來預測2019 年的細菌耐藥率，并采用絕對誤差和相對誤差對預測值進行評價，時間序列的平穩性運用EViews10.0，通過擴展迪基-福勒檢驗（augmented dickey-fuller test，ADF）［8］對單位根進行檢驗。

2 結果

2.1 全球耐藥趨勢及醫院CRKP 耐藥形勢

據筆者調研，美國非營利組織“動態疾病、經濟和政策研究中心”（center for disease dynamics,economics &policy，CDDEP）對科學文獻以及國家和地區監測系統的分析數據提示，中國肺炎克雷伯菌的碳青霉烯類耐藥形勢已較為嚴峻，見圖1～3。結合院細菌耐藥監測情況，近6 年院CRKP 耐藥率呈現持續上升趨勢，2018 年院美羅培南耐藥率為22.62%，稍低于全國年度平均水平28.60%，但該菌2018 年在院的碳青霉烯耐藥率較之2017 年出現了大幅度提高，引起了臨床的高度關注，特別是碳青霉烯類使用頻率較高的重點科室，如重癥監護室、血液科等。MDRO千日檢出率近6 年來總體變化趨勢不大，雖然MDRPA 和MRSA 的千日檢出率亦相對穩定，2018 年分別為0.15‰和0.31‰，但CRE 千日檢出率出現大幅度抬頭趨勢，臨床面臨著較大的壓力與挑戰。

圖1 肺炎克雷伯菌對碳青霉烯類藥物的耐藥性

圖2 肺炎克雷伯菌的細菌耐藥率

2.2 模型處理與檢驗

將耐藥率原始數據作時序圖，可觀察到時間序列的波動比較大，并不能滿足時間序列分析的平穩化要求，見圖4。該時間序列經ADF 檢驗，P=0.482，P＞0.05，提示序列存在單位根。對該序列進行一階差分后序列圖顯示基本達到平穩，經ADF檢驗P＜0.05，提示序列無單位根，符合序列平穩化要求，見圖5。通過譜密度查看周期性變化，顯示序列的季節性波動不明顯，因此該時間序列可選擇ARIMA（p，1，q）模型，初步篩選ARIMA（1，1，1）、IMA（1，1）、ARI（1，1）三個時間序列模型。初選模型均通過了Ljung-Box Q 檢驗，提示上述模型觀測值是相互獨立的，殘差是白噪聲序列，各備選模型參數信息，見表1。通過方差、AIC、SBC、R 方及均方根誤差（root mean square error，RMSE）和平均絕對誤差（mean absolute error，MAE）評估模型擬合優度，顯示ARIMA（1，1，1）模型權重較高，其Ljung-Box Q 統計量Q=21.348，P=0.166，模型參數估計結果P＞0.05，其殘差ACF 和PACF 圖，見圖6。

圖3 MDRO千日檢出率變化

圖4 原始數據時序圖

圖5 一階差分時序圖

圖6 模型ARIMA（1，1，1）的殘差ACF和PACF圖

2.3 模型的預測

盡管ARIMA（1，1，1）模型參數權重更高，但與IMA（1，1）相比殘差無明顯差異，見圖7。選取2018 年1 至12 月的真實值分別與兩個模型的預測值進行比較，計算模型預測值的相對誤差和平均相對誤差。結果顯示，IMA（1,1）模型的平均相對誤差相對較小，IMA（1，1）的擬合效果更接近于真實值，與可視化觀察結果一致。基于上述的分析流程和推斷，IMA（1，1）模型的擬合度最優。

圖7 ARIMA（1，1，1）和IMA（1，1）模型的擬合曲線

3 討論

近年來，多重耐藥菌導致的感染是造成醫療保健相關發病率、死亡率和醫療成本支出的主要來源［9-11］。特別是肺炎克雷伯菌作為腸桿科兼性厭氧菌，是醫院感染的重要病原菌［12］。中國的細菌耐藥監測網CHINET 顯示CRKP 檢出率逐年上升［13］。細菌耐藥問題已經被G20 峰會和聯合國大會正式納入政治議題。2017 年3 月國衛辦醫發（2017）10號文出臺，進一步強化了碳青霉烯類抗菌藥物的管理。因此，制定有效的戰略可以對這類病原菌實施監測和預測，同時可以有助于采取有效地干預措施。本研究數學模型的識別方式，將臨床收集到以時間為序的細菌耐藥數據進行未來值探索與預測，可為我們的公共衛生部門提供重要的決策信息。目前，時間序列模型已廣泛應用于患病率的預測，但基于細菌耐藥率的預測研究卻不多。通過本研究可以看到，貴州省人民醫院CRKP 檢出率近年大幅增高趨勢，預測提示存在持續走高的風險。結合衛健委發布的《碳青霉烯類抗菌藥物臨床應用評價細則》以及醫院實際應用情況分析，需要密切關注以下問題：（1）手衛生不到位使耐藥菌株病區內或病區間傳播，必然導致碳青霉烯類用量增加，出現耐藥菌篩選壓力；（2）部分臨床醫生對經驗治療把握不到位，導致療程過長，或者對定植菌未仔細評估；（3）不少醫務人員存在重癥就是耐藥菌感染的錯誤認識，以及碳青霉烯類使用頻次不足等問題；（4）高級職稱的臨床醫生對抗菌藥物培訓參與度不夠，持續學習不足，對碳青霉烯類藥物應用指征了解不充分；（5）醫院抗菌藥物TDM 監測手段不完善，對感染性疾病的療效評估缺失個體化評估手段，一定程度上影響了臨床用藥的決策。

表1 備選模型的估計參數

表2 2018年1至12月耐藥率觀測值與預測值的比較

本研究通過原始數據的參數估計和檢驗，獲得了模型IMA（1，1），應用該模型預測的結果顯示，擬合值的動態趨勢與真實值的升降規律較為相似，提示時間序列模型可較好的分析和預測細菌耐藥率的變化趨勢。不過，該模型雖然具有不受數據類型束縛、短期預測良好的特點，但考慮到細菌耐藥規律的復雜性，模型的長期應用需要保證有足夠的數據，并及時補充新的數據，同時對模型類型和參數不斷完善修正，才可能確保預測值的精準度，便于“早發現、早干預”，為進一步加強抗菌藥物的管理，有效控制細菌耐藥作出信息預警。