[統合分析] 如何生成single arm meta-analysis 的圖表

在統合分析裡面, 通常會有兩個arm, 不過有的時候會需要對症狀的發生率, 疾病的發生率進行統合, 此時就需要做所謂的single arm meta-analysis, R 裡面有非常好用的套件, 以下就由我來示範!

首先假設假想例子:

Study 1: 30個病人有5個中風, 15男有2個中風, 15女有3個中風

Study 2: 40個病人有7個中風, 20男有3個中風, 20女有4個中風

Study 3: 20個病人有9個中風, 10男有4個中風, 10女有5個中風

Study 4: 70個病人有10個中風, 30男有5個中風, 40女有5個中風

則可以使用R 的 metaprop 函式來生成森林圖(forest plot), 其程式碼如下:

library(meta)
library(metafor)
library("grid")

studies <- c("Study 1","Study 2","Study 3","Study 4")
obs <- c(5,7,9,10)
denom <- c(30,40,20,70)


m1 <- metaprop(obs, denom, studies, comb.fixed=T) #com.fixed 表示fixed effect model


forest(m1, leftlabs=c("Study","Stroke event","Total"), 
       ref=0.2, colgap=unit(3,"mm"), col.by="black",fontsize=10,
       rowgap=unit(3,"mm"))

生成圖片如下:

                                                                                                                                           

其中圖裡面的 I^2 就是統合分析中拿來評估異質性(heterogenity)的指標, 值越低表示異質性越低, 意指為統合分析中研究間的差異性越小

另外就是統合分析裡面, 經常要看所謂的publication bias, 程式碼也是非常簡單, 如下

funnel(m1)

metabias(m1) #Egger's test

如圖可以看到, 明顯兩邊不對稱, 雖然這麼少的文獻要看publication bias不太夠(依照目前cochrane handbook的建議, 至少要十篇論文來生成funnel plot會比較有意義), 不過一般來說不對稱就代表可能有publication bias!

subgroup analysis, 一般中文翻作次群體分析, 以這個例子為例, 上面有提供男女資訊, 那我們就會想知道男女是否會有差異, 一樣用程式碼做示範,

studies2 <- c("Study 1","Study 2","Study 3","Study 4",
              "Study 1","Study 2","Study 3","Study 4")

obs2 <- c(2,3,4,5,3,4,5,5)

denom2 <- c(15,20,10,30,15,20,10,40)

setting<-c("M", "M", "M","M",
           "F", "F", "F", "F")

m2 <- metaprop(obs2, denom2, studies2, comb.fixed=F,byvar=setting,
) #com.fixed 表示fixed effect model

forest(m2, leftlabs=c("Study","Stroke event","Total"), 
       ref=0.1, colgap=unit(5,"mm"), col.by="black",fontsize=10,
       rowgap=unit(3,"mm"))

之後再利用meta-regression, 替男女兩組做比較, 看是否有差異!

df1 <- structure(list(study_nr = c(1,2,3,4,1,2,3,4),
                     group = c("A", "A", "A", "A",
                               "B", "B", "B", "B"), 
                     n_tot = c(15,20,10,30,15,20,10,40), 
                     n_event = c(2,3,4,5,3,4,5,5)), 
                class = c("tbl_df", "tbl", "data.frame"),
                row.names = c(NA, -8L))
df1

# 計算 log(odds) 以及 corresponding sampling variances
df1 <- escalc(measure="PLO", xi=n_event, ni=n_tot, data=df1)
df1

# fit meta-regression model with 'group' as moderator
res1 <- rma(yi, vi, mods = ~ group, data=df1, method="DL")
res1

# 得到 odd ratio 和 95% CI
predict(res1, newmods=1, intercept=FALSE, transf=exp, digits=2)

Test for Residual Heterogeneity:

QE(df = 6) = 9.3173, p-val = 0.1565

p value >=0.05, 沒有顯著意義!

reference: 

1. metafor package

2. stackoverflow

3. https://rpubs.com/pekong/532068

[攝護腺癌] 術後放療 versus 救援性放療 -1 (Radicals-RT trial)

要現在做放射治療呢? 還是之後PSA高起來再做呢? 好難啊   (圖片取自吉卜力工作室)

攝護腺癌手術後如果有危險因子的話, 應該直接做術後輔助放療(adjuvant RT)或是等到之後PSA 高起來再做救援性放療(salvage RT), 一直是一個臨床上難解的問題.

在西方國家, 由於飲食的關係, 攝護腺癌的人數非常多, 這次的實驗族群主要是位於丹麥, 英國, 愛爾蘭以及加拿大, 共收案1396個病人, 為第三期隨機分派試驗(RCT), 其收案條件為接受攝護腺手術(radical prostectomy), 以及有一個以上的危險因子(pT3-4, Gleason score of 7–10, positive margins, or 手術前 PSA ≥10 ng/mL)

Biochemical failure 定義成 連續兩次PSA >= 0.1ng/ml 或是連續三次上升

放射治療技術有兩種, 52.5Gy/20Fr or 66Gy/33Fr（prostate RT +/-WPRT）

primary endpoints 定義為 freedom from distant metastases;

Secondary outcome 有看存活分析, 治療毒性, 開始使用賀爾蒙(non-protocol use),或是病人自述結果, 文章裡面有說本來是沒有要看biochemical progression的, 後來為了跟下面P.S. 提到的兩位兄弟一起做統合, 所以才列進來, 順帶一提這篇收案一千三百多人, 是三兄弟裡面最多

P.S. 其實這個主題目前至少還有兩篇trials, 不過primary endpoint 都是看biochemical progression, 分別是  RAVES and GETUG-AFU 1, 之後再來整理! 這些文章最後都變成一篇統合, 叫做ARTISTIC meta-analysis, 也是一樣之後來介紹! 這裡偷偷劇透一下, 統合分析的結果出來也是支持early salvage RT, 所以目前證據看起來early salvage RT在臨床上比較可行!

結果:

Adjuvant RT 和 salvage RT 在 5-year biochemical progression-free survival 跟 Freedom from non-protocol hormone therapy at 5 years 沒有統計上差別, 副作用方面adjuvant RT 還比較多

reference:

1. Parker, Christopher C., et al. "Timing of radiotherapy after radical prostatectomy (RADICALS-RT): a randomised, controlled phase 3 trial." The Lancet (2020).

2. Kneebone, A., et al. "A Phase III Multi-Centre Randomised Trial comparing adjuvant versus early salvage Radiotherapy following a Radical Prostatectomy: Results of the TROG 08.03 and ANZUP “RAVES” Trial." International Journal of Radiation Oncology• Biology• Physics 105.1 (2019): S37-S38.

3. Carrie C, Magné N, Burban-Provost P, et al. Short-term androgen deprivation therapy combined with radiotherapy as salvage treatment after radical prostatectomy for prostate cancer (GETUG-AFU 16): a 112-month follow-up of a phase 3, randomised trial. Lancet Oncol. 2019;20(12):1740-1749. doi:10.1016/S1470-2045(19)30486-3

[輻射生物] 輻射傷害機轉

在輻射生物學裡面, 有一個很重要個概念就是輻射造成的傷害, 應該怎樣來做操作型定義(operational definition), 一般輻射生物的書會分成三種, 以下一一來介紹:

1. lethal damage: 直接導致細胞死亡

2. potential lethal damage: 照射後細胞是否死亡要看照射後的環境, 如果照射後把細胞放在抑制生長的環境, 一段時間後再把細胞種在充滿養分適合生長的培養皿, 則細胞比較不會死亡, 如果直接把照射後的細胞種在充滿養分的培養皿, 則細胞容易死亡, 因為在充滿養分的培養皿, 細胞容易進行分裂, 則損傷就無法修補, 進而死亡!

3. sublethal damage: 如果把劑量拆成兩等分,而且中間隔一段時間, 則存活率會上升! 因為中間間格的時間會進行所謂的sublethal damage repair

了解了輻射傷害的定義後, 在來就是可以針對所謂的劑量率效應(dose-rate effect)來做描述, dose rate effect 的定義是當劑量率位在 1 Gy per minute 到 0.3 Gy per hour 的時候, 且暴露輻射的時間足夠長的時候, sublethal damage repair 的效應會出現, 因此存活率的圖可以看成無數個初始線段做相加! 

這裡可以稍微注意一下劑量率效應在臨床上的應用, 我們一般子宮頸癌近接治療通常是在十分鐘左右給腫瘤的GTV 大概 650cGy, 換算下來大概, 0.65Gy per minute 左右, 因此也會有所謂的劑量率效應(dose-rate effect)

ICRU 38 報告定義 HDR(high dose rate) 必須每分鐘至少給20cGy, 換算下來至少要 12Gy/Hr, HDR 跟LDR相比, 其針對不同cell line 的作用比較一致, 因此臨床上適合用來對付輻射抗性(radio-resistant) 較強的細胞

reference:

1. Radiobiology Eric and Hall 8th  

2. Khan's physics of radiation therapy 6th

[人工智慧] 如何比較模型間的預測能力

在人工智慧的論文裡面, 最重要的就是比較模型之間的預測能力, 模型的預測能力常用所謂的ROC-AUC (Area under curve) 來做表示, 一提到AUC, 有修過統計的應該會有模糊的印象吧(笑), 簡單來說就是利用調整threshold, 來構做ROC(Receiver operating Curve), 並求取底下的面積, 以下為示範案例!

因為本身是癌症科醫師, 因此當然使用癌症相關資料集, 以下用R語言來做示範!

library(xlsx)
library("glmnet")
library(pROC)
library("ggcorrplot")
library("readxl")
library("randomForest")
library(mlbench)
library(Hmisc)

data("BreastCancer")
cohort <- BreastCancer #cohort from NCI library
cohort <- as.data.frame(cohort)      

先來看一下資料分布, 其資料內部主要是乳房腫瘤的切片資料, 以及最後是良性或是惡性, 在許多人工智慧場合, 理論事都應該先做所謂的EDA(exploratory data analysis)來看資料分布, 不過很多人都沒做XD

這裡主要是預測二元資料(binary data), 但是如果多考慮時間的話, 在生物統計上經常性使用所謂的存活分析(survival analysis), 其中Cox proportional hazards model 經常性地被使用, 評估Cox model 的預測能力時, 常用所謂Concordance index, 其發明者就是就是Hmisc包的作者, Frank Harrell, 所以Concordance index 又被叫做 Harrell's concordance index.

這裡借用Hmisc 包的 describe 函數來做EDA

describe(BreastCancer)

其結果如下:

這張圖縱軸是代表每個項目,  橫軸是把連續值分成十個部分, 每個部份去看分布, 做完EDA後, 可以知道Bare Nuclei 的缺失值比較多, 為了避免後續處理出問題, 因此進行缺失值處理

 

cohort <- cohort[complete.cases(cohort), ] # remove NA value

人工智慧最重要的部分就是把資料分成所謂的training set 以及 testing set, 隨機分法非常多種, 這裡就不贅述, 以下是程式碼

sample_size = floor(0.8*nrow(cohort)) 

train_index = sample(seq_len(nrow(cohort)),size = sample_size) 
train =cohort[train_index,] 
test=cohort[-train_index,]

true_values <- as.factor(ifelse(testy=="malignant",1,0))
true_valuestr <- as.factor(ifelse(trainy=="malignant",1,0))

上面這段程式碼就是把資料分成training set: 八成, testing set: 兩成; 以及把良惡性的結果變成factor

之後利用logistic regression 來看跑看看模型:

logitMod <- glm(Class~ Cl.thickness+Cell.size+Cell.shape+Marg.adhesion
                +Epith.c.size+Bare.nuclei
                , data=train, family=binomial(link="logit"))
summary(logitMod)

final_y <- predict(logitMod,testx,type="response")
roc1<-roc(true_values,final_y,direction = "<")
ggroc(roc1,legacy.axes="TRUE")
auc1<-auc(roc1)


final_ytr <- predict(logitMod,trainx,type="response")
roctr<-roc(true_valuestr,final_ytr,direction = "<")
ggroc(roctr,legacy.axes="TRUE")
auctr1<-auc(roctr)

下面是先放testing cohort 的 ROC curve (AUC= 0.907)

可以注意到橫軸是1-specificity, 縱軸是所謂的sensitivity, 這張圖就是利用調整threshold, 來改變sensitivity以及spectivity來完成的

接下來是看training cohort (AUC=1)

雖然logistic regression 表現已經相當驚人(AUC =0.907), 但還是想跟別的模型比較, 因此使用常見的機器學習模型, 隨機森林, random forest 來做示範, 程式碼如下:

#random forest
rf <- randomForest(Class~ Cl.thickness+Cell.size+Cell.shape+Marg.adhesion
                   +Epith.c.size+Bare.nuclei
                   , data=train,ntree=50,sampsize=100)

final_yrf <- predict(rf,testx,type="prob")
roc1rf<-roc(true_values,final_yrf[,2],direction = "<")
ggroc(roc1rf,legacy.axes="TRUE")
aucrf<-auc(roc1rf)


final_ytrrf <- predict(rf,trainx,type="prob")
roctrrf<-roc(true_valuestr,final_ytrrf[,2],direction = "<")
ggroc(roctrrf,legacy.axes="TRUE")
auctrrf<-auc(roctrrf)

下面還是先放模型在 testing cohort 的表現 (AUC=0.987)

在下面才放模型在training cohort的表現(AUC=0.997)

最後用Delong test 對兩個模型在testing cohort的表現進行比較 

roc.test(roc1,roc1rf,method="delong")

p-value = 0.005132,明顯random forest表現比較好
 

reference:

1. https://rpubs.com/skydome20/R-Note10-Missing_Value