Curs3: Compararea probabilitatilor de supravietuire intre doua grupuri
Conf.dr. Maria Jivulescu
Breviar teoretic*
Ne propunem compararea eficacitatii unui medicament, prin compararea probabilitii de supravietuire dintre cele doua grupe urmarite intr-un studiu clinic: grupul de tratament si grupul de placebo.
Vom urmari ca obiectiv compararea celor doua curbe de suparvietuire, penru a observa diferentele, raspunzand la intrebarea:
Cele doua functii de supravietuire sunt diferite din punct de vedere statistic?
Vom formula urmatoarele ipoteze statistice: \[H_0: S_1(t)=S_2(t), \text{pentru toti t}\] \[H_1: S_1(t)\neq S_2(t), \text{pentru anumiti t}\]
Pentru a testa aceste ipoteze se foloseste testul Log-Rank, ce consta in parcurgerea urmatorilor pasi:
se noteaza \(t_1<t_2<\ldots <t_k\) momentele de timp in care are loc evenimentul urmarit;
\(d_{1i}\) si \(d_{2i}\) numarul de persoane decedate la momentul \(t_i\) din grupul 1, respectiv 2; avem \(d_i=d_{1i}+d_{2i}\)
\(n_{1i}\) si \(n_{2i}\) numarul de persoane aflate la risc la momentul \(t_i\) din grupul 1, respectiv 2; avem \(n_i=n_{1i}+n_{2i};\)
Situatia la momentul \(t_i\) se descrie prin urmatorul tabel de contingenta;
as.table(matrix(c(1,"d_{1i}","n_{1i}-d_{1i}","n_{1i}",2,"d_{2i}","n_{2i}-d_{2i}","n_{2i}","Total","d_{i}","n_{i}-d_{i}","n_{i}"), nrow=3, byrow=TRUE, dimnames=list(c(),
c("Grup", "Decedati", "Supravietuitori", "Total"))))## Grup Decedati Supravietuitori Total
## A 1 d_{1i} n_{1i}-d_{1i} n_{1i}
## B 2 d_{2i} n_{2i}-d_{2i} n_{2i}
## C Total d_{i} n_{i}-d_{i} n_{i}Observatii:
- Sub ipoteza \(H_0\) variabila \(d_{1i}\) este o variabila distribuita hipergeometric de parametrii \(n_i\), \(n_{1i}\) si \(d_i\).
-\(E(d_{1i})=\frac{n_{1i}d_i}{n_i}\)
-\(\sigma^2(d_{1i})=\frac{n_{1i}n_{2i}(n_i-d_i)d_i}{n_i^2(n_i-1)}\)
-Se calculeaz o statistica \(U=\sum\limits_{i=1}^k(d_{1i}-E(d_{1i}))\)
-Avem \(E(U)=0, \sigma^2(U)=\sum\limits_{i=1}^k\frac{n_{1i}n_{2i}(n_i-d_i)d_i}{n_i^2(n_i-1)}\)
-prin standardizare se obine \[Z=\frac{U}{\sqrt{\sigma^2(U)}}\sim N(0,1)\]
-Alternativ, \(Z^2\sim \chi^2(1)\)
Aplicatia1 (problema 3.4, pag 41)
Un studiu clinic isi propune studiul eficacitatii unui nou tip de plasture folosit impotriva dependentei de nicotina. Subiectii sunt aleator impartiti in grupul de tratament (cu noul tip de plasture) si control (cu cel mai bun produs de acelasi tip). Se masoara (in luni) timpul in care o persoana din studiu il petrece fara consum de tutun. Datele inregistrate sunt:
Tratament: 3.4, 3.6+, 4.1, 4.9+, 5.8+
Control: 2.0, 3.7+, 4.3, 4.9+
Rezolvare analitica
Se ordoneaza crescator momentele de timp in care a avut loc evenimentul urmarit: \[2<3.4<4.1<4.3\] Pentru fiecare moment de timp \(t_i\) se va centraliza situatia evenimentelor ce au loc, intr-un tabel de contingenta (asa cum apare in Breviarul teoretic).
Astfel, la momentul \(t_1=2\), avem
as.table(matrix(c(1,0,5,5,2,1,3,4,"Total",1,8,9), nrow=3, byrow=TRUE, dimnames=list(c(),
c("Grup", "Decedati", "Supravietuitori", "Total"))))## Grup Decedati Supravietuitori Total
## A 1 0 5 5
## B 2 1 3 4
## C Total 1 8 9Asadar, \(d_{11}=0, E(d_{11})=\frac{5\cdot 1}{9}=\frac{5}{9}, \sigma^2(d_{11})=\frac{20}{81}\).
La momentul \(t_2=3.4\) situatia este:
as.table(matrix(c(1,1,4,5,2,0,3,3,"Total",1,7,8), nrow=3, byrow=TRUE, dimnames=list(c(),
c("Grup", "Decedati", "Supravietuitori", "Total"))))## Grup Decedati Supravietuitori Total
## A 1 1 4 5
## B 2 0 3 3
## C Total 1 7 8Avem \(d_{12}=1, E(d_{12})=\frac{5\cdot 1}{8}=\frac{5}{8}, \sigma^2(d_{12})=\frac{15}{64}\).
La momentul \(t_3=4.1\) situatia este:
as.table(matrix(c(1,1,2,3,2,0,2,2,"Total",1,4,5), nrow=3, byrow=TRUE, dimnames=list(c(),
c("Grup", "Decedati", "Supravietuitori", "Total"))))## Grup Decedati Supravietuitori Total
## A 1 1 2 3
## B 2 0 2 2
## C Total 1 4 5\(d_{13}=1, E(d_{13})=\frac{3\cdot 1}{5}=\frac{3}{5}, \sigma^2(d_{13})=\frac{6}{25}\).
La momentul \(t_4=4.3\) situatia este:
as.table(matrix(c(1,0,2,2,2,1,1,2,"Total",1,3,4), nrow=3, byrow=TRUE, dimnames=list(c(),
c("Grup", "Decedati", "Supravietuitori", "Total"))))## Grup Decedati Supravietuitori Total
## A 1 0 2 2
## B 2 1 1 2
## C Total 1 3 4Avem \(d_{14}=0, E(d_{14})=\frac{2\cdot 1}{4}=\frac{1}{2}, \sigma^2(d_{14})=\frac{1}{4}\).
Pentru a lua o decizie legata de testul facut, vom calcula
\(U=\sum\limits_{i=1}^4(d_{1i}-E(d_{1i}))=-0.2806\)
\(\sigma^2(U)=0.9713\)
\(z=-\frac{0.2806}{\sqrt{0.9713}}=-0.2847\).
\(2P(Z>0.2847)=0.7759\) sau \(z^2=0.081\) i \(P(\chi^2(1)>0.081)=0.7759\)
Deci, $H_0$ nu se respinge la un nivel de semnificatie de $0.05$.
Concluzie: Cele doua curbe nu difera semnificativ din punct de vedere statistic.
Rezolvare in R
Se incarca pachetul folosit in Analiza Supravietuirii
library(survival)Vom organiza datele din problema intr-un tabel:
t1<-c(3.4,3.6,4.1,4.9,5.8) #vectorul timpiilor inregistrati, din grupul tratament
t2<-c(2,3.7,4.3,4.9)#vectorul timpiilor inregistrati, din grupul control
t<-c(t1,t2) #vectorul timpiilor inregistrati, din ambele grupuri
s1<-c(1,0,1,0,0) #vectorul tipului de observatie din grupul de trataent
s2<-c(1,0,1,0) #vectorul tipului de observatie din grupul de control
s<-c(s1,s2) #vectorul tipului de observatie
rx<-c(1,1,1,1,1,0,0,0,0) #vectorul de grup (grup1=tratament, sau grupul 2=control)
g<-data.frame(cbind(t,s,rx)) #tabelul de date
gVom face o scurta analiza statistica a timpiilor de suparvietuire dintre cele doua grupuri.
KM1<-survfit(Surv(t1, s1)~1, data=g)
KM1## Call: survfit(formula = Surv(t1, s1) ~ 1, data = g)
##
## n events median 0.95LCL 0.95UCL
## [1,] 5 2 NA 4.1 NAsummary(KM1)## Call: survfit(formula = Surv(t1, s1) ~ 1, data = g)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3.4 5 1 0.800 0.179 0.516 1
## 4.1 3 1 0.533 0.248 0.214 1Se observa ca nu se poate determina mediana timpului de supravietuire in grupul de tratament, conform rezultatelor probabilitatile de supravietuire sunt peste 0.5, indiferent de momentul de timp.
KM2<-survfit(Surv(t2, s2)~1, data=g)
KM2## Call: survfit(formula = Surv(t2, s2) ~ 1, data = g)
##
## n events median 0.95LCL 0.95UCL
## [1,] 4 2 4.3 2 NAsummary(KM2)## Call: survfit(formula = Surv(t2, s2) ~ 1, data = g)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2.0 4 1 0.750 0.217 0.4259 1
## 4.3 2 1 0.375 0.286 0.0839 1Pentru grupul de controla timpul median are valoarea \(t_m=4.3\).
Vom reprezenta curbele de supravietuire pentru cele doua grupuri.
KM<-survfit(Surv(t, s)~rx, data=g)
KM## Call: survfit(formula = Surv(t, s) ~ rx, data = g)
##
## n events median 0.95LCL 0.95UCL
## rx=0 4 2 4.3 2.0 NA
## rx=1 5 2 NA 4.1 NA summary(KM)## Call: survfit(formula = Surv(t, s) ~ rx, data = g)
##
## rx=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2.0 4 1 0.750 0.217 0.4259 1
## 4.3 2 1 0.375 0.286 0.0839 1
##
## rx=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3.4 5 1 0.800 0.179 0.516 1
## 4.1 3 1 0.533 0.248 0.214 1plot(KM,lty=1:2,mark.time=TRUE, main="Curbe de supravietuire", col="red",
xlab="Timpi inregistrati",
ylab="Probabilitatea de supravietuire")
abline(h=0.5)
Din aceste date se observa ca timpul de supravietuire este mai bun in grupul de tratament(1). Pentru a face o diferentiere statistica justificata riguros, vom aplica un test de diferentiere, cum ar fi testul Log-Rank.
Vom aplica testul Log-Rank pentru a testa diferenta dintre cele doua tipuri de tratament. Pentru acest lucru vom folosi comanda survdiff
survdiff(Surv(t, s)~rx, data=g)## Call:
## survdiff(formula = Surv(t, s) ~ rx, data = g)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=0 4 2 1.72 0.0458 0.081
## rx=1 5 2 2.28 0.0345 0.081
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.8Se observa ca se obtine p-value=0.8, deci nu exista diferenta semnificativa intre cele doua curbe de supravietuire.
Generalizarea teoriei
- situatii cand sunt mai multe grupuri;
- teste alternative: Wilcoxen, Tarrone-Ware, Peto, Flemington-Harrinton;
- aceste teste difera prin faptul ca pun anumite “greutati” pe diverse momente de timp;
Testul Wilcox de diferentiere dintre cele doua curbe de supravietuire
wilcox.test(rx, data=g)## Warning in wilcox.test.default(rx, data = g): cannot compute exact p-value with
## ties## Warning in wilcox.test.default(rx, data = g): cannot compute exact p-value with
## zeroes##
## Wilcoxon signed rank test with continuity correction
##
## data: rx
## V = 15, p-value = 0.03689
## alternative hypothesis: true location is not equal to 0Aplicatia2 (problema 2, pag 75, Kleibaum)
Setul de date Anderson descrie inregistrarile legate de timpul pana la remisie a 42 de pacienti bolnavi de leucemie. Mai multe info despre acest set de date la
https://cardiomoon.github.io/autoReg/reference/anderson.html
Ne prropunem sa descriem curbele KM ale variabilei logwbc= log al nr de celule albe. Variabila logwbc fiind continua se va imparti in 3 categorii, dupa operatia de logaritmare. Astfel, se definesc urmatoarele categorii:low, medium, and high, dupa cum urmeaza: low (0–2.30), n = 11;
medium (2.31–3.00), n = 14; high (>3.00), n = 17. Sa se studieze diferenta dintre cele trei curbe KM asociate celor trei categorii.
Rezolvare in R
Vom apela setul de date anderson (salvat pe calculator)
q<- read.table("anderson2.txt",header=TRUE)Vom crea cele 3 categorii din variabila \(log_{WBC}\), conform cerintei problemei.
q$log_WBCcat[0<q$log_WBC & q$log_WBC<=2.3] <- " 1low"
q$log_WBCcat[2.31<=q$log_WBC & q$log_WBC<=3] <- "2medium"
q$log_WBCcat[q$log_WBC>3] <- "3high"Vom repreeznta garfic ceel 3 curbe de supravietuire, pentru fiecare categorie a variabilei \(log_{WBC}\)
KM<-survfit(Surv(st, status)~q$log_WBCcat, data=q)
KM## Call: survfit(formula = Surv(st, status) ~ q$log_WBCcat, data = q)
##
## n events median 0.95LCL 0.95UCL
## q$log_WBCcat= 1low 11 4 NA 23 NA
## q$log_WBCcat=2medium 14 10 17 10 NA
## q$log_WBCcat=3high 17 16 6 5 11summary(KM)## Call: survfit(formula = Surv(st, status) ~ q$log_WBCcat, data = q)
##
## q$log_WBCcat= 1low
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 11 1 0.909 0.0867 0.754 1
## 12 10 1 0.818 0.1163 0.619 1
## 15 9 1 0.727 0.1343 0.506 1
## 23 5 1 0.582 0.1687 0.330 1
##
## q$log_WBCcat=2medium
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 1 0.929 0.0688 0.8030 1.000
## 4 13 1 0.857 0.0935 0.6921 1.000
## 6 12 1 0.786 0.1097 0.5977 1.000
## 8 11 1 0.714 0.1207 0.5129 0.995
## 10 9 1 0.635 0.1308 0.4240 0.951
## 13 6 1 0.529 0.1457 0.3085 0.908
## 17 5 1 0.423 0.1501 0.2112 0.848
## 22 4 2 0.212 0.1297 0.0636 0.704
## 23 2 1 0.106 0.0990 0.0169 0.662
##
## q$log_WBCcat=3high
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 17 1 0.9412 0.0571 0.8357 1.000
## 2 16 2 0.8235 0.0925 0.6609 1.000
## 3 14 1 0.7647 0.1029 0.5875 0.995
## 4 13 1 0.7059 0.1105 0.5194 0.959
## 5 12 2 0.5882 0.1194 0.3952 0.876
## 6 10 2 0.4706 0.1211 0.2842 0.779
## 7 7 1 0.4034 0.1210 0.2241 0.726
## 8 6 3 0.2017 0.1022 0.0747 0.544
## 11 3 1 0.1345 0.0875 0.0376 0.481
## 12 2 1 0.0672 0.0646 0.0102 0.442
## 16 1 1 0.0000 NaN NA NAplot(KM,lty=1:3,mark.time=TRUE, main="Curbe de supravietuire",
xlab="Timpi inregistrati",
ylab="Probabilitatea de supravietuire")
Vom face testul de diferentiere, pentru a testa daca exista diferente intre categorii.
survdiff(Surv(st, status)~q$log_WBCcat, data=q) ## Call:
## survdiff(formula = Surv(st, status) ~ q$log_WBCcat, data = q)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## q$log_WBCcat= 1low 11 4 13.06 6.2880 12.7695
## q$log_WBCcat=2medium 14 10 10.72 0.0489 0.0809
## q$log_WBCcat=3high 17 16 6.21 15.4173 23.1040
##
## Chisq= 26.4 on 2 degrees of freedom, p= 2e-06Se observa ca se obtine un p-value foarte apropiat de 0, deci exista diferenta semnificativa intre timpii de supravietuire dintre cele 3 curbe.
wilcox.test(q$log_WBC, data=q)## Warning in wilcox.test.default(q$log_WBC, data = q): cannot compute exact
## p-value with ties##
## Wilcoxon signed rank test with continuity correction
##
## data: q$log_WBC
## V = 903, p-value = 1.707e-08
## alternative hypothesis: true location is not equal to 0- Referinta: Clinical Statistics, Olga Korosteleva