基于多重填補的廣義線性模型在腎臟疾病研究中的應用

王 威 綜述 楊 帆 審校

在基于臨床數(shù)據(jù)的研究中,無論是回顧性還是前瞻性的研究,常常會出現(xiàn)數(shù)據(jù)缺失的情況。當出現(xiàn)了缺失值,我們應給予妥當?shù)奶幚?以便做出可靠的統(tǒng)計推斷。目前,應對缺失值的方法是多種多樣的,主要有刪除法和填補法。刪除法又稱完整案例分析(complete case analysis)或列表刪除(listwise deletion),即刪除存在缺失值的所有觀測,對保留下來的無缺失值的觀測進行統(tǒng)計分析。此方法的操作是最簡單易行的,因此為大多數(shù)研究者所采用,也是多數(shù)統(tǒng)計分析軟件所默認采用的方法。然而,此方法在實際運用中可能會存在兩點不足:一是缺失值的存在可能不是隨機的,即存在缺失值的觀測與完整觀測之間存在某些方面的差異,且這些差異會給后續(xù)的參數(shù)估計帶來偏倚;二是在研究變量數(shù)較多的情況下,所有研究變量均無缺失的觀測可能會比較少,運用此法會舍棄過多觀測,不僅會嚴重降低投入產(chǎn)出比,而且會增大參數(shù)估計的標準差和置信區(qū)間,降低統(tǒng)計功效。只有在存在缺失值的觀測占比很小的時候(比如<5%),所造成的參數(shù)估計的偏倚和統(tǒng)計功效的降低才可以近似忽略[1]。因此,多建議對缺失值進行填補后再進行分析。填補又稱插補,可分為單一填補法和多重填補法。單一填補包括均值填補,啞變量填補和基于回歸模型的單一填補等。單一填補因為會降低被填補變量的不確定性,縮小參數(shù)估計的標準差,已逐漸被淘汰,本文不作詳細介紹[2]。多重填補法因其考慮了缺失的不確定性等優(yōu)點越來越受到大家的推崇[3]。另外,廣義線性模型(generalized linear models, GLM),囊括線性回歸模型、二元Logistic回歸模型、Poisson回歸模型等,在臨床數(shù)據(jù)分析中占據(jù)著重要地位,其中二元Logistic回歸常常作為主要的多因素模型,用以產(chǎn)出因果推斷的統(tǒng)計學依據(jù)[4-5]。目前,雖然二元Logistic回歸與多重填補技術(shù)的結(jié)合在危險因素分析的醫(yī)學研究中的應用越來越多,如Foerster等[6]通過二元logistic回歸和多重填補方法建立了一個風險分層模型以更好的識別適合上尿路上皮癌內(nèi)鏡下保留腎臟手術(shù)的患者,但以實例分析介紹兩者如何結(jié)合應用的文獻較少。本文就數(shù)據(jù)缺失的模式和比例、多重填補的流程進行簡單的梳理,并使用急性腎損傷的數(shù)據(jù)演示多重填補(mice程序包)與廣義線性模型的結(jié)合應用,以期為含缺失值的臨床數(shù)據(jù)的統(tǒng)計分析提供參考依據(jù)。

數(shù)據(jù)缺失的模式與比例

多重填補的創(chuàng)始人Little和Rubin教授將數(shù)據(jù)缺失的原因分為以下三種模式:完全隨機缺失(missing completely at random, MCAR)、隨機缺失(missing at random, MAR)和非隨機缺失(missing not at random, MNAR)。完全隨機缺失指的是數(shù)據(jù)的缺失與否既不能歸因于已觀測的變量,也不能歸因于未觀測的變量,即存在缺失的觀測與完整觀測來自于同一個分布的總體,完整案例分析方法僅適用于這種缺失模式[7]。但有學者認為這種缺失模式只是理論上存在的,是隨機缺失的一種特殊情況[8]。隨機缺失指的是數(shù)據(jù)的缺失與否取決于已觀測的變量,即存在缺失的觀測和完整觀測的差別可以被其他觀測變量解釋。有研究者認為,“隨機缺失”一詞的命名有歧義,其 “隨機”一詞與“完全隨機缺失”中的“隨機”意義是不同的,前者是一種有條件的、可控制的缺失,更精確的表達方式為隨機條件缺失(conditionally missing at random)[7-9]。非隨機缺失指的是數(shù)據(jù)的缺失與否取決于未觀測到的變量,即存在缺失的觀測和完整觀測的差別無法用所有已知變量解釋。將多重填補技術(shù)應用于含有前兩種缺失(即完全隨機缺失和隨機缺失)的數(shù)據(jù),可得到近似無偏的參數(shù)估計結(jié)果,而非隨機缺失則不在多重填補技術(shù)的應用范圍。在實際應用中,事先幾乎不可能知道數(shù)據(jù)的缺失究竟屬于哪種類型,那么在應用多重填補法對缺失數(shù)據(jù)進行處理之時,需要先假定數(shù)據(jù)缺失符合隨機缺失。只要構(gòu)建填補模型時納入足夠多的變量,數(shù)據(jù)缺失的模式就會非常接近隨機缺失,可以使用多重填補法進行處理[10-11]?；蛘呖梢赃x擇多個可能適用的多重填補模型分別填補,對各個模型所得完整數(shù)據(jù)集進行相同的統(tǒng)計分析,將各自的結(jié)果匯總,進行敏感性分析[12]。

在應用多重填補法處理缺失數(shù)據(jù)時,數(shù)據(jù)缺失的比例對多重填補的影響目前尚存在爭論[13]。人們所推薦的多重填補最多可接受的數(shù)據(jù)缺失的比例在5%～50%[14-16]。近年有學者對模擬數(shù)據(jù)進行了大量的填補運算,發(fā)現(xiàn)相比于數(shù)據(jù)缺失的比例,納入填補模型的輔助變量所起的作用更加明顯。即便是在缺失數(shù)據(jù)比例高達90%的情況下,只要在填補模型中納入合適的輔助變量,填補后的數(shù)據(jù)依然可以進行無偏倚的統(tǒng)計推斷[13]。然而此研究所填補變量均為符合正態(tài)分布的連續(xù)變量,且為模擬數(shù)據(jù),其結(jié)論可能不適用于分類變量和其他分布類型的變量。

因此,在進行多重填補之前,應對所需填補數(shù)據(jù)的缺失值進行多方面的評估,必要時需了解各個變量收集、整理、清洗的全過程,盡可能納入更多輔助變量進入填補模型,并使用多個填補模型,進行敏感性分析,方能盡可能減少偏倚。

多重填補的流程

以R語言程序包mice(multivariate imputation by chained equation)為例,由該程序包及相應函數(shù)實現(xiàn)數(shù)據(jù)多重填補的原理和過程可知 (圖1),一個完整的數(shù)據(jù)多重填補過程應包括缺失值的填補和使用所得數(shù)據(jù)進行統(tǒng)計分析兩大部分[12,17]。因此,當我們使用多重填補法對含缺失值的數(shù)據(jù)進行回歸分析時,建議按照如下流程進行數(shù)據(jù)填補和回歸分析:

圖1 基于多重填補的廣義線性模型的實現(xiàn)流程

查看數(shù)據(jù)缺失情況當我們拿到一個數(shù)據(jù)集后,首先應當明確數(shù)據(jù)的缺失模式,即前文中提到的MCAR、MAR和 MNAR,但這一點往往難以做到。不過,我們能做的是查看數(shù)據(jù)集的缺失情況,包括明確含缺失值的變量、變量類型和缺失比例,這對是否進行填補和多重填補中一些命令行參數(shù)的設定很重要。比如,若因變量有缺失,一般我們會選擇刪除相應的觀測行,而不對因變量進行填補;若缺失的變量既有連續(xù)變量,也有分類變量,那我們可以通過函數(shù)factor()先對分類變量進行定義,或在多重填補時為每個需要填補的變量設定正確的參數(shù)defaultMethod;若某個或某些變量的缺失比例較大(>50%),我們一般不建議進行數(shù)據(jù)填補,而建議放棄使用該變量。

缺失值的多重填補當缺失比例較低時,我們便可以采用函數(shù)mice()嘗試缺失值的多重填補。在這個過程中,有3個重要的參數(shù)需要設定。第一個為參數(shù)m,用來設定多重填補生成的完整數(shù)據(jù)集的數(shù)量,默認為5。第二個為method或defaultMethod,用來指定變量的填補方法,比如pmm對應定量變量、logreg對應二分類變量、polyreg對應無序多分類變量、polr對應有序多分類變量。第三個為seed,設置種子數(shù),使得下次再運行代碼時所得到的數(shù)據(jù)集與前次一致,保證可重復性。

完整數(shù)據(jù)集的統(tǒng)計分析在完成數(shù)據(jù)的填補之后,我們會得到一個mids對象,該對象包含了通過鏈式多元填補得到的多個完整數(shù)據(jù)集。由于函數(shù)glm()中的參數(shù)data所識別的數(shù)據(jù)集類型為data frame (數(shù)據(jù)框)、list(列表)或環(huán)境 (environment),那么我們就無法直接使用函數(shù)glm()對mids對象進行回歸分析。解決的方法是,聯(lián)用函數(shù)with()和glm(),對mids對象中的每個數(shù)據(jù)集均進行回歸分析。

統(tǒng)計分析結(jié)果的匯總在使用with()完成了每個數(shù)據(jù)集的回歸分析后,我們會得到每個數(shù)據(jù)集對應的回歸分析結(jié)果,即多套回歸分析結(jié)果。接下來,我們需要使用函數(shù)pool()對多個回歸分析結(jié)果進行匯總,綜合為包含回歸系數(shù)、標準誤、P值等在內(nèi)的一套回歸分析結(jié)果。

在廣義線性模型中的應用

案例情況如下:某腎內(nèi)科醫(yī)師回顧性收集了109例在院內(nèi)發(fā)生急性腎損傷和231例未發(fā)生急性腎損傷患者的相關(guān)信息,包括一般特征(性別、年齡、BMI等級)和入院時的實驗室檢查結(jié)果(中性粒細胞和淋巴細胞比值、尿素氮、尿肌酐、尿酸、鈉離子、鉀離子、超敏反應蛋白、血清乳酸),擬通過二元logistic回歸分析急性腎損傷的影響因素。

在廣義線性模型中,將glm()函數(shù)的famliy參數(shù)設置為binomial時,所擬合的模型便是二元logistic回歸模型。接下來,我們以此類回歸模型為例,展示基于多重填補數(shù)據(jù)集、缺失數(shù)據(jù)集和原始數(shù)據(jù)集的廣義線性模型的結(jié)果,并將基于多重填補的回歸分析結(jié)果與基于無缺失值原始數(shù)據(jù)的回歸分析結(jié)果進行比較(本文所用的R代碼均可通過github網(wǎng)站查看https://github.com/hamody316/mice/blob/main/README.md#mice)。

由于本次所使用的的數(shù)據(jù)集并無缺失值,因此我們使用R程序包simFrame中的函數(shù)setNA()為性別、年齡、體質(zhì)量指數(shù)(BMI)等級、尿肌酐、鈉離子和血清乳酸構(gòu)造了缺失值,使得他們均有5%的隨機缺失。另外,在該數(shù)據(jù)集中,性別為二分類變量,BMI等級為有序多分類變量(等級變量),其他均為連續(xù)變量。

假定我們已確定了納入到多因素回歸分中的變量為性別、年齡、BMI等級、尿肌酐、鈉離子和血清乳酸。在設定了填補數(shù)據(jù)集數(shù)量(m=10)、變量的填補方法和種子數(shù)后,我們得到了10個完整的數(shù)據(jù)集,使用這10個數(shù)據(jù)集進行回歸分析的匯總結(jié)果如表1。

表1 基于多重填補的回歸分析結(jié)果

為了能夠說明多重填補的效果,我們用含缺失值的數(shù)據(jù)集和無缺失值的原始數(shù)據(jù)集分別進行了回歸分析,所得結(jié)果見表2和表3;同時,我們還對表3和表1的回歸系數(shù)進行了差值計算(表4)。

表2 基于含缺失值數(shù)據(jù)集的回歸分析結(jié)果

表3 基于原始數(shù)據(jù)的回歸分析結(jié)果

表4 回歸系數(shù)的差值

由表3的結(jié)果可知,所納入的多個變量在基于多重填補的回歸分析中具有統(tǒng)計學意義,與急性腎損傷有統(tǒng)計學關(guān)聯(lián),這與表1中的結(jié)論一致。由表2的結(jié)果可知,在基于含缺失值的回歸分析結(jié)果中,年齡和BMI等級與急性腎損傷無統(tǒng)計學關(guān)聯(lián),這與基于原始數(shù)據(jù)的回歸分析的結(jié)果不一致,說明表2中的結(jié)果發(fā)生了偏倚。另外,由表4可知,基于多重填補數(shù)據(jù)得到的變量系數(shù)與基于原始數(shù)據(jù)得到的變量系數(shù)的差值很小,進一步說明多重填補得到結(jié)果的穩(wěn)健性和可靠性。當然,本文多重填補的結(jié)果并不意味著多重填補對所有含缺失值的數(shù)據(jù)集而言都是最佳的填補方法,我們應該根據(jù)數(shù)據(jù)類型、缺失情況等選擇合適的填補方法,比如隨機森林填補法、K鄰近值法等[18-19]。

小結(jié):本文對缺失值的常見處理方法、缺失值模式、缺失值比例、多重填補的流程進行了簡單總結(jié),并通過急性腎損傷的影響因素分析的實例,展示了基于多重填補的廣義線性模型的分析過程及結(jié)果。對比基于多重填補的回歸分析結(jié)果、基于原始數(shù)據(jù)的回歸分析結(jié)果和基于缺失數(shù)據(jù)的回歸分析結(jié)果可知,在缺失比例較低時,前兩者的回歸分析結(jié)果(如各變量的回歸系數(shù))雖然有一定的差值,但整體結(jié)論并未出現(xiàn)偏倚,即在回歸方程中各變量的統(tǒng)計學意義是一致的。這也表明,在條件適當?shù)那闆r下,可以在腎病相關(guān)數(shù)據(jù)的統(tǒng)計分析中使用多重填補。本文結(jié)合腎病相關(guān)案例數(shù)據(jù)展示了如何使用R語言實現(xiàn)該方法,并對該方法所得結(jié)果的穩(wěn)健性進行了驗證,希望能為廣大醫(yī)護或科研人員在處理缺失數(shù)據(jù)的思路和實踐上提供參考。