《Introductory Statistics With R》閱讀筆記 ????????
壽命分析是生物學(xué)和醫(yī)學(xué)領(lǐng)域的重要課題，尤其是在工程應(yīng)用中的可靠性分析中。這樣的數(shù)據(jù)通常高度非正態(tài)分布，因此使用標(biāo)準(zhǔn)線性模型是有問題的。生命周期數(shù)據(jù)經(jīng)常被審查（censored），因為不知道確切的生命周期，只知道它比給定的值長。例如，在一項癌癥試驗中，有些人無法進(jìn)行隨訪，或者干脆活過了研究期。在統(tǒng)計分析中忽略審查是錯誤的，有時會造成極端的后果。
12.1 Essential concepts
- X :true lifetime
-
T :censoring time
T可以是一個隨機變量，也可以是一個固定時間，這取決于上下文，但如果它是隨機的，那么它通常應(yīng)該是非信息的，我們在這里描述的方法是適用的。
有時，“其他原因造成的死亡”被認(rèn)為是給定疾病死亡率的censoring event，在這些情況下，確保這些其他原因與疾病狀態(tài)無關(guān)尤為重要。 -
S(t):survival function
測量在給定時間存活的概率。 實際上，它只是1減去X的累積分布函數(shù)，1- F(t)。 -
h(t):hazard function or force of mortality
假設(shè)被試在t時刻還活著，測量在短時間t內(nèi)死亡的(無窮小的)風(fēng)險。 -
如果壽命分布具有密度f，則h(t) = f(t)/S(t)
這通常被認(rèn)為是比(例如)生存分布的平均中值更基本的數(shù)量，并被用作建模的基礎(chǔ)。
12.2 Survival objects
survival
這個包實現(xiàn)了大量的先進(jìn)技術(shù)。對于目前的目的，只使用它的一個小子集。
> library(survival) #載入survival
- 生存與Surv類的對象一起工作，該類是結(jié)合時間和檢查信息的數(shù)據(jù)結(jié)構(gòu)。
- 類對象使用Surv函數(shù)構(gòu)造，帶兩個參數(shù)：observation time和event indicator。后者可以編碼為邏輯變量，0/1變量或1/2變量。實際上，Surv還可以與三個參數(shù)一起使用，用于處理具有開始時間和結(jié)束時間的數(shù)據(jù)(“交錯輸入”)和間隔審查數(shù)據(jù)(知道事件發(fā)生在兩個日期之間，例如在疾病的重復(fù)測試中)。
> library(ISwR)
> data(melanom) #載入生存數(shù)據(jù)
> attach(melanom) #把數(shù)據(jù)集放在檢索路徑上
> names(melanom)
[1] "no" "status" "days" "ulc" "thick" "sex"
> head(melanom)
no status days ulc thick sex
1 789 3 10 1 676 2
2 13 3 30 2 65 2
3 97 2 35 2 134 2
- status變量是研究結(jié)束時患者狀態(tài)的一個指標(biāo):1表示“死于惡性黑色素瘤”，2表示“在1978年1月1日還活著”，3表示“死于其他原因”；
- days為觀察時間，以天為單位；
- ulc為(1/2)有/無腫瘤有無潰瘍；
- thick為厚度,單位1/100 mm；
- sex為患者性別，1表示女性，2表示男性。
#創(chuàng)建一個Surv對象，將狀態(tài)變量的值2和3視為檢查
> Surv(days, status==1)
[1] 10+ 30+ 35+ 99+ 185 204 210 232 232+ 279 295 355+
[13] 386 426 469 493+ 529 621 629 659 667 718 752 779
- 帶有“ +”標(biāo)記檢查觀察結(jié)果。
- 注意：帶有“+”的是狀態(tài)1的，剩下的是狀態(tài)2和3的。
12.3 Kaplan–Meier estimates
Kaplan-Meier estimator 允許在右截尾情況下計算估計的生存函數(shù)。它也被稱為product-limit estimator（乘積極限估計），因為描述這個過程的一種方法是，它將沒有截尾觀察或沒有死亡的區(qū)間內(nèi)的條件生存曲線相乘。
使用稱為 survfit 的函數(shù)來計算生存函數(shù)的 Kaplan-Meier 估計量。它只需要一個參數(shù)，即Surv對象。 它返回一個survfit對象。 如上所述，我們認(rèn)為“其他原因?qū)е碌乃劳觥笔且环N檢查，并執(zhí)行以下操作：
> survfit(Surv(days,status == 1)~1)
Call: survfit(formula = Surv(days, status == 1) ~ 1)
n events median 0.95LCL 0.95UCL
205 57 NA NA NA
- 得到一些簡要的統(tǒng)計數(shù)據(jù)和中值生存的估計，在這種情況下，后者甚至都不有趣，因為估計是無限的。
要查看實際的 Kaplan-Meier 估計，請在 survfit 對象上使用 summary。將 survfit 對象保存到一個名為 surv.all 變量中。這是因為它包含了所有病人的原始生存功能，而不考慮病人的特征。
> surv.all <- survfit(Surv(days,status==1)~1)
> summary(surv.all,censored=F)
Call: survfit(formula = Surv(days, status == 1) ~ 1)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
185 201 1 0.995 0.00496 0.985 1.000
204 200 1 0.990 0.00700 0.976 1.000
210 199 1 0.985 0.00855 0.968 1.000
- 它包含生存函數(shù)在事件時刻的值。censored=T 可以顯示審查時間。
- KM是一個階躍函數(shù)，它的跳躍點是在 time 上給出的，跳躍后的值是在 survival 上給出的。另外，給出了曲線的標(biāo)準(zhǔn)誤差估計和真曲線的點態(tài)置信區(qū)間。
> plot(surv.all)
在同一圖上繪制兩個或多個生存函數(shù)通常很有用，以便可以直接比較它們（圖12.2）。 要獲得按性別劃分的生存功能，請執(zhí)行以下操作：
> plot(surv.bysex, conf.int=T, col=c("black","red"))
- conf.int=T 打開置信區(qū)間，如果想要達(dá)到99% 的置信水平，設(shè)置 conf.int=0.99。
- col 給每組賦予不同的顏色。
12.4 The log-rank test
對數(shù)秩檢驗用于檢驗兩條或多條生存曲線是否相同。它是基于觀察每個死亡時間的人口，并根據(jù)每一組中處于危險中的個人的數(shù)量，計算預(yù)期死亡人數(shù)。然后將其與所有死亡時間相加，并通過與χ2檢驗相似（但不相同）的程序與觀察到的死亡人數(shù)進(jìn)行比較。
利用函數(shù) survdiff 計算對數(shù)秩檢驗。這實際上實現(xiàn)了一系列由參數(shù) rho 指定的測試，允許對零假設(shè)進(jìn)行各種非比例的風(fēng)險替代，但是 rho = 0 的默認(rèn)值給出了 log-rank 測試。
> survdiff(Surv(days,status==1)~sex)
Call:
survdiff(formula = Surv(days, status == 1) ~ sex)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 126 28 37.1 2.25 6.47
sex=2 79 29 19.9 4.21 6.47
Chisq= 6.5 on 1 degrees of freedom, p= 0.01
- 該規(guī)范使用了線性和廣義線性模型的模型公式。 然而，這個測試只能處理分組的數(shù)據(jù)，所以如果你在右邊指定多個變量，它將處理所有預(yù)測變量組合產(chǎn)生的數(shù)據(jù)分組。 它也不區(qū)分因子和數(shù)字代碼。 survfit 也是如此。
也可以指定分層分析，即在數(shù)據(jù)集的分層中分別進(jìn)行觀測值和預(yù)測值的計算。 例如，可以計算按潰瘍分層的性別效應(yīng)的對數(shù)等級檢驗如下:
> survdiff(Surv(days,status==1)~sex+strata(ulc))
Call:
survdiff(formula = Surv(days, status == 1) ~ sex + strata(ulc))
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 126 28 34.7 1.28 3.31
sex=2 79 29 22.3 1.99 3.31
Chisq= 3.3 on 1 degrees of freedom, p= 0.07
- 請注意，這會使性的影響顯得不那么顯著。 一種可能的解釋可能是，男性在疾病處于比女性更嚴(yán)重的狀態(tài)時尋求治療，因此當(dāng)根據(jù)疾病進(jìn)展的衡量標(biāo)準(zhǔn)進(jìn)行調(diào)整時，性別差異會縮小。
12.5 The Cox proportional hazards model
作為第一個例子，考慮一個單回歸變量的模型:
> summary(coxph(Surv(days,status==1)~sex))
Call:
coxph(formula = Surv(days, status == 1) ~ sex)
n= 205, number of events= 57
coef exp(coef) se(coef) z Pr(>|z|)
sex 0.6622 1.9390 0.2651 2.498 0.0125 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 1.939 0.5157 1.153 3.26
Concordance= 0.59 (se = 0.034 )
Likelihood ratio test= 6.15 on 1 df, p=0.01
Wald test = 6.24 on 1 df, p=0.01
Score (logrank) test = 6.47 on 1 df, p=0.01
- coef 指兩組之間危險比的估計對數(shù),也將其作為實際危險比 exp(coef) 給出。
- 下面的一行還給出了危險比的倒置比率(交換組)和置信區(qū)間。 最后，對模型中的顯著效應(yīng)進(jìn)行了三次全面的檢驗。
- 這些在大樣本中都是等價的，但在小樣本中可能有所不同。 注意，Wald 檢驗與基于估計系數(shù)除以標(biāo)準(zhǔn)誤差的 z 檢驗相同，而得分檢驗等同于對數(shù)秩檢驗(只要模型只涉及一個簡單的分組)。
更詳細(xì)的例子，涉及連續(xù)協(xié)變量和分層變量：
> summary(coxph(Surv(days,status==1)~sex+log(thick)+strata(ulc)))
Call:
coxph(formula = Surv(days, status == 1) ~ sex + log(thick) +
strata(ulc))
n= 205, number of events= 57
coef exp(coef) se(coef) z Pr(>|z|)
sex 0.3600 1.4333 0.2702 1.332 0.1828
log(thick) 0.5599 1.7505 0.1784 3.139 0.0017 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 1.433 0.6977 0.844 2.434
log(thick) 1.750 0.5713 1.234 2.483
Concordance= 0.673 (se = 0.039 )
Likelihood ratio test= 13.3 on 2 df, p=0.001
Wald test = 12.88 on 2 df, p=0.002
Score (logrank) test = 12.98 on 2 df, p=0.002
- 可以看出，性別變量的重要性已經(jīng)進(jìn)一步降低。
Cox模型假設(shè)了一個潛在的基線危險函數(shù)，以及一個對應(yīng)的生存曲線。在分層分析中，每一層都有一條相同的曲線。它們可以通過使用 coxph 輸出的 survfit 來提取，當(dāng)然也可以使用針對 survfit 對象的 plot 方法來繪制(圖12.3):
> plot(survfit(coxph(Surv(days,status==1)~log(thick)+sex+strata(ulc))))
- 注意，survfit 的默認(rèn)設(shè)置是為協(xié)變量在其平均值處的假個體生成曲線。 本例腫瘤厚度為1.86 mm，性別為1.39(!) 。
- 可以使用 survfit 的 newdata 參數(shù)來指定要為其計算生存曲線的數(shù)據(jù)框。
