晨晰統計部落格新站（統計、SPSS、BIG DATA討論園地）

筆者在之前的文章中（https://reurl.cc/OXXpR9），淺白地介紹新冠疫苗研究的「有效率」與「安全性數據」，本篇文章繼續介紹新聞稿（https://reurl.cc/gmO2vV）中的專有名詞釋義。

三、如何解讀次群體分析的結果

上單元我們已實例示範RCS在存活分析及邏輯斯迴歸分析的應用，本單位我們繼續介紹RCS在線性迴歸分析的應用

範例三、分析醫師的執刀經驗(累積執刀數目)情形與病患住院日數的相關性

傳統迴歸分析模型的一重要假設為自變項和依變項呈線性關係，然實際收案資料在分析時甚少能符合此假設條件，若能以模式模擬自變量與依變項之間的非線性關係，將可更忠實呈現自變項和依變項的相關。其中Restricted cubic spline (RCS) 可清楚描述自變項與依變項之間的關係（無論線性或非線性關係），且適用於各分析模式（包括線性迴歸、邏輯斯迴歸、存活COX迴歸分析等），而成為近年最常見的非線性相關分析方法。本單元我們依序實例示範RCS在各分析的應用(使用R statistic)

首先請先安裝及執行相關package

在前一篇文章中，我們示範了如何在二階層橫斷面資料中（人/時間點）點選 GLMM 的操作，本篇文章將仔細介紹報表的解讀。

下圖為報表最先出現的「Case Processing Summary」，。首先要注意資料有沒有遺漏值，有可能結果變項（Y）或解釋變項（X）具有遺漏值，此時會刪除在任何一個變項具有遺漏值的樣本（listwise deletion），本例中剛好有一筆遺漏值。下方「Model Summary」則顯示機率分佈為常態分佈以及連結函數為 identity，此即為標準的 linear regression。

筆者在不久之前示範了線性混合模式（Linear mixed model, LMM）以SPSS 25版應用於二個階層橫斷面資料的操作與解讀（請見：https://reurl.cc/9X6qrV、https://reurl.cc/m9R2LM、https://reurl.cc/ygD0EO），本系列文章則介紹二個階層縱貫面資料（longitudinal）的操作與應用。

我們使用一個真實的耳鼻喉科的資料作為示範，此資料一共包括 122 名阻塞睡眠呼吸中止症（Obstructive sleep apnea, OSA）的未成年患者，他們全部都接受扁桃腺線樣體切除手術治療，且術前及術後第三個月、術後第六個月皆有完整的 24 小時血壓測量數據，因此該 SPSS 資料一共 122*3 = 366 筆的橫列資料筆數。

新冠肺炎噬虐全球，已造成數以百萬計的民眾死亡，以及目前為止在許多國家仍有嚴重災情。而新冠肺炎疫苗的研發、試驗與上市，則成為了近期最為最大的事件，例如國際大藥廠輝瑞（Pfizer）德國生技公司 BioNTech 的 BNT162b2 疫苗試驗最終結果顯示「疫苗有效率」為 95％，而此結果也已經於 2020/12/10 發表在最頂尖醫學雜誌「The New England Journal of Medicine; NEJM」（新英格蘭醫學雜誌），全文 PDF 檔可以免費下載（https://reurl.cc/6l1ZVk）。

最近幾週（筆者撰寫文章的時間為 2020/12/20）也有許多關於疫苗研究的新聞，但畢竟不是每位民眾都具有醫學研究的背景知識，因此新聞稿上面會有各樣的醫學研究的專有名詞，本篇文章就略微解釋常見的幾個名詞，以幫助廣大民眾在閱讀相關新聞時，能有正確的觀念。

筆者前一陣子去參加「2020科技部生科司舉辦的心臟學門學術研討會：運用資料庫解決臨床研究問題的心法寶典」，會議中聽到了一些新的觀念以及未來的展望想跟各位讀者分享，首先真實世界的證據(RWE, Real World Evidence)相較於過去變得越來越重要，因為現在資料取得比起過去更加的便利與快速，科技的進步也使得這些大量的資料能被使用者分析與應用，以新藥的開發為例，新藥的開發到上市，都會透過進行嚴格的臨床試驗研究(Randomized controlled trial, RCT)來去評估藥物安全性及有效性，但因為執行一個RCT需要花費大量的金錢與人力，而且往往需要一段很長的時間，且RCT狹隘的收案條件往往無法使試驗結果外推至更廣的族群(缺乏真實世界的情況)，長期的療效及安全性也無法從有限的追蹤期間得知。因此，運用真實世界數據 (Real-world data; RWD)進行研究所產生的真實世界證據 (Real-world evidence; RWE) 可提供較為貼近真實世界的情況。

在醫學統計的領域中，針對信度的分析，若研究者有一個新的測量工具想與黃金指標(Gold standard)或是過去常用的測量工具進行比較時，我們一般希望新的方法與黃金指標(Gold standard) 或是過去常用的測量工具兩者之間具有相當程度的一致性(agreement)，這樣才能證明新的測量工具具有一定的準確性。當針對類別型資料(categorical data)，實務上最常使用的是 kappa 統計量，計算方式如下:

如表格所示，假設今天研究者想比較新工具與黃金指標工具在測量上的一致性，藉以用來作為新工具是否能用於實務上的標準，因此蒐集了兩種工具使用在同樣條件(協助測量的人、測量的時間等等非隨機因素)下，針對同一群人進行測量的結果。

筆者前一陣子有遇到一個分析的案子，這個案子很特別，收集資料的方式是以季為單位，每一個有乾癬診斷的病人，會去收集以乾癬診斷的那一天開始起算，往後追蹤每一季的資料，因此每個人的資料筆數為每三個月一筆，直到2013/12/31或是退保日期(死亡)為止。因此這個資料的筆數相當的驚人，假設某個病人從我們資料中最早的2001/1/1開始往後追蹤，中間都沒發生死亡，因此這個病人最多會有52筆資料(一年有四季*13年的追蹤)，實際上，我們這個資料最後的筆數將近400萬筆資料，這應該也稱得上是所謂的大數據(Big data)了吧!

為了幫助讀者理解資料的長相，如下圖所示，每一筆的資料其實是一季的收集，以第一個人為例，這個人在13年的追蹤過程中，共有8筆(8季的資料)。

在前一篇文章中，我們示範了如何在二階層橫斷面資料中（人/眼睛）點選 GLMM 的操作，本篇文章將仔細介紹報表的解讀。

下圖為報表最先出現的「Case Processing Summary」，。首先要注意資料有沒有遺漏值，有可能結果變項（Y）或解釋變項（X）具有遺漏值，此時會刪除在任何一個變項具有遺漏值的樣本（listwise deletion）。下方「Model Summary」則顯示機率分佈為二項式及連結函數為 logit，此即為標準的binary logistic regression。

接著我們使用一個真實的眼科資料作為示範，此資料集一共 1000 筆整，來自於 500 名接受兩種不同近視雷射手術的收案對象，每位患者都同時接受左右眼的手術。「ID」為患者編號，「Eyeball」為左右眼的編號（1 為左眼、2 為右眼），「Group」為雷射方法（1 為 LASEK，0 為 FSBK）。結果變項為「Outcome」，1 為兩年內發生近視度數回歸（myopia regression），可視為手術失敗，0 則為手術成功。

在筆者之前的文章中，曾經介紹線性混合模式（Linear mixed model, LMM）的概念與應用場景（請見：https://reurl.cc/8GdGr4、https://reurl.cc/pdOdqx、https://reurl.cc/20p0xE或是https://reurl.cc/WdlL3y），在本系列文章中，將介紹以 SPSS 25 版進行實際的操作示範與結果解釋。

LMM 可被應用於分析多層次資料（Multilevel data），因此也稱為多層次模型/迴歸（Multilevel model / regression）。所謂多層次的定義，只要資料結構為至少「2 層」以上即可稱為多層次的資料。下圖示範一個「4 層」的資料結構，第四層為醫院，第三層為每家醫院的醫生、第二層為每位醫生治療的病人，第一層為每位病人的多筆資料（例如同一位個案的左右眼、多個不同病灶、重複測量等）。

在本系列前一篇文章中，筆者介紹了傳統方法的基本作法，接著在本文中將介紹當代的主流方法（最主要會介紹兩種：雙變量模型與階層模型），以及說明這些主流方法可以克服哪些傳統作法不合理的假設。

統合分析（meta-analysis）普遍被應用在各種學科與領域，然而其中診斷型統合分析（meta-analysis for diagnosis）這個分支，無論在方法學的開展或普遍性的應用都是以醫學領域為主，主要用在於一個「標記」（marker）對於一個疾病或狀態（通常是二元變項）的區別/鑑別能力（discrimination）的評估。

下圖列出在 PubMed 搜尋診斷型統合分析的相關論文，很明顯發表的文章越來越多，顯示此方法的重要性。

在各類的醫學量化研究中，資料集中含有遺漏值（missing value）幾乎是無可避免的狀況。關於遺漏值的機制（mechanism），可分成三類：（一）Missing completely at random (MCAR)、（二）Missing at random (MAR) 以及（三）Missing not at random（MNAR）。

八、使用「rms」得到 Resampling calibration 及 Validation 指標

在「rms package」中，只提供 k-fold 跟 bootstrap 兩種驗證模式的方式，首先我們先來看 bootstrap validation 的指令非常的簡單，首先先固定種子數，否則每一次重複抽樣的樣本組合都會不一樣，於是就會得到不一樣的結果。

                  欲驗證的是 fit2 這個物件，亦即包含「Proteinuria」的完整模型；B=2000 代表一共要 2000 個 bootstrap samples。除非是樣本數很大（例如數萬甚至數十萬），否則建議 bootstrap samples 設至少一千次。

七、Validation 原理說明

在現今以預測模型做為主要目標的論文，幾乎都要再報告 Resampling calibration 及 Validation 的結果，才有機會可以通過統計審查這一關，特別是發在優質的雜誌時。

雖然在上述第五點，雖然有提到以「PredictABEL」得到一些其他 Calibration 指標，但須注意，這些指標都是利用該模式的「預測機率」算出來的。然而其實這個模式的「類推性」（generalizability）是受限制的，因為一個模式通常都是在我們同一群樣本之下，用精挑細選的方式揀擇最後留在模式中的變項，因此通常會傾向得到一個高估模式表現的結果（例如比較高的解釋力、比較低的預測誤差）。

五、使用「PredictABEL」得到其他 Calibration 及 Discrimination 指標

             在「PredictABEL」套件中，有許多關於 Calibration 及 Discrimination 的指標。下方為常用的 Hosmer-Lemeshow goodness-of-fit，其中「predRisk= predRisk2」使用完整模型（預測變項包括Proteinuria + SOFA + Age + Male + DM + SBP）的預測機率。

三、使用「pROC」比較兩個預測變項的曲線下面積

接著安裝並呼叫「pROC」（compare ROC）套件，由於更換了套件，因此此時要再 attach 一次資料集「proteinuria」。

再來創造「roc1」物件，括號前面的「proteinuria$AKI」表示資料集「proteinuria」中的「AKI」這個變項是結果變項；括號後面的「proteinuria$SOFA」則是資料集「proteinuria」中的「SOFA」這個變項是預測變項。

筆者曾在之前的一系列文章，介紹風險預測模型（risk prediction model）使用 area under the curve（AUC）的限制（https://reurl.cc/GkX1Vy），以及應該如何評估預測模型的表現（https://reurl.cc/RdON4x、https://reurl.cc/GkX1rZ、https://reurl.cc/W4qnLL），進而介紹預測模型的新指標－IDI 與 NRI（https://reurl.cc/M7Xedk、https://reurl.cc/qDZE8N），本系列文章則以 R 軟體進行操作示範。

一、安裝軟體及套件，並匯入資料

首先讀者可先自行安裝最新版 R 軟體（https://cran.r-project.org/bin/windows/base/），並且安裝操作介面較為友善的免費版 R studio（https://rstudio.com/products/rstudio/download/）。接著在 R studio 右下角的『Install』點擊一下，並陸續安裝以下套件（package）：『rms』、『Hmisc』、『PredictABEL』以及『pROC』。