Hirosawa T, Harada Y, Morinaga K, et al. Eosinopenia as a diagnostic marker of bloodstream infection in a general internal medicine setting: a cohort study[J]. BMC Infectious Diseases, 2020, 20(1): 1-7.
BMC Infectious Diseases: 2018 Journal Impact Factor:?? 2.565
文章地址：https://bmcinfectdis.biomedcentral.com/articles/10.1186/s12879-020-4814-5
嗜酸性粒細(xì)胞減少的情況：
1. Cusing綜合征或糖皮質(zhì)醇激素使用
2. 應(yīng)激導(dǎo)致內(nèi)源性激素產(chǎn)生增加從
3. 炎性介質(zhì)釋放導(dǎo)致（急性感染）
4. 藥物導(dǎo)致Krause JR, Boggs DR (1987). "Search for eosinopenia in hospitalized patients with normal blood leukocyte concentration".?Am. J. Hematol.?24?(1): 55–63.?doi:10.1002/ajh.2830240108.?PMID?3799595.
Gil H, Magy N, Mauny F, Dupond JL (2003). "[Value of eosinopenia in inflammatory disorders: an "old" marker revisited]".?Rev Méd Interne?(in French).?24?(7): 431–5.?doi:10.1016/S0248-8663(03)00138-3.?PMID?12829215.
Bass DA, Gonwa TA, Szejda P, Cousart MS, DeChatelet LR, McCall CE:?Eosinopenia of acute infection: production of eosinopenia by chemotactic factors of acute inflammation.?J Clin Invest?1980,?65:?1265-1271.?10.1172/JCI109789
摘要
背景：嗜酸性粒細(xì)胞在血流感染的應(yīng)用價(jià)值很少研究。本研究目的為評估嗜酸性細(xì)胞減少是否能作為預(yù)測血流感染的標(biāo)志物。
方法：本研究為門診部和普通住院部的回顧性隊(duì)列研究。納入189個至少有2套血培養(yǎng)結(jié)果的病人，排除血培養(yǎng)前2周使用過抗生素、激素治療、血液系統(tǒng)腫瘤的病人。單因素的診斷價(jià)值使用AUROCs來評價(jià)。主要終點(diǎn)為血培養(yǎng)陽性。
結(jié)果：重度嗜酸性粒細(xì)胞減少（<24.4cells/mm3）的診斷價(jià)值較小但有統(tǒng)計(jì)學(xué)意義（AUROC：0.648,95%CI 0.547-0.748），靈敏度和特異度也不高。模型的變量包括基線變量（年齡、性別），CRP，中性粒細(xì)胞計(jì)數(shù)，AUROC為0.729，加入嗜酸性粒細(xì)胞減少后AUROC有少許增加，為0.758，同時(shí)net reclassification improvement (NRI)有顯著增加(P?=?0.003). 然而integrated discrimination index (IDI) 無明顯意義(P?=?0.284)
結(jié)論：嚴(yán)重嗜酸性粒細(xì)胞減少可以被考慮為一個廉價(jià)的血流感染標(biāo)志物，但其診斷可靠性并不高。
介紹
血培養(yǎng)診斷血流感染非常重要，但是其陽性率并不高（除了特定條件下如感染性心內(nèi)膜炎、腦膜炎），同時(shí)也有假陽性的問題。
迄今為止，尚無研究發(fā)現(xiàn)與感染嚴(yán)重程度和預(yù)后相關(guān)的高度敏感，特異性強(qiáng)，易于測量，快速且廉價(jià)的血流感染標(biāo)志物。雖然寒戰(zhàn)、CRP、qSOFA被認(rèn)為是潛在標(biāo)志物，但其靈敏度和特異度均不足。
嗜酸性粒細(xì)胞減少，被之前的研究鑒定為好的感染診斷標(biāo)志。盡管一些研究報(bào)告說，外周血嗜酸性粒細(xì)胞的缺乏不能用作醫(yī)院住院環(huán)境中菌血癥的臨床可靠標(biāo)志物。但這些研究患者數(shù)量較少，并不局限于普通內(nèi)科。我們假設(shè)嗜酸性粒細(xì)胞減少時(shí)血流感染的可靠指標(biāo)
方法學(xué)
研究設(shè)計(jì)及病人選擇
回顧性、單中心隊(duì)列包括門診和一般住院病人（排除急診和ICU），血培養(yǎng)時(shí)間:1 January to 31 December, 2018.This study was conducted in accordance with the current version of the Declaration of Helsinki. The study protocol was approved by the institutional ethics committee of Dokkyo Medical University (No. R-20-18?J).
病人
399 adult patients (age?>?15?years)在研究時(shí)間段內(nèi)有血培養(yǎng)結(jié)果，排除178個培養(yǎng)前2周內(nèi)使用抗生素的，25個使用激素的，2個血液系統(tǒng)腫瘤的，還有5個病人因?yàn)槿狈Y料而排除。最終納入189個病人。
結(jié)局指標(biāo)和定義
血培養(yǎng)提示的血流感染是主要結(jié)局指標(biāo)，血流感染的定義為至少一次血培養(yǎng)陽性。單次常見背景污染菌陽性定義為陰性，常見背景菌包括：coagulase-negative Staphylococcus species（凝固酶陰性的葡萄球菌）, Bacillus species（芽孢桿菌屬）, Propionibacterium acnes（痤瘡丙酸桿菌）or Corynebacterium species（棒狀桿菌屬）不同研究對嗜酸性粒細(xì)胞下降的定義不同，在該研究中，嗜酸性粒細(xì)胞減少通過單因素分析定義<?24.3 cells/mm3，qSOFA的定義respiratory rate of ≥22/min, altered mentation（GCS<15）, or systolic blood pressure of ≤100?mmHg 血培養(yǎng)步驟 10ml血分為需氧和厭氧培養(yǎng)為1套，在另一側(cè)的靜脈再采10ml分為無氧和需氧培養(yǎng)為第2套。接種進(jìn)行菌種鑒定（were confirmed using matrix-assisted laser desorption ionisation–time of flight mass spectrometry.）和藥敏。
數(shù)據(jù)收集
2個主治醫(yī)生確認(rèn)培養(yǎng)出細(xì)菌的致病性，體查數(shù)據(jù)在采集血液時(shí)記錄，其他marker與血液采集的時(shí)間在一天內(nèi)測量。
分析
Continuous variables are presented as medians and interquartile ranges [25th–75th percentiles] and were compared using the Mann–Whitney U test. Categorical or binary variables are presented as numbers (percentages) and were compared using the chi-squared test or Fisher’s exact test.?連續(xù)變量使用中位數(shù)和四分位間距，并用MW U檢驗(yàn)比較，分類變量用百分?jǐn)?shù)描述和卡方或fisher比較，單因素和多因素logistic回歸的準(zhǔn)確性用AUROC評價(jià)。
之前的研究鑒定出CRP是血流感染的有效預(yù)測指標(biāo)，計(jì)算integrated discrimination index (IDI) and net reclassification improvement (NRI)評估加入嗜酸性粒細(xì)胞后是否提高診斷價(jià)值。內(nèi)部驗(yàn)證：ordinary nonparametric bootstrapping with 1000 bootstrap samples and bias-corrected, accelerated 95% CIs
結(jié)果
189個病人納入研究，在研究期間有25個血培養(yǎng)陽性，80個嗜酸性粒細(xì)胞降低。4個病人培養(yǎng)出多種菌，12個培養(yǎng)出G+菌，9個培養(yǎng)出G-菌。培養(yǎng)陰性和陽性的病人基線資料見表格，其中WBC和CRP有顯著差異。
其他可能影響嗜酸性粒細(xì)胞的資料見表S1
結(jié)果
表2位單因素分析結(jié)果，嗜酸性粒細(xì)胞計(jì)數(shù)（both before and after bootstrapped multiple regression analysis (1000 bootstrap replicates)）、中性粒細(xì)胞計(jì)數(shù)、CRP均有統(tǒng)計(jì)學(xué)意義，而白細(xì)胞計(jì)數(shù)、qSOFA和寒戰(zhàn)均無統(tǒng)計(jì)學(xué)意義。Cut-off后，CRP靈敏度高特異度低，中性粒細(xì)胞計(jì)數(shù)靈敏度低特異度高
表3展示了預(yù)測血流感染模型的AUROCs，基線模型加入中性粒細(xì)胞計(jì)數(shù)和CRP后提高了AUROC（0.650到0.729），IDI?和NRI均有統(tǒng)計(jì)學(xué)意義；在次基礎(chǔ)上加上嗜酸性粒細(xì)胞計(jì)數(shù)后，AUROC輕度增加到0.758，NRI有統(tǒng)計(jì)學(xué)意義，IDI無統(tǒng)計(jì)學(xué)意義。Bootstrap多次后最終模型的AUROC仍相同。圖1展示了三個模型的ROC曲線。
討論
總結(jié)發(fā)現(xiàn)：According to our findings, 預(yù)測能力有統(tǒng)計(jì)學(xué)意義但l靈敏度和特異度不高，比qSOFA及寒戰(zhàn)診斷價(jià)值高（該隊(duì)列不包括ICU和急診的病人），診斷模型有輕度改良，可能可以作為預(yù)測標(biāo)志物，然而還需要更多研究。
本研究的穩(wěn)健性，與其他研究的不同（可能的原因?yàn)閝SOFA和寒戰(zhàn)僅可用于ICU及急診病人）
研究的局限性：
1.單中心、單科室不能推廣到ICU和急診。
2.排除了178個2周內(nèi)使用過抗生素的病人（如此多病人可能以為本研究在三甲醫(yī)院，很多病人在就診前就已經(jīng)去過其他醫(yī)院看?。?，可能有選擇偏倚
3.對于血培養(yǎng)沒有嚴(yán)格的指征，由醫(yī)生決定
4.排除了血液系統(tǒng)腫瘤、嗜酸性粒細(xì)胞增多、激素使用的病人，因?yàn)檫@些病人極少見。因此該結(jié)論也不適用與這些疾病流行度高的病人。
5.僅輕度和中度寒戰(zhàn)的病人納入了研究，之前研究表明寒戰(zhàn)程度越高尤其是伴有發(fā)抖時(shí)，菌血癥可能性越大，這可能是本研究中寒戰(zhàn)沒有統(tǒng)計(jì)學(xué)差異的原因。
6.本研究沒有納入PCT（很多病人沒有做該檢測），盡管其預(yù)測培養(yǎng)陽性的能力不足，但降鈣素原被認(rèn)為是血流感染的預(yù)測標(biāo)志物。
審稿意見與答復(fù)：
Statistical method description is also not clear. They did no attempt to test validity of the prediction model using some methods such as?cross-validation and bootstrapping. Correcting optimism problem in study with small sample size have been more emphasized. The conclusions drawn inadequately supported by the data shown.?
Response:?As reviewer 1 pointed out,?we tested validity of the prediction model using bootstrapping and cross-validation.?We added the result of bootstrapping in the manuscript. Eosinopenia alone yielded almost same overall predictive ability (AUROC: 0.648, 95%CI: 0.491–0.805, cut-off = 24.4 cells/mm3) from 10-fold cross validation procedure. The addition of eosinopenia to the baseline variables (age, sex) and CRP similarly improved the AUROC (from 0.650 to 0.753, 95%CI 0.662-0.851, cut-off = 24.4 cells/mm3) from 10-fold cross validation procedure. We also changed the conclusions supported by the data shown.
統(tǒng)計(jì)學(xué)方法疑惑：
logistic回歸cross-validation and bootstrapping如何實(shí)現(xiàn)（from 10-fold cross validation procedure.?Following bootstrapped multiple regression analysis (1000 bootstrap replicates)）Internal validation of the prediction models was conducted using ordinary nonparametric bootstrapping with 1000 bootstrap samples and bias-corrected, accelerated 95% CIs.[Ref:Efron B TR. An introduction to the bootstrap: CRC press; 1994.]
基本理論
即把所有樣本分為10份（彼此不交叉），其中9份作為訓(xùn)練集，1份作為測試集（共10組模型）（https://www.zhihu.com/question/39259296）
上圖中訓(xùn)練生成10組模型，10組模型又分別對自己的測試集產(chǎn)出測試結(jié)果，測試結(jié)果的平均值就是返回結(jié)果。10組模型的預(yù)測結(jié)果差距越小，表明模型參數(shù)魯棒性越好，出現(xiàn)過擬合的可能性也越小，因?yàn)榭梢宰C明模型在不同數(shù)據(jù)集上表現(xiàn)穩(wěn)定。cross-validation的最大好處是所有數(shù)據(jù)都會參與到訓(xùn)練和預(yù)測中。
cross-validation有一個問題，就是在數(shù)據(jù)分組的時(shí)候缺乏隨機(jī)性，而且每個數(shù)據(jù)只能固定屬于K組中的一組，可能會造成對于最終結(jié)果的影響。Booststrap的思想是完全隨機(jī)放回性采樣。Booststrap的用法如下：
1.數(shù)據(jù)分組
假設(shè)要分成10組，則先設(shè)置一個采樣比例，比如設(shè)置成70%。則10組數(shù)據(jù)是每次從數(shù)據(jù)集D中隨機(jī)采樣總數(shù)70%的數(shù)據(jù)構(gòu)成數(shù)據(jù)集D1，然后把數(shù)據(jù)放回，再隨機(jī)采樣總數(shù)70%的數(shù)據(jù)構(gòu)成數(shù)據(jù)集D2........以此類推，放回式采樣10份。
2.訓(xùn)練生成10組模型
3.效果評估（取平均值）
如果用來做模型的選擇和預(yù)測誤差的估計(jì)，cross-validation的特點(diǎn)是bias小，方差大，方差基本不會因?yàn)閒old數(shù)目提高有明顯減小，Booststrap的特點(diǎn)是bias大但是方差小。bias描述的是根據(jù)樣本擬合出的模型的輸出預(yù)測結(jié)果的期望與樣本真實(shí)結(jié)果的差距，簡單講，就是在樣本上擬合的好不好。
評價(jià)診斷指標(biāo)的參數(shù)NRI和IDI
關(guān)于net reclassification improvement (NRI)【新模型對舊模型的重新分類導(dǎo)致準(zhǔn)確度的改善，重分類改善指標(biāo)】和integrated discrimination index (IDI)【綜合判別改善指數(shù)】的詳細(xì)解讀見：
https://cjasn.asnjournals.org/content/7/8/1355
NRI：The net proportion of patients with events reassigned to a higher risk category (NRIevents) and of patients without events reassigned to a lower risk category (NRInonevents) is calculated (Appendix). The NRI is the sum of NRIevents and NRInonevents. It is interpreted as the proportion of patients reclassi?ed to a more appropriate risk category.
The IDI is the sum of IDIevents and IDInonevents. The IDIevents is also equal to the difference in the average sensitivity (normally termed the integrated sensitivity [IS]) of the two models across all risk thresholds and the IDI for those without events (IDInonevents) is equal to the difference in the average 1-speci?city (integrated 1-speci?city [IP]).?
以下內(nèi)容來源 "解螺旋·臨床醫(yī)生科研成長平臺" 原創(chuàng)文章：怎樣判斷一個診斷（風(fēng)險(xiǎn)預(yù)測）模型的好壞？
模型的區(qū)分度和校準(zhǔn)度
一個好的預(yù)測模型應(yīng)該能針對某個結(jié)局，把一群人的風(fēng)險(xiǎn)高低區(qū)分開來，這就是區(qū)分度（discrimination）。它跟患者在人群中的分布特征有關(guān)，模型中納入的自變量（如性別、年齡、某些實(shí)驗(yàn)室檢查等）如果有異質(zhì)性，那么模型的區(qū)分度就好，否則就差。區(qū)分度一般以我們熟悉的ROC曲線下面積（AUC）來評價(jià)，或稱C統(tǒng)計(jì)量（C-statistics）。AUC越高，模型對高低風(fēng)險(xiǎn)人群的區(qū)分度越好。一個廣泛接受的評判標(biāo)準(zhǔn)是，AUC在0.6以下為低區(qū)分度，0.6 – 0.75 是中區(qū)分度，0.75以上為高區(qū)分度。
但一個區(qū)分度很好的模型，卻可能有著較差的校準(zhǔn)度（calibration），或稱擬合優(yōu)度（Goodness?of fit）。比如它能判斷一個人發(fā)生某疾病的風(fēng)險(xiǎn)是另一個人的5倍，它判斷這兩人的風(fēng)險(xiǎn)分別為5%和1%，而實(shí)際上兩人的風(fēng)險(xiǎn)是50%和10%，那這個模型也挺離譜的，這就是校準(zhǔn)度不好。模型的校準(zhǔn)度可以用Hosmer-Lemeshow檢驗(yàn)（H-L檢驗(yàn)），若結(jié)果得到顯著統(tǒng)計(jì)學(xué)意義，則表明預(yù)測值與觀測值之間（即模型和真實(shí)之間）有差異。但H-L檢驗(yàn)有其不足。有統(tǒng)計(jì)學(xué)差異只能提示這個模型整體上跟觀測值有差異，卻無法展示更多細(xì)節(jié)。做出校準(zhǔn)曲線來有可能會看到這樣一種情況：?
這是MAGGIC心衰風(fēng)險(xiǎn)評分量表的校準(zhǔn)曲線，其用于評價(jià)心衰1年死亡風(fēng)險(xiǎn)。黃色曲線為模型預(yù)測值，藍(lán)色帶有誤差線的點(diǎn)為觀測值?？梢娝鼈兦懊嬉欢蔚惋L(fēng)險(xiǎn)時(shí)吻合得不錯，而有輕微高估；但30%開始，模型預(yù)測的風(fēng)險(xiǎn)低于實(shí)際觀察值10%左右，自此拉開差距。這個模型在外部驗(yàn)證集中檢測的時(shí)候，AUC為0.77。像這種情況的模型是否毫無用處呢？如果有其他研究表明（或行業(yè)內(nèi)達(dá)成共識），風(fēng)險(xiǎn)高于20%的時(shí)候就一定要采取某種干預(yù)，那么后面的差異其實(shí)影響不大，這都是可以在討論中Argue一下的，深入分析自己的模型的價(jià)值和適用范圍，也可圈可點(diǎn)。區(qū)分度和校準(zhǔn)度對一個模型來說都是很重要的評價(jià)，許多新開發(fā)的模型沒有得到充分的評價(jià)，就只好流散于蒼茫文獻(xiàn)海，泯然一滴水，后人遇到了還要辛苦甄別。有一項(xiàng)關(guān)于心血管系統(tǒng)風(fēng)險(xiǎn)預(yù)測模型研究的系統(tǒng)綜述發(fā)現(xiàn)，只有63%的模型匯報(bào)了區(qū)分度，而匯報(bào)校準(zhǔn)度的更少，才36%。
你的模型要和別人的模型來一場PK，如果你開發(fā)的風(fēng)險(xiǎn)預(yù)測模型是為了解決一個全新的問題，展示自己好也就夠了。但更多時(shí)候是希望能改進(jìn)現(xiàn)有的解決方案，那么當(dāng)然還要表明它比現(xiàn)有的好（至少某一方面），才能有發(fā)表的機(jī)會呀。區(qū)分度比較兩者的AUC是可以的，但沒什么細(xì)節(jié)。校準(zhǔn)度的比較則可用赤池信息準(zhǔn)則（Akaike Index Criterion, AIC）或貝葉斯信息準(zhǔn)則（Bayesian Index Criterion, BIC），其值越低，校準(zhǔn)度越好，但同樣沒什么細(xì)節(jié)，可能相對較好的那個模型離真實(shí)情況還是挺遠(yuǎn)的。但AIC和BIC見得更少了。這里再重點(diǎn)介紹一種近年來應(yīng)用愈發(fā)增多的指標(biāo)，聽說有的小伙伴已經(jīng)遇到過審稿人要求返修做這個分析了，就是凈重新分類指數(shù)（Net?Reclassification Index, NRI）。聽名字，這意思就是舊模型把一群人分為高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)，新模型會把其中幾個人挪挪窩，有的挪對了有的挪錯了，“挪對的 – 挪錯的”就是凈重分類。指數(shù)嘛，再除以總?cè)藬?shù)的百分比咯。但還要分成兩撥來看，即觀測發(fā)生結(jié)局事件的和未發(fā)生結(jié)局事件的，因?yàn)樗鼈兎謩e代表假陰性和假陽性的概率。還是在一個情景中說明吧。
相加NRI（additive NRI）的值可以從200 ~ -200，即新模型完全分類正確（所有低風(fēng)險(xiǎn)的人都未發(fā)生結(jié)局事件，且所有高風(fēng)險(xiǎn)的人都發(fā)生了結(jié)局事件）為200，反之完全分類錯誤為-200。它主要的局限性就是沒有考慮發(fā)生與不發(fā)生事件的兩撥人在總體中的分布情況。
絕對NRI（absolute NRI）以總體為分母，則回避了這個局限。
相加NRI用得多一點(diǎn)，但其實(shí)它們各有意義，應(yīng)該兼看。再看兩個例子來理解它倆的差別吧。
例1：?樣本20000人，觀察下來，發(fā)生與未發(fā)生結(jié)局事件的兩撥人各有10000。下圖綠色表示正確重新分類，粉色表示錯誤重新分類，灰色表示沒有重新分類。
計(jì)算下來，相加NRI為12，是正值；絕對NRI是6%，也是正值，如圖下方的表格所示。也就是說，雖然未發(fā)生結(jié)局事件的有300人被錯誤重新分類，但總體上來說還是利大于弊的。
?例2：樣本是11000人，觀察下來，發(fā)生與未發(fā)生結(jié)局事件的人比例不再是1:1了，分別為1000人與10000人。
如表所示，這回相加NRI和上邊一樣是12，還是正值。如果光看相加NRI，似乎新模型還是個好模型。但看絕對NRI，原來是個負(fù)值。
再看看細(xì)節(jié)，新模型對發(fā)生了結(jié)局事件的人做了更好的重分類，可遺憾的是，這撥人所占比例更?。欢?b>對于未發(fā)生結(jié)局事件的人，其重分類效果并不好，這部分人占的比例又比較大，最終絕對NRI是負(fù)的，暴露了新模型的局限性。
下面是一個研究中的真實(shí)案例，情況跟上邊的例2差不多。評價(jià)術(shù)后發(fā)生心?；蛐脑葱运劳龅娘L(fēng)險(xiǎn)，原來有一個評分量表（Revised Cardiac Risk Index, RCRI），此為舊模型。研究者想看看加上冠脈造影（CCTA）能不能更好地預(yù)測風(fēng)險(xiǎn)，即CCTA + RCRI為新模型。
觀察下來，結(jié)局發(fā)生與不發(fā)生在總體中的分布為7.7% : 92.3%，統(tǒng)計(jì)重新分類的情況，計(jì)算兩個NRI，出現(xiàn)了相加NRI為正值、絕對NRI為負(fù)值的情況。
發(fā)生這種情況，又要拋棄新模型咯？粗看或許的是，但其實(shí)還有另一個討論思路，這個思路在許多研究中都適用，就是假陰性與假陽性哪個后果更嚴(yán)重。
發(fā)生結(jié)局事件的人，當(dāng)初被錯誤地分到低風(fēng)險(xiǎn)組，就是假陰性，讓病人和醫(yī)生都放松警惕，可能錯過了及時(shí)干預(yù)的機(jī)會。未發(fā)生結(jié)局事件的人，當(dāng)初被錯誤地分到高風(fēng)險(xiǎn)組，就是假陽性，讓人提心吊膽草木皆兵，可能過度干預(yù)造成額外的負(fù)擔(dān)。
如果結(jié)局事件挺嚴(yán)重的，比如本例中的心?；蛐脑葱运劳?，就算風(fēng)險(xiǎn)低我也想要極力避免，要干預(yù)，那么假陽性較多的模型也有其價(jià)值。如果結(jié)局事件不那么嚴(yán)重，或者即使發(fā)生了還有較好的補(bǔ)救方案，那么當(dāng)初測試時(shí)它低估了我的風(fēng)險(xiǎn)我也不是太介意，假陰性多一點(diǎn)的模型也能接受?？傊绻愕姆治龇椒芸吹礁嗉?xì)節(jié)，也就給了你更多更深入分析討論的機(jī)會，為自己的研究成果找到亮點(diǎn)。?
Alba, A. C. et al. Discrimination and Calibration of Clinical Prediction Models: Users’ Guides to the Medical Literature. JAMA 318, 1377–1384 (2017).
關(guān)于假陰性與假陽性的權(quán)衡，還有其他指標(biāo)，比如凈受益率，可以做決策曲線來分析。
R語言計(jì)算NRI和IDI（來源：解螺旋文章）
R里有2個包專事計(jì)算NRI，分別為nricens和PredictABEL。從最后結(jié)果來說，nricens計(jì)算出來的是絕對NRI，PredictABEL則為相對NRI。但我們已經(jīng)知道計(jì)算原理了呀，而且它們都能生成新舊模型分類的對照表，所以其實(shí)只用其中一個包就都可以計(jì)算了。
先用一份示例數(shù)據(jù)做個模型。這是survival包里帶有的一份梅奧診所的數(shù)據(jù)，記錄了418位患者的臨床指標(biāo)，觀察這些因素與原發(fā)性膽汁性肝硬化（PBC）的關(guān)系。其中前312位參加了RCT，其他的只參加了觀察隊(duì)列。
library(survival)
### logistic回歸
egData <- pbc[1:312, ]
這個表中的結(jié)局是status，0 = 截尾（刪失），1 = 接受移植，2 = 死亡。研究目的“死亡與否”是個二分類變量，所以要做些變換。再看time一欄，有的不夠2000天，這些樣本要么是沒到2000天就死亡了，要么是刪失了。我們要刪掉2000天內(nèi)刪失的數(shù)據(jù)。
egData = egData[egData$time > 2000 | (egData$time < 2000 & egData$status == 2), ]
選好后做一個event向量，把status的三個狀態(tài)變成死亡 = 1, 未死亡 = 0。
#?當(dāng)time<2000，且status為2（死亡）時(shí)，記為1，否則為0。
event = ifelse(egData$time < 2000 & egData$status == 2, 1, 0)
# 合并數(shù)據(jù)
egData = cbind(egData,event)
然后選擇模型的自變量（predictors）。太多了，選取其中幾個示例。就以年齡、膽紅素、白蛋白為舊模型（standard），三者加上一個凝血酶原時(shí)間為新模型（new）。
一般做logistic回歸是用glm(因變量 ~ 自變量1 + 自變量2 + …… +自變量n,?family = binomial('logit'),?data = 數(shù)據(jù)表)，但如果自變量較多的話，前面那個運(yùn)算式就會很長很長，萬一這些自變量還是基因名或編號，就很想死了。所以順便講一個簡化的辦法，即把那一串先寫成formula。
# 舊模型
fml.std <- as.formula(paste0('event~', paste(colnames(egData)[c(5,11,13)], collapse = '+')))
#?event~age+bili+albumin
# 新模型
fml.new <- as.formula(paste0('event~', paste(colnames(egData)[c(5,11,13,19)], collapse = '+')))
#?event~age+bili+albumin+protime
# 模型擬合，x=TRUE是將來用nricens包計(jì)算時(shí)要求用到的，表示輸出結(jié)果中是否包含所用到的數(shù)據(jù)表，平時(shí)可以不寫。
mstd = glm(fml.std, family = binomial('logit'), data = egData, x=TRUE)
mnew = glm(fml.new, family = binomial('logit'), data = egData, x=TRUE)
模型就這樣做完了~ 先不急著計(jì)算NRI，先看看它的總體情況。
summary(mstd)
# 結(jié)果描述
Call: glm(formula = fml.std, family = binomial("logit"), data = egData, x = TRUE)
Deviance Residuals: Min 1Q Median 3Q Max
-2.78723 -0.57296 -0.39507 0.09268 2.82804
Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.01472 1.79654 -0.008 0.993461
age 0.05847 0.01600 3.655 0.000258 ***
bili 0.39126 0.06522 5.999 1.98e-09 ***
albumin -1.49447 0.42575 -3.510 0.000448 ***
?---
?Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 371.21 on 311 degrees of freedom
Residual deviance: 243.77 on 308 degrees of freedom
AIC: 251.77
Number of Fisher Scoring iterations: 6
殘差、相關(guān)系數(shù)、各個自變量的統(tǒng)計(jì)顯著性等，注意倒數(shù)第二行的AIC，就是上一期提到的赤池信息準(zhǔn)則，表示模型校準(zhǔn)度，很少有人匯報(bào)呢。
NRI的計(jì)算
# install.packages("nricens")
library(nricens)
NRI <- nribin(mdl.std = mstd, mdl.new = mnew, updown = 'category',cut = c(0.3,0.6), niter = 10000, alpha = 0.05)
填上新舊兩個模型。Cut是判斷風(fēng)險(xiǎn)高低的臨界值，現(xiàn)在我們寫了2個，也就是0~29%為低風(fēng)險(xiǎn)，30%~59%為中風(fēng)險(xiǎn)，60%~100%為高風(fēng)險(xiǎn)。現(xiàn)實(shí)中可以查閱相關(guān)文獻(xiàn)進(jìn)行設(shè)置，預(yù)測風(fēng)險(xiǎn)達(dá)到多少需要怎樣干預(yù)之類的。
Updown為定義一個樣本的風(fēng)險(xiǎn)是否變動的方式，category是指分類值，即我就熟悉的低、中、高風(fēng)險(xiǎn)，另有一種diff，為連續(xù)值。選diff時(shí)，cut就設(shè)1個值，比如0.02，即認(rèn)為當(dāng)預(yù)測的風(fēng)險(xiǎn)在新舊模型中相差2%時(shí)，即被認(rèn)為是重新分類了。這種方法用的比較少。
后面幾個參數(shù)就比較有意思了，niter為重復(fù)取樣的次數(shù)，即boostrap方法，不做的話將其設(shè)為0就好了；做的話建議至少1000次，這也是默認(rèn)值，但我（讀書少）見過的研究都10000次。然后將統(tǒng)計(jì)顯著性alpha設(shè)為0.05。
如果不做bootstrap，就是這個結(jié)果。有重新分類情況的詳表，最后是NRI和各種變動的概率。第一個NRI如前所述，是絕對NRI，大家可以根據(jù)之前的知識和上邊的詳表自己計(jì)算驗(yàn)證一下，此時(shí)可手動計(jì)算出相對NRI。其他指標(biāo)隨便看看。
如果做了bootstrap，就會多出一個表：
因?yàn)樽隽?0000次重復(fù)取樣，相當(dāng)于有10000個NRI，于是就有了標(biāo)準(zhǔn)誤和置信區(qū)間，剛才我們設(shè)alpha = 0.05，所以后面的Lower和Upper就是95%置信區(qū)間的下界和上界。
同時(shí)，做不做bootstrap都會得到一張圖，表示各數(shù)據(jù)點(diǎn)在新舊模型中的分布。
PredictABEL包的做法
# install.packages("PredictABEL")
library(PredictABEL)
pstd <- mstd$fitted.values
pnew <- mnew$fitted.values
先把兩個模型中的預(yù)測風(fēng)險(xiǎn)值提出來，也就是模型中的fitted.value。這個包只能從預(yù)測風(fēng)險(xiǎn)計(jì)算，剛才的nricens包可以用模型，也可以用預(yù)測風(fēng)險(xiǎn)（把mdl.std和mdl.new參數(shù)換成p.std和p.new）。
reclassification(data = egData, cOutcome = 21, predrisk1 = pstd,predrisk2 = pnew, cutoff = c(0,0.30,0.60,1))
cOutcome是結(jié)局事件的列序號，剛才我們不是把event放到最后了么，即第21列，填上。兩個預(yù)測風(fēng)險(xiǎn)值也相應(yīng)填上。這里的cutoff跟剛才的不一樣，還要填上前面的0和后面的1，成為完整的0~100%的區(qū)間。
然后得到一個重新分類表：
跟上邊nricens做的差不多了。不過這個包沒有bootstrap的選項(xiàng)。
接著看下面的結(jié)果，這里的分類NRI是咱們上回說的相加NRI，同樣可以根據(jù)上一期的知識手動計(jì)算一下。記得咱們并沒有設(shè)置bootstrap吧？可這里也有個95%置信區(qū)間，只是內(nèi)部調(diào)用了一個更為簡陋的只能計(jì)算連續(xù)NRI的improveProb()函數(shù)做的，而且連續(xù)變化的臨界值也不太透明，遂不管了。
最后還有個IDI是指，發(fā)生和未發(fā)生結(jié)局事件樣本的平均預(yù)測風(fēng)險(xiǎn)差異，在新模型中提高了0.44%。
這兩個包當(dāng)然各有優(yōu)劣。nricens計(jì)算時(shí)可控制的參數(shù)較多，匯報(bào)起來顯得華麗一些。PredictABEL計(jì)算結(jié)果則多了個IDI和大家喜聞樂見的p值。但也有學(xué)者表示，兩個模型的差異未必要求p<0.05。
PredictABEL包還有很多有用的功能，比如可以做Hosmer- Lemeshow校準(zhǔn)曲線，當(dāng)然也附送p值
plotCalibration(data = egData, cOutcome = 21, predRisk = pnew, groups=6, rangeaxis=c(0,1), plottitle="Calibration plot", xlabel="Predicted risk", ylabel="Observed risk", filename="calibration_table.txt")
兩個模型的ROC曲線
labels <- c("standard model", "new model")
plotROC(data = egData, cOutcome = 21, predrisk = cbind(pstd,pnew), labels=labels)
區(qū)分度箱形圖：
