Análise de Sobrevivência Aplicada - 2a. edição
Colosimo e Giolo (2024)
Comandos em Linguagem R
Capítulo 1
Conceitos básicos e exemplos
Seção 1.3.1 - Figura 1.1
options(OutDec = ",") # decimais com vírgula
par(mar=c(5,4,1.5,1)) # c(bottom,left,top,right)
par(mfrow=c(1,2))
id<-c(1:10) # pacientes de 1 a 10
Tseg<-c(25,12,17,15,28,20,13,15,20,17) # tempo de seguimento
Tent<-c(0,2,3,4,6,9,12,12,14,17) # tempo de entrada
cens<-c(1,1,1,1,0,1,1,1,0,0)
plot(Tseg,Tseg*0,pch=" ",ylim=range(c(0,10)), xlim=range(c(0,38)),xlab="Calendário",
ylab="Indivíduos",bty="n",cex.lab=1.2,cex.axis=1.0,xaxt="n")
abline(v=0, lty=3,col="black")
abline(v=34,lty=3,col="black")
abline(v=17,lty=3,col="black")
for(i in 1:10){ lines(c(Tent[i],(Tent[i]+Tseg[i])),c(id[i],id[i]))
points((Tent[i]+Tseg[i]), id[i], pch=ifelse(cens[i]==1,16,1)) }
axis(1,at=c(0,12,24,36),labels = c("2016","2017","2018","2019"),cex.axis=1.0)
plot(Tseg,Tseg*0,pch=" ",ylim=range(c(0,10)), xlim=range(c(0,35)),
xlab="Tempo (meses)",ylab="Indivíduos",bty="n",cex.lab=1.2,cex.axis=1.0)
for(i in 1:10){
lines(c(0,(Tseg[i])),c(id[i],id[i]))
points((Tseg[i]),id[i],pch=ifelse(cens[i]==1,16,1))
} Seção 1.5.1 - Figura 1.3
par(mfrow=c(1,1))
t<-seq(1,35,0.5)
st1<-exp(-t/14.5)
st2<-exp(-((t/22.5)^3.0))
plot(c(0,t),c(1,st2),xlab="Tempo (anos)",ylab="S(t)",type="l",lty=1)
lines(c(0,t),c(1,st1),type="l",lty=4)
points(20,0.495,pch=16)
points(10,0.5,pch=16)
points(10,0.916,pch=16)
lines(c(0,20),c(0.5,0.5),lty=3,type="l")
lines(c(0,10),c(0.916,0.916),lty=3,type="l")
lines(c(10,10),c(0.0,0.916),lty=3,type="l")
lines(c(20,20),c(0,0.495),lty=3,type="l")
legend(25,0.90,lty=1,"Grupo 1",bty="n")
legend(25,0.83,lty=4,"Grupo 2",bty="n")Seção 1.5.3 - Figura 1.5
par(mfrow=c(1,2))
time<-c(0:30)
Ht<-(time/60)^1 # taxa de falha acumulada constante
Ht1<-(time/45)^1.8 # taxa de falha acumulada crescente
Ht2<-(time/120)^0.5 # taxa de falha acumulada decrescente
plot(time,time*0,pch=" ",ylim=range(c(0,0.5)),xlim=range(c(0,27)),xlab="Tempo",
ylab="Taxa de falha acumulada", bty="n", cex.lab=1.0, cex.axis=1.0)
lines(time[time<24],Ht[time<24], type="l", lty=1)
lines(time[time<25],Ht1[time<25], type="l", lty=2)
lines(time[time<23],Ht2[time<23], type="l", lty=4)
legend(0,0.50, expression(paste(Lambda[1](t))), lty=1, bty="n", cex=1.2)
legend(0,0.46, expression(paste(Lambda[2](t))), lty=2, bty="n", cex=1.2)
legend(0,0.42, expression(paste(Lambda[3](t))), lty=4, bty="n", cex=1.2)
time<-c(0:30)
ht<-rep((1/60),30) # taxa de falha constante
ht1<-(1/45)*(time/45)^(1.8-1) # taxa de falha crescente
ht2<-(1/120)*(time/120)^(0.5-1) # taxa de falha decrescente
plot(time,time*0,pch=" ",ylim=range(c(0,0.05)),xlim=range(c(0,27)),xlab="Tempo",
ylab="Taxa de falha", bty="n", cex.lab=1.0, cex.axis=1.0)
lines(time[time<22],ht[time<22], type="l", lty=1)
lines(time[time<22],ht1[time<22], type="l", lty=2)
lines(time[time<22],ht2[time<22], type="l", lty=4)
legend(16,0.050, expression(paste(lambda[1](t))), lty=1, bty="n", cex=1.2)
legend(16,0.046, expression(paste(lambda[2](t))), lty=2, bty="n", cex=1.2)
legend(16,0.042, expression(paste(lambda[3](t))), lty=4, bty="n", cex=1.2)Capítulo 2
Técnicas não paramétricas
Seção 2.3 - Estimador de Kaplan-Meier
library(survival)
tempos<-c(1,2,3,3,3,5,5,16,16,16,16,16,16,16,16,1,1,1,1,4,5,7,8,10,10,12,16,16,16)
cens<-c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,0)
grupos<-c(rep("Controle",15),rep("Esteroide",14))
dados<-data.frame(tempos,cens,grupos)
ekm<-survfit(Surv(tempos,cens)~grupos,dados) # Kaplan-Meier
summary(ekm)
par(mfrow=c(1,1))
plot(ekm,lty=c(2,1),xlab="Tempo (em semanas)",ylab="S(t) estimada", mark.time=F)
legend(1,0.3,lty=c(2,1),c("Controle","Esteroide"),lwd=1,bty="n") ### Figura 2.3 ###
dados_est<-subset(dados, grupos=="Esteroide")
ss<-survfit(coxph(Surv(tempos,cens)~1,method="breslow",dados_est)) # Nelson-Aalen
summary(ss)
dados_est1<-subset(dados, grupos=="Controle")
ss1<-survfit(coxph(Surv(tempos,cens)~1,method="breslow",dados_est1)) # Nelson-Aalen
summary(ss1)
par(mfrow=c(1,1))
plot(ekm, lty=c(1,1), xlab="Tempo (em semanas)",ylab="S(t) estimada",
mark.time=F, cex.lab=1.1, cex.axis=1.1)
lines(c(0,ss$time), c(1,ss$surv), type="s", lty=4)
lines(c(0,3,16),c(1,0.857,0.857), type="s", lty=4)
legend(13,0.87,c("Controle"), bty="n")
legend(13,0.45,c("Esteroide"), bty="n")
legend(1,0.3,lty=c(1,4),c("Kaplan-Meier","Nelson-Aalen"),lwd=1, bty="n", cex=1.1) ### Intervalos de Confiança ###
ekm<- survfit(Surv(tempos,cens)~grupos, conf.type="plain")
summary(ekm)
ekm<- survfit(Surv(tempos,cens)~grupos, conf.type="log-log")
summary(ekm)
ekm<- survfit(Surv(tempos,cens)~grupos, conf.type="log")
summary(ekm)
ekm<- survfit(Surv(tempos,cens)~grupos) # Obs: default é conf.type="log"
summary(ekm)Seção 2.4.1 - Estimador de Nelson-Aalen
tempos<-c(1,2,3,3,3,5,5,16,16,16,16,16,16,16,16,1,1,1,1,4,5,7,8,10,10,12,16,16,16)
cens<-c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,0)
grupos<-c(rep("Controle",15),rep("Esteroide",14))
dados<-data.frame(tempos,cens,grupos)
dados_est<-subset(dados, grupos=="Esteroide")
ss<-survfit(coxph(Surv(tempos,cens)~1,method="breslow",dados_est))
summary(ss)
racum<- -log(ss$surv)
cbind(ss$time,racum) ekm<-survfit(Surv(tempos,cens)~grupos,conf.type="log-log",dados); print(ekm)Seção 2.5.1 - Exemplo: reincidência de tumor sólido
library(survival)
tempos<-c(3,4,5.7,6.5,6.5,8.4,10,10,12,15)
cens<-c(1,0,0,1,1,0,1,0,1,1)
ekm<-survfit(Surv(tempos,cens)~1,conf.type="plain")
summary(ekm)
plot(ekm,conf.int=T,xlab="Tempo(em meses)",ylab="S(t) estimada",mark.time=T)
ekm<-survfit(Surv(tempos,cens)~1,conf.type="plain")
aux<-summary(ekm)
tm <-aux$time[1] + sum(diff(aux$time)*aux$surv[-length(aux$surv)])
tm
source("https://docs.ufpr.br/~giolo/asa/arquivos/var_tm.txt")
var.tmedio<-var_tm(ekm)
var.tmedioSeção 2.6 - Comparação de curvas de sobrevivência
library(survival)
tempos<- c(1,2,3,3,3,5,5,16,16,16,16,16,16,16,16,1,1,1,1,4,5,7,8,10,10,12,16,16,16)
cens<-c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,0)
grupos<-c(rep(1,15), rep(2,14))
survdiff(Surv(tempos,cens)~grupos, rho=0)Seção 2.6.1 - Exemplo: Eficácia da Imunização pela Málaria
tempos<-c(7,8,8,8,8,12,12,17,18,22,30,30,30,30,30,30,8,8,9,10,10,14,
15,15,18,19,21,22,22,23,25,8,8,8,8,8,8,9,10,10,10,11,17,19)
cens<-c(rep(1,10), rep(0,6), rep(1,15), rep(1,13))
grupo<-c(rep(1,16), rep(2,15), rep(3,13))
library(survival)
ekm<-survfit(Surv(tempos,cens)~grupo)
summary(ekm)
plot(ekm, lty=c(1,4,2), xlab="Tempo(em dias)",ylab="S(t) estimada")
legend(1,0.3,lty=c(1,4,2),c("Grupo 1","Grupo 2","Grupo 3"),lwd=1,bty="n",cex=0.8)
survdiff(Surv(tempos,cens)~grupo, rho=0)
dados<-data.frame(tempos,cens,grupo)
survdiff(Surv(tempos,cens)~grupo, data=subset(dados, grupo %in% c(1,2)), rho=0)
survdiff(Surv(tempos,cens)~grupo, data=subset(dados, grupo %in% c(2,3)), rho=0)
survdiff(Surv(tempos,cens)~grupo, data=subset(dados, grupo %in% c(1,3)), rho=0)Seção 2.6.2 - Outros testes não paramétricos
tempos<-c(1,2,3,3,3,5,5,16,16,16,16,16,16,16,16,1,1,1,1,4,5,7,8,10,10,12,16,16,16)
cens<-c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,0)
grupos<-c(rep("Controle",15),rep("Esteroide",14))
library(coin)
logrank_test(Surv(tempos,cens)~factor(grupos),type="logrank")
logrank_test(Surv(tempos,cens)~factor(grupos),type="Gehan-Breslow")
logrank_test(Surv(tempos,cens)~factor(grupos),type="Tarone-Ware") survdiff(Surv(tempos,cens)~grupos, rho=0)
survdiff(Surv(tempos,cens)~grupos, rho=1) Capítulo 3
Modelos probabilísticos
Seção 3.6.1 - Exemplo 1: pacientes com câncer de bexiga
library(survival)
tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
ajust1<-survreg(Surv(tempos,cens)~1, dist='exponential')
ajust1
alpha<-exp(ajust1$coefficients[1])
alpha
ajust2<-survreg(Surv(tempos,cens)~1, dist='weibull')
ajust2
alpha<-exp(ajust2$coefficients[1])
gama<-1/ajust2$scale
cbind(gama,alpha)
ajust3<-survreg(Surv(tempos,cens)~1, dist='lognorm')
ajust3
ekm<-survfit(Surv(tempos,cens)~1)
time<-ekm$time
st<-ekm$surv
ste<- exp(-time/20.41)
stw<- exp(-(time/21.34)^1.54)
stln<- pnorm((-log(time)+ 2.72)/0.76)
cbind(time,st,ste,stw,stln)
par(mfrow=c(1,3))
plot(st, ste, pch=16, ylim=c(0,1), xlim=c(0,1),
xlab = "S(t): Kaplan-Meier", ylab="S(t): exponencial",
cex.lab=1.5,cex.axis=1.3)
abline(a=0,b=1)
plot(st, stw, pch=16, ylim=c(0,1), xlim=c(0,1),
xlab = "S(t): Kaplan-Meier", ylab="S(t): Weibull",
cex.lab=1.5,cex.axis=1.3)
abline(a=0,b=1)
plot(st, stln, pch=16, ylim=c(0,1), xlim=c(0,1),
xlab = "S(t): Kaplan-Meier", ylab="S(t): log-normal",
cex.lab=1.5,cex.axis=1.3)
abline(a=0,b=1)
par(mfrow=c(1,3), mar=c(5,5,2,2))
invst<-qnorm(st)
plot(time, -log(st), pch=16, xlab="Tempo", ylab="-log(S(t))",
cex.lab=1.5,cex.axis=1.3)
plot(log(time), log(-log(st)), pch=16, xlab="log(Tempo)", ylab="log(-log(S(t)))",
cex.lab=1.5,cex.axis=1.3)
plot(log(time), invst, pch=16, xlab="log(Tempo)", ylab=expression(Phi^-1*(S(t))),
cex.lab=1.5,cex.axis=1.3)
ajust1$loglik[2]
ajust2$loglik[2]
ajust3$loglik[2]
par(mfrow=c(1,2))
plot(ekm, conf.int=F, xlab="Tempo (em meses)", ylab="S(t)")
lines(c(0,time), c(1,stw), lty=2)
legend(25,0.8,lty=c(1,2),c("Kaplan-Meier","Weibull"),bty="n",cex=0.9)
plot(ekm, conf.int=F, xlab="Tempo (em meses)", ylab="S(t)")
lines(c(0,time), c(1,stln), lty=2)
legend(25,0.8,lty=c(1,2),c("Kaplan-Meier","Log-normal"), bty="n",cex=0.9)Ajuste do modelo gama generalizada utilizando o pacote "flexsurv"
install.packages("flexsurv") # instalação do pacote
library(flexsurv)
tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
ajust4<-flexsurvreg(Surv(tempos,cens)~1, dist='gengamma')
ajust4$loglik
ajust4
summary(ajust4)Obtenção de var(tm) e var (tmed) via reamostragem bootstrap. Clique ➔ Função (Giolo, 2010)
Seção 3.6.2 - Exemplo 2: pacientes em quimioterapia
library(survival)
tempos<-c(7,8,10,12,13,14,19,23,25,26,27,31,31,49,59,64,87,89,107,117,
119,130,148,153,156,159,191,222,200,203,210,220,220,228,230,
233,235,240,240,240,241,245,247,248,250)
cens<-c(1,1,1,1,1,0,1,1,0,1,1,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0)
ekm<-survfit(Surv(tempos,cens)~1,type=c("kaplan-meier"))
summary(ekm)
ajust1<-survreg(Surv(tempos,cens)~1,dist='exponential')
ajust1
alpha<-exp(ajust1$coefficients[1])
alpha
ajust2<-survreg(Surv(tempos,cens)~1,dist='weibull')
ajust2
alpha<-exp(ajust2$coefficients[1])
gama<-1/ajust2$scale
cbind(gama, alpha)
ajust3<-survreg(Surv(tempos,cens)~1,dist='lognorm')
ajust3
library(flexsurv) ## instale previamente o pacote
ajust4<-flexsurvreg(Surv(tempos,cens)~1, dist='gengamma')
ajust4
ajust1$loglik[1]
ajust2$loglik[1]
ajust3$loglik[1]
ajust4$loglik[1]
time<-ekm$time
st<-ekm$surv
ste<- exp(-time/262)
stw<- exp(-(time/284.6594)^0.8347)
stln<- pnorm((-log(time)+ 5.181)/1.724)
cbind(time,st,ste,stw,stln) ### Método 1 ###
par(mfrow=c(1,3))
plot(ste,st,pch=16,ylim=range(c(0.3,1)), xlim=range(c(0.3,1)),
ylab = "S(t): Kaplan-Meier", xlab="S(t): exponencial")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(stw,st,pch=16,ylim=range(c(0.3,1)), xlim=range(c(0.3,1)),
ylab = "S(t): Kaplan-Meier", xlab="S(t): Weibull")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(stln,st,pch=16,ylim=range(c(0.3,1)), xlim=range(c(0.3,1)),
ylab = "S(t): Kaplan-Meier", xlab="S(t): log-normal")
lines(c(0,1), c(0,1), type="l", lty=1)
### Método 2 ###
par(mar=c(5,5,1.5,1))
par(mfrow=c(1,3))
invst<-qnorm(st)
plot(time, -log(st),pch=16,xlab="Tempos",ylab="-log[S(t)]", cex.lab=1.5, cex.axis=1.3)
plot(log(time),log(-log(st)),pch=16,xlab="log(Tempos)",ylab="log[-log(S(t))]",
cex.lab=1.5, cex.axis=1.3)
plot(log(time),invst, pch=16,xlab="log(Tempos)",
ylab=expression(Phi^-1*(S(t))), cex.lab=1.5, cex.axis=1.3) ### Figura 3.10 ###
par(mfrow=c(1,1))
plot(ekm, conf.int=F, xlab="Tempos (em dias)", ylab="S(t)", cex.lab=1.2, cex.axis=1.1)
lines(c(0,time),c(1,stln), lty=2)
legend(25,0.4, lty=c(1,2), c("Kaplan-Meier","Log-normal"), bty="n", cex=1.0)Capítulo 4
Modelos de regressão paramétricos
Seção 4.6.1 - Sobrevida de pacientes com leucemia aguda
rm(list=ls(all=TRUE))
temp<-c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
cens<-rep(1,17)
wbc<-c(2300,750,4300,2600,6000,10500,10000,17000,5400,7000,9400,32000,35000,
100000,100000,52000,100000)
lwbc<-c(3.36,2.88,3.63,3.41,3.78,4.02,4.00,4.23,3.73,3.85,3.97,4.51,4.54,
5.00,5.00,4.72,5.00)
dados<-cbind(temp,cens,wbc,lwbc)
library(survival)
dados<-data.frame(dados)
i<-order(dados$temp)
dados<-dados[i,]
ekm<- survfit(Surv(dados$temp,dados$cens)~1)
summary(ekm)
st<-ekm$surv
temp<-ekm$time
invst<-qnorm(st)
### Figura 4.1 ###
par(mfrow=c(1,3))
plot(temp, -log(st), pch=16, xlab="Tempos", ylab="-log(S(t))")
plot(log(temp), log(-log(st)), pch=16, xlab="log(Tempos)", ylab="log(-log(S(t))")
plot(log(temp), invst, pch=16, xlab="log(Tempos)", ylab=expression(Phi^-1*(S(t)))) ### Figura 4.2 ###
ajust1<-survreg(Surv(dados$temp, dados$cens)~1, dist='exponential')
t<-dados$temp
bo<- ajust1$coefficients[1]
res<- t*exp(-bo)
resm<-dados$cens-res
par(mar=c(5,5,2,2)) # c(bottom,left,top,right)
par(mfrow=c(1,2))
plot(dados$wbc,resm, xlab = "WBC", ylab="Resíduos martingal", pch=20)
lines(smooth.spline(dados$wbc,resm),col=1)
title("(a)", cex.main=1.0)
plot(dados$lwbc,resm, xlab = "log(WBC)", ylab="Resíduos martingal", pch=20)
lines(smooth.spline(dados$lwbc,resm),col=1)
title("(b)", cex.main=1.0) ### Ajuste dos modelos ###
ajust1<-survreg(Surv(dados$temp, dados$cens)~dados$lwbc, dist='exponential')
ajust1
ajust1$loglik
ajust2<-survreg(Surv(dados$temp, dados$cens)~dados$lwbc, dist='weibull')
ajust2
ajust2$loglik
gama<-1/ajust2$scale; gama ### Figura 4.3 ###
ajust1<-survreg(Surv(dados$temp, dados$cens)~dados$lwbc, dist='exponential')
x<-dados$lwbc
t<-dados$temp
bo<- ajust1$coefficients[1]
b1<- ajust1$coefficients[2]
res<- t*exp(-bo-b1*x)
ekm <- survfit(Surv(res,dados$cens)~1, type=c("kaplan-meier"))
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, conf.int=F, lty=c(1,1), xlab="Resíduos de Cox-Snell",
ylab="S(ei) estimada", cex.lab=1.2, cex.axis=1.2)
res<-sort(res)
exp1<-exp(-res)
lines(res,exp1,lty=2)
legend(1.1,0.9, lty=c(1,2), c("Kaplan-Meier","Exponencial padrão"),
lwd=1, bty="n", cex=1.0)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, ylim=c(0,1), xlim=c(0,1), xlab="S(ei): Kaplan-Meier",
ylab= "S(ei): exponencial padrão", pch=16, cex.lab=1.2, cex.axis=1.2)
abline(a=0,b=1, col="gray", lty=2)
points(st, sexp1, pch=16) ### Figura 4.4 ###
ajust1<-survreg(Surv(temp, cens)~lwbc, dist='exponential', data=dados)
ajust1
x1<-4.0; temp1<-0:150
ax1<-exp(ajust1$coefficients[1]+ajust1$coefficients[2]*x1)
ste1<-exp(-(temp1/ax1))
x1<-3.0; temp2<-0:150
ax2<-exp(ajust1$coefficients[1]+ajust1$coefficients[2]*x1)
ste2<-exp(-(temp2/ax2))
par(mfrow=c(1,1))
plot(temp1,temp1*0,pch=" ",ylim=range(c(0,1)), xlim=range(c(0,150)),
xlab="Tempos (em semanas)",ylab="S(t) estimada",bty="n")
lines(temp1,ste1,lty=2)
lines(temp2,ste2,lty=5)
lines(c(100,100),c(0.05,0.70), col="gray")
legend(5,0.3,lty=c(2,5),c("log(wbc) = 4","log(wbc) = 3"),lwd=1, bty="n")
points(100,0.172, pch=20)
points(100,0.559, pch=20)Seção 4.6.2 - Grupos de pacientes com leucemia aguda
temp<-c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65,56,65,17,7,16,22,
3,4,2,3,8,4,3,30,4,43)
cens<-c(rep(1,17),rep(1,16))
lwbc<-c(3.36,2.88,3.63,3.41,3.78,4.02,4,4.23,3.73,3.85,3.97,4.51,4.54,5,5,4.72,
5,3.64,3.48,3.6,3.18,3.95,3.72,4,4.28,4.43,4.45,4.49,4.41,4.32,4.90,5,5)
grupo<-c(rep(0,17),rep(1,16))
library(survival)
ekm1<-survfit(Surv(temp,cens)~grupo)
summary(ekm1)
st1<-ekm1[1]$surv
time1<-ekm1[1]$time
invst1<-qnorm(st1)
st2<-ekm1[2]$surv
time2<-ekm1[2]$time
invst2<-qnorm(st2)
### Figura 4.5 ###
par(mfrow=c(1,3))
plot(time1, -log(st1), pch=16, xlab="Tempo", ylab="-log[S(t)]",
cex.lab=1.5, cex.axis=1.3)
points(time2, -log(st2))
legend(100, 0.7, pch=c(16,1), c("Ag+","Ag -"), bty="n", cex=1.2)
plot(log(time1), log(-log(st1)), pch=16, xlab="log(Tempo)",
ylab="log[-log(S(t))]", cex.lab=1.5, cex.axis=1.3)
points(log(time2), log(-log(st2)))
legend(3.5, -1.5, pch=c(16,1), c("Ag+","Ag -"), bty="n", cex=1.2)
plot(log(time1), invst1,pch=16, xlab="log(Tempo)", ylab=expression(Phi^-1*(S(t))),
cex.lab=1.5, cex.axis=1.3)
points(log(time2), invst2)
legend(0.5,-1, pch=c(16,1), c("Ag+","Ag -"), bty="n", cex=1.2) dados<-data.frame(cbind(temp,cens,lwbc,grupo))
attach(dados)
library(survival)
ajust1<-survreg(Surv(temp,cens)~1, dist='exponential')
ajust1
ajust2<-survreg(Surv(temp,cens)~lwbc, dist='exponential')
ajust2
ajust3<-survreg(Surv(temp,cens)~grupo, dist='exponential')
ajust3
ajust4<-survreg(Surv(temp,cens)~lwbc + grupo, dist='exponential')
ajust4
ajust5<-survreg(Surv(temp,cens)~lwbc + grupo + lwbc*grupo, dist='exponential')
ajust5
summary(ajust4) ### Figura 4.6 ###
t<-temp
x1<-lwbc
x2<-grupo
bo<-6.83
b1<- -0.7
b2<- -1.02
res<- t*exp(-bo-b1*x1-b2*x2)
ekm <- survfit(Surv(res,dados$cens)~1, type=c("kaplan-meier"))
par(mfrow=c(1,2))
plot(ekm, conf.int=F, lty=c(1,1), xlab="Resíduos de Cox-Snell",
ylab="S(ei) estimada", cex.lab=1.4, cex.axis=1.3)
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=2)
legend(1.7, 0.9, lty=c(1,2), c("Kaplan-Meier","Exponencial padrão"),
lwd=1, bty="n", cex=1.2)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(ei): Kaplan-Meier", xlim=c(0,1), ylim=c(0,1),
ylab= "S(ei: exponencial padrão", pch=16, cex.lab=1.4, cex.axis=1.3)
abline(a=0,b=1, col="gray", lty=2)
points(st, sexp1, pch=16) ### Figura 4.7 ###
x1<-4.0; x2<-0.0
temp1<-0:150
ax1<-exp(6.83-0.70*x1-1.02*x2)
ste1<-exp(-(temp1/ax1))
x1<-3.0; x2<-0.0
temp2<-0:150
ax2<-exp(6.83-0.70*x1-1.02*x2)
ste2<-exp(-(temp2/ax2))
par(mfrow=c(1,2))
plot(temp1, temp1*0, pch=" ", ylim=range(c(0,1.09)), xlim=range(c(0,150)),
xlab="Tempos (em semanas)", ylab="S(t) estimada", bty="n",
cex.lab=1.2, cex.axis=1.2)
lines(temp1, ste1, lty=1)
lines(temp2, ste2, lty=2)
legend(75, 0.8, lty=c(1,2), c("log(wbc) = 4","log(wbc) = 3"),
lwd=1, bty="n", cex=1.0)
legend(50, 1.16, c("Ag+"), bty="n", cex=1.15)
x1<-4.0; x2<-1.0
temp1<-0:150
ax1<-exp(6.83-0.70*x1-1.02*x2)
ste1<-exp(-(temp1/ax1))
x1<-3.0; x2<-1.0
temp2<-0:150
ax2<-exp(6.83-0.70*x1-1.02*x2)
ste2<-exp(-(temp2/ax2))
plot(temp1, temp1*0, pch=" ", ylim=range(c(0,1.09)), xlim=range(c(0,150)),
xlab="Tempos (em semanas)", ylab="S(t) estimada",
bty="n", cex.lab=1.2, cex.axis=1.2)
lines(temp1, ste1, lty=1)
lines(temp2, ste2, lty=2)
legend(75, 0.8, lty=c(1,2), c("log(wbc) = 4","log(wbc) = 3"),
lwd=1, bty="n", cex=1.0)
legend(50, 1.16, c("Ag-"), bty="n", cex=1.15) ### Figura 4.8 ###
x1<-4.0; x2<-0.0
temp1<-0:150
risco1<-1/(exp(6.83-0.70*x1-1.02*x2))
risco1<-rep(risco1, 151)
x1<-3.0; x2<-0.0
temp2<-0:150
risco2<-1/(exp(6.83-0.70*x1-1.02*x2))
risco2<-rep(risco2, 151)
plot(temp1,temp1*0, pch=" ", ylim=range(c(0,0.11)), xlim=range(c(0,150)),
xlab="Tempos (em semanas)",ylab="Taxa de falha estimada", bty="n",
cex.lab=1.2, cex.axis=1.2)
lines(temp1, risco1, lty=1)
lines(temp2, risco2, lty=2)
legend(75, 0.09, lty=c(1,2), c("log(wbc) = 4","log(wbc) = 3"),
lwd=1, bty="n", cex=1.0)
legend(50, 0.118, c("Ag+"), bty="n", cex=1.15)
x1<-4.0; x2<-1.0
temp1<-0:150
risco1<-1/(exp(6.83-0.70*x1-1.02*x2))
risco1<-rep(risco1,151)
x1<-3.0; x2<-1.0
temp2<-0:150
risco2<-1/(exp(6.83-0.70*x1-1.02*x2))
risco2<-rep(risco2, 151)
plot(temp1, temp1*0, pch=" ", ylim=range(c(0,0.11)), xlim=range(c(0,150)),
xlab="Tempos (em semanas)", ylab="Taxa de falha estimada", bty="n",
cex.lab=1.2, cex.axis=1.2)
lines(temp1, risco1, lty=1)
lines(temp2, risco2, lty=2)
legend(75, 0.09, lty=c(1,2), c("log(wbc) = 4","log(wbc) = 3"),
lwd=1, bty="n", cex=1.0)
legend(50, 0.118, c("Ag-"), bty="n", cex=1.15)Seção 4.6.3 - Estudo sobre aleitamento materno
rm(list=ls(all=TRUE))
desmame<-read.table("https://docs.ufpr.br/~giolo/asa/dados/desmame.txt",h=T)
attach(desmame)
library(survival)
ekm<- survfit(Surv(tempo,cens)~V4)
summary(ekm)
survdiff(Surv(tempo,cens)~V4, rho=0)
### Figura 4.9 ###
par(mfrow=c(1,1))
plot(ekm, lty=c(1,4), mark.time=F, xlab="Tempo até o desmame (em meses)", ylab="S(t)")
text(18.5, 0.93, c("Dificuldades para amamentar"), bty="n", cex=0.9)
legend(15.5, 0.9, lty=c(4), c("Sim"), bty="n", cex=0.9)
legend(18.5, 0.9, lty=c(1), c("Não"), bty="n", cex=0.9) ### Tabela 4.9 ###
library(flexsurv)
fitg <- flexsurvreg(Surv(tempo, cens) ~ 1, data = desmame, dist="gengamma"); fitg
fitg <- flexsurvreg(Surv(tempo, cens) ~ V1, data = desmame, dist="gengamma"); fitg
# similar para os demais # ajust1<-survreg(Surv(tempo,cens)~V1+V3+V4+V6, dist='lognorm')
ajust1
summary(ajust1)
### Figura 4.10 ###
xb<-ajust1$coefficients[1] + ajust1$coefficients[2]*V1 + ajust1$coefficients[3]*V3 +
ajust1$coefficients[4]*V4 + ajust1$coefficients[5]*V6
sigma<-ajust1$scale
res<-(log(tempo)-(xb))/sigma # resíduos padronizados
resid<-exp(res) # exponencial dos resíduos padronizados
ekm<- survfit(Surv(resid,cens)~1)
resid<-ekm$time
sln<-pnorm(-log(resid))
par(mfrow=c(1,2))
plot(ekm$surv,sln, xlab="S(ei*): Kaplan-Meier",
ylab="S(ei*): log-normal padrão", pch=16)
plot(ekm, conf.int=F, mark.time=F, xlab="Resíduos (ei*)",
ylab="Sobrevivência estimada", pch=16)
lines(resid, sln, lty=2)
legend(1.3,0.8,lty=c(1,2),c("Kaplan-Meier","Log-normal padrão"),cex=0.8,bty="n")
### Figura 4.11 ###
ei<- -log(1-pnorm(res)) # resíduos de Cox-Snell
ekm1<-survfit(Surv(ei,cens)~1)
t<-ekm1$time
st<-ekm1$surv
sexp<-exp(-t)
par(mfrow=c(1,2))
plot(st,sexp,xlab="S(ei): Kaplan-Meier",ylab="S(ei): exponencial padrão",pch=16)
plot(ekm1,conf.int=F, mark.time=F, xlab="Resíduos de Cox-Snell",
ylab="Sobrevivência estimada")
lines(t, sexp, lty=4)
legend(1.0,0.8,lty=c(1,4),c("Kaplan-Meier","Exponencial padrão"),cex=0.8,bty="n")Capítulo 5
Modelo de regressão de Cox
Seção 5.7.1 - Estudo sobre câncer de laringe
laringe<-read.table("https://docs.ufpr.br/~giolo/asa/dados/laringe.txt", h=T)
attach(laringe)
library(survival)
laringe$est<-factor(estadio)
fit2<-coxph(Surv(tempos,cens)~est, data=laringe, method="breslow")
summary(fit2); fit2$loglik
fit3<- coxph(Surv(tempos,cens)~est+idade, data=laringe, method="breslow")
summary(fit3); fit3$loglik
fit4<-coxph(Surv(tempos,cens)~est+idade+est*idade, data=laringe, method="breslow")
summary(fit4); fit4$loglik resid(fit4,type="scaledsch")
cox.zph(fit4, transform="identity", terms=TRUE, global=TRUE) # identity: g(t)=t
### Figura 5.4 ###
par(mfrow=c(2,4))
plot(cox.zph(fit4,transform="identity",terms=FALSE)) Ht<-basehaz(fit4, centered=F)
tempos<-Ht$time
H0<-Ht$hazard
S0<- exp(-H0)
round(cbind(tempos,S0,H0), digits=3) ### Figura 5.5 ###
resm<-resid(fit4,type="martingale")
cens<-laringe$cens
res<-cens - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, cens)~1)
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, mark.time=F, conf.int=F, xlab="Resíduos de Cox-Snell", ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(0.7, 0.9, lty=c(1,3), c("Kaplan Meier","Exponencial padrão"),
lwd=1, bty="n", cex=0.9)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab="S(e): exponencial padrão", pch=16) ### Figura 5.6 ###
tt<-sort(tempos)
aux1<-as.matrix(tt)
n<-nrow(aux1)
aux2<-as.matrix(cbind(tempos,S0))
S00<-rep(max(aux2[,2]),n)
for(i in 1:n){
if(tt[i]> min(aux2[,1])){
i1<- aux2[,1]<= tt[i]
S00[i]<-min(aux2[i1,2])}}
ts0<-cbind(tt,S00)
ts0
b<-fit4$coefficients
id<-50
st1<-S00^(exp(b[4]*id)) # S(t|x) estádio I e idade = 50 anos
st2<-S00^(exp(b[1]+((b[4]+b[5])*id))) # S(t|x) estádio II e idade = 50 anos
st3<-S00^(exp(b[2]+((b[4]+b[6])*id))) # S(t|x) estádio III e idade = 50 anos
st4<-S00^(exp(b[3]+((b[4]+b[7])*id))) # S(t|x) estádio IV e idade = 50 anos
id<-65
st11<-S00^(exp(b[4]*id)) # S(t|x) estádio I e idade = 65 anos
st21<-S00^(exp(b[1]+((b[4]+b[5])*id))) # S(t|x) estádio II e idade = 65 anos
st31<-S00^(exp(b[2]+((b[4]+b[6])*id))) # S(t|x) estádio III e idade = 65 anos
st41<-S00^(exp(b[3]+((b[4]+b[7])*id))) # S(t|x) estádio IV e idade = 65 anos
par(mfrow=c(1,2))
plot(tt,st1, type="s", ylim=range(c(0,1)), xlab="Tempo (em meses)",
ylab="S(t|x)", lty=1, cex.lab=0.9)
lines(tt,st2, type="s", lty=2)
lines(tt,st3, type="s", lty=3)
lines(tt,st4, type="s", lty=4)
legend(0, 0.25, lty=c(1,2,3,4), c("Estádio I","Estádio II","Estádio III","Estádio IV"),
lwd=1, bty="n", cex=0.8)
title("Idade = 50 anos", cex.main=0.9)
plot(tt,st11, type="s", ylim=range(c(0,1)), xlab="Tempo (em meses)",
ylab="S(t|x)", lty=1, cex.lab=0.9)
lines(tt,st21,type="s",lty=2)
lines(tt,st31,type="s",lty=3)
lines(tt,st41,type="s",lty=4)
legend(0, 0.25, lty=c(1,2,3,4), c("Estádio I","Estádio II","Estádio III","Estádio IV"),
lwd=1, bty="n", cex=0.8)
title("Idade = 65 anos", cex.main=0.9) ### Figura 5.7 ###
Ht1<- -log(st1); Ht2<- -log(st2); Ht3<- -log(st3); Ht4<- -log(st4)
Ht11<- -log(st11); Ht21<- -log(st21); Ht31<- -log(st31); Ht41<- -log(st41)
par(mfrow=c(1,2))
plot(tt, Ht1, type="s", ylim=range(c(0,4)), xlim=c(0,8), xlab="Tempo (em meses)",
ylab="Risco acumulado", lty=1, cex.lab=0.9)
lines(tt, Ht2, type="s", lty=2)
lines(tt, Ht3, type="s", lty=3)
lines(tt, Ht4, type="s", lty=4)
legend(0,3.5, lty=c(1,2,3,4), c("Estádio I","Estádio II","Estádio III","Estádio IV"),
lwd=1, bty="n", cex=0.8)
title("Idade = 50 anos", cex.main=0.9)
plot(tt, Ht11, type="s", ylim=range(c(0,4)), xlim=c(0,8), xlab="Tempo (em meses)",
ylab="Risco acumulado", lty=1, cex.lab=0.9)
lines(tt, Ht21, type="s", lty=2)
lines(tt, Ht31, type="s", lty=3)
lines(tt, Ht41, type="s", lty=4)
legend(0,3.5, lty=c(1,2,3,4), c("Estádio I","Estádio II","Estádio III","Estádio IV"),
lwd=1, bty="n", cex=0.8)
title("Idade = 65 anos", cex.main=0.9)Seção 5.7.2 - Análise dos dados de aleitamento materno
rm(list=ls(all=TRUE))
library(survival)
desmame<-read.table("https://docs.ufpr.br/~giolo/asa/dados/desmame.txt",h=T)
attach(desmame)
fit<-coxph(Surv(tempo,cens)~V1+V3+V4+V6, data=desmame, method="breslow")
summary(fit)
fit$loglik### Figura 5.8 ###
par(mar=c(5,5,2,2)) # c(bottom,left,top,right)
par(mfrow=c(2,2))
fit1<-coxph(Surv(tempo[V1==0],cens[V1==0])~1, data=desmame, method="breslow")
summary(fit1)
fit1$loglik
ss<- survfit(fit1)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
plot(ss$time[ss$time<18],log(H0)[ss$time <18], ylim=c(-4,0.3), xlim=range(c(0,20)),
xlab="Tempo",ylab = expression(log(Lambda[0]* (t))),bty="n",type="s")
fit2<-coxph(Surv(tempo[V1==1],cens[V1==1]) ~ 1, data=desmame, method="breslow")
ss<- survfit(fit2)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
lines(ss$time[ss$time <=20],log(H0)[ss$time <=20],type="s",lty=2)
legend(7,-2.5,lty=c(2,1),c("V1 = 1 (não)","V1 = 0 (sim)"),lwd=1,bty="n",cex=0.9)
title("V1: Experiência de amamentação", cex.main=1.0)
fit1<-coxph(Surv(tempo[V3==0],cens[V3==0]) ~ 1, data=desmame, method="breslow")
summary(fit1)
fit1$loglik
ss<- survfit(fit1)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
plot(ss$time[ss$time<=18],log(H0)[ss$time<=18], ylim=c(-4,0.3), xlim=c(0,20),
xlab="Tempo",ylab = expression(log(Lambda[0] * (t))),bty="n",type="s")
fit2<-coxph(Surv(tempo[V3==1],cens[V3==1]) ~ 1, data=desmame, method="breslow")
ss<- survfit(fit2)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
lines(ss$time[ss$time<=17],log(H0)[ss$time<=17],type="s",lty=2)
legend(6.5,-2.5,lty=c(2,1),c("V3 = 1 (< 6 meses)","V3 = 0 (> 6 meses)"),lwd=1,bty="n",cex=0.9)
title("V3: Conceito de amamentação", cex.main=1.0)
fit1<-coxph(Surv(tempo[V4==0],cens[V4==0]) ~ 1, data=desmame, method="breslow")
summary(fit1)
fit1$loglik
ss<- survfit(fit1)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
plot(ss$time[ss$time<=18],log(H0)[ss$time<=18], ylim=c(-4,0.3), xlim=c(0,20),
xlab="Tempo", ylab = expression(log(Lambda[0]* (t))),bty="n",type="s")
fit2<-coxph(Surv(tempo[V4==1],cens[V4==1]) ~ 1, data=desmame, method="breslow")
ss<- survfit(fit2)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
lines(ss$time[ss$time<=17],log(H0)[ss$time<=17],type="s",lty=2)
legend(7,-2.5,lty=c(2,1),c("V4 = 1 (sim)","V4 = 0 (não)"),lwd=1,bty="n",cex=0.9)
title("V4: Dificuldades de amamentar", cex.main=1.0)
fit1<-coxph(Surv(tempo[V6==0],cens[V6==0]) ~ 1 , data=desmame, method="breslow")
summary(fit1)
fit1$loglik
ss<- survfit(fit1)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
plot(ss$time[ss$time<=18],log(H0)[ss$time<=18], ylim=c(-4,0.3), xlim=c(0,20),
xlab="Tempo",ylab = expression(log(Lambda[0]* (t))),bty="n",type="s")
fit2<-coxph(Surv(tempo[V6==1],cens[V6==1]) ~ 1, data=desmame, method="breslow")
ss<- survfit(fit2)
s0<-round(ss$surv,digits=5)
H0<- -log(s0)
lines(ss$time[ss$time<=17],log(H0)[ss$time<=17],type="s",lty=2)
legend(6.5,-2.5,lty=c(2,1),c("V6 = 1 (não)","V6 = 0 (sim)"),lwd=1,bty="n",cex=0.9)
title("V6: Leite materno na maternidade", cex.main=1.0) ### Figura 5.9 ###
resid(fit, type="scaledsch")
cox.zph(fit, transform="identity") # g(t) = t
par(mfrow=c(2,2))
plot(cox.zph(fit))
### Correlação: rho ###
i<-order(desmame$tempo)
desmame1<-desmame[i,]
falhas<-desmame1[desmame1$cens==1,]
rps<-resid(fit, type="scaledsch")
cor(falhas[,2],rps[,1])
cor(falhas[,2],rps[,2])
cor(falhas[,2],rps[,3])
cor(falhas[,2],rps[,4]) ### Figura 5.10 ###
resm<-resid(fit,type="martingale")
res<-cens - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, cens)~1)
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, mark.time=F, conf.int=F, xlab="Resíduos de Cox-Snell", ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(1,0.8, lty=c(1,3), c("Kaplan Meier","Exponencial padrão"), lwd=1, bty="n", cex=1.0)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st,sexp1, ylim= c(0,1), xlab="S(e): Kaplan-Meier", ylab="S(e): exponencial padrão", pch=16)
abline(a=0,b=1, lty=3)Seção 5.7.3 - Análise dos dados de leucemia pediátrica
rm(list=ls())
leuc<-read.table("https://docs.ufpr.br/~giolo/asa/dados/leucemia.txt", h=T)
attach(leuc)
idadec<-ifelse(idade>96,1,0)
leuinic<-ifelse(leuini>75,1,0)
zpesoc<-ifelse(zpeso>-2,1,0)
zestc<-ifelse(zest>-2,1,0)
vacc<-ifelse(vac>15,1,0)
pasc<-ifelse(pas>5,1,0)
riskc<-ifelse(risk>1.7,1,0)
r6c<-r6
leucc<-as.data.frame(cbind(leuinic,tempos,cens,idadec,zpesoc,zestc,pasc,vacc,riskc,r6c))
detach(leuc); attach(leucc)
library(survival)
fit<-coxph(Surv(tempos,cens)~leuinic + idadec + zpesoc + zestc + pasc + vacc +
riskc + r6c, data=leucc, method="breslow")
summary(fit); -2*fit$loglik[2]
fit3<-coxph(Surv(tempos,cens)~leuinic + idadec + zpesoc + pasc + vacc,
data=leucc, method="breslow")
summary(fit3)
-2*fit3$loglik[2] ### Figura 5.11 ###
par(mfrow=c(2,3))
fit<-coxph(Surv(tempos[leuinic==1],cens[leuinic==1]) ~ 1, data=leucc, method="breslow")
ss<- survfit(fit)
s0<-round(ss$surv, digits=5)
H0<- -log(s0)
plot(ss$time, log(H0), xlab="Tempos", ylim=range(c(-5,1)),
ylab=expression(log(Lambda[0]*(t))), bty="n", type="s")
fit<-coxph(Surv(tempos[leuinic==0],cens[leuinic==0]) ~ 1, data=leucc, method="breslow")
ss<- survfit(fit)
s0<-round(ss$surv, digits=5)
H0<- -log(s0)
lines(ss$time, log(H0), type="s", lty=4)
legend(1.5, -4, lty=c(4,1), c("leuini < 75","leuini > 75"), lwd=1, bty="n", cex=0.9)
title("LEUINI")
# Obs: Análogo para as demais covariáveis. ### Figura 5.12 ###
resid(fit3, type="scaledsch")
cox.zph(fit3, transform="identity") # g(t) = t
par(mfrow=c(2,3))
plot(cox.zph(fit3)) ### Figura 5.13 ###
par(mfrow=c(1,2))
rd<-resid(fit3, type="deviance") # resíduos deviance
rm<-resid(fit3, type="martingale") # resíduos martingal
pl<-fit3$linear.predictors
plot(pl, rm, xlab="Preditor linear", ylab="Resíduos martingal", pch=16)
plot(pl, rd, xlab="Preditor linear", ylab="Resíduos deviance" , pch=16) ### Figura 5.14 ###
par(mfrow=c(2,3))
dfbetas<-resid(fit3, type="dfbeta")
plot(leuinic,dfbetas[,1], xlab="Leuini", ylab="Influência para leuini")
plot(idadec, dfbetas[,2], xlab="Idade", ylab="Influência para idade")
plot(zpesoc, dfbetas[,3], xlab="Zpeso", ylab="Influência para zpeso")
plot(pasc, dfbetas[,4], xlab="Pas", ylab="Influência para pas")
plot(vacc, dfbetas[,5], xlab="Vac", ylab="Influência para vac") ### Figura 5.15 ###
resm<-resid(fit3,type="martingale")
res<-leucc$cens - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, leucc$cens)~1)
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, mark.time=F, conf.int=F, xlab="resíduos", ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res,exp1, lty=3)
legend(1,0.8, lty=c(1,3), c("Kaplan Meier","Exponencial padrão"), lwd=1, bty="n", cex=0.9)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st,sexp1, xlab="S(e): Kaplan-Meier", ylab="S(e): exponencial padrão", pch=16)
abline(a=0,b=1,lty=3)Capítulo 6
Extensões do modelo de Cox
6.4 - Análise dos dados de pacientes HIV
rm(list=ls())
aids<-read.table("https://docs.ufpr.br/~giolo/asa/dados/aids.txt", h=T)
aids1<-subset(aids, ti<tf)
aids2<-na.omit(aids1[,1:7])
attach(aids2); aids2
library(survival)
fit1<- coxph(Surv(ti,tf,cens)~id + factor(grp), method="breslow", data=aids2)
summary(fit1) ### Figura 6.1 ###
resm<-resid(fit1, type="martingale")
res<-aids2$cens - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, aids2$cens)~1)
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, mark.time=F, conf.int=F, xlab="Resíduos de Cox-Snell",
ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(0.6,1.0,lty=c(1,3),c("Kaplan-Meier","Exponencial padrão"),
lwd=1, bty="n", cex=0.9)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st,sexp1,xlab="S(e): Kaplan-Meier",ylab="S(e): exponencial padrão",pch=16)
abline(a=0, b=1, lty=3) ### Figura 6.2 ###
par(mar=c(5,4.5,2,2))
par(mfrow=c(1,1))
Ht<-basehaz(fit1, centered=F)
tempos<-Ht$time
H0<-Ht$hazard
H01<-H0*exp(fit1$coef[1]*30) # soronegativo 30 anos
H02<-H0*exp(fit1$coef[1]*30+fit1$coef[3]) # ARC 30 anos
H03<-H0*exp(fit1$coef[1]*30+fit1$coef[4]) # AIDS 30 anos
H04<-H0*exp(fit1$coef[1]*30+fit1$coef[2]) # assintomático 30 anos
plot(tempos, H01, type="s", lty=4,col=1, ylim=c(0,2), xlab="Tempo (em dias)",
ylab=expression(hat(Lambda)(t)))
points(tempos[tempos<=130],H04[tempos<=130], type="s",lty=1,col=1)
points(tempos[tempos>=130&tempos<=247],H02[tempos>=130&tempos<=247],
type="s",lty=1,col=1)
points(tempos[tempos>=247],H03[tempos>=247], type="s", lty=1, col=1)
lines(c(130,130), c(0.01,0.24), lty=3, col=1)
lines(c(247,247), c(0.49,0.72), lty=3, col=1)
axis(1, c(130, 247), c(expression(t[1]), expression(t[2])))
legend(0,1.8, lty=c(4,1), c("30 anos e grupo HIV soronegativo",
"30 anos com mudança de grupo"), bty="n")Seção 6.5 - Análise dos dados de leucemia pediátrica
rm(list=ls())
library(survival)
leucc<-read.table("https://docs.ufpr.br/~giolo/asa/dados/leucc.txt",h=T)
attach(leucc)
fit1<-coxph(Surv(tempos,cens)~idadec+zpesoc+pasc+vacc+strata(leuinic),
data=leucc, method="breslow")
summary(fit1) leucc1<-subset(leucc,leuinic==0); leucc2<-subset(leucc,leuinic==1)
fit2<-coxph(Surv(tempos,cens)~idadec+zpesoc+pasc+vacc,leucc1,method="breslow")
fit3<-coxph(Surv(tempos,cens)~idadec+zpesoc+pasc+vacc,leucc2,method="breslow")
summary(fit2); summary(fit3)
trv<-2*(-fit1$loglik[2]+fit2$loglik[2]+fit3$loglik[2])
trv; 1-pchisq(trv,4)
cox.zph(fit1,transform="identity")
### Figura 6.3 ###
par(mfrow=c(1,4))
plot(cox.zph(fit1)) ### Figura 6.4 ###
H0<-basehaz(fit1, centered=F); H0 # risco acumulado de base
H01<-as.matrix(H0[1:81,1]) # risco acumulado de base: estrato 1
H02<-as.matrix(H0[82:101,1]) # risco acumulado de base: estrato 2
tempo1<- H0$time[1:81] # tempos: estrato 1
S01<-exp(-H01) # sobrevivência de base: estrato 1
round(cbind(tempo1,S01,H01),digits=5) # funções de base: estrato 1
tempo2<- H0$time[82:101] # tempos: estrato 2
S02<-exp(-H02) # sobrevivência de base: estrato 2
round(cbind(tempo2,S02,H02),digits=5) # funções de base: estrato 2
par(mfrow=c(1,2))
plot(tempo2[tempo2<=4], H02[tempo2<=4], lty=4, type="s", xlab="Tempo (anos)",
xlim=range(c(0,4)), ylab=expression(Lambda[0][j]* (t)))
lines(tempo1[tempo1<=4], H01[tempo1<=4], type="s", lty=1)
legend(0.0, 9.5, lty=c(1,4), c("Leuini < 75000","Leuini > 75000"),
lwd=1, bty="n", cex=0.8)
plot(c(0, tempo1), c(1,S01), lty=1, type="s", xlab="Tempo (anos)",
xlim=range(c(0,4)), ylab=expression(S[0][j]*(t)))
lines(c(0, tempo2), c(1,S02), lty=4, type="s")
legend(1.8, 0.9, lty=c(1,4), c("Leuini < 75000","Leuini > 75000"),
lwd=1,bty="n", cex=0.8)Seção 6.6 - Estudo sobre hormônio de crescimento
library(survival)
hg2<-read.table("https://docs.ufpr.br/~giolo/asa/dados/hg2.txt", h=T)
attach(hg2)
fit1<-coxph(Surv(tempos,cens)~factor(raca) + factor(trauma) + factor(recemnas) +
factor(renda) + ialtura + factor(trauma)*factor(recemnas),
data=hg2, method="breslow")
summary(fit1)
rendac<-ifelse(renda < 4, 1, 2)
fit2<-coxph(Surv(tempos,cens)~factor(raca) + factor(trauma) + factor(recemnas) +
factor(rendac) + ialtura + factor(trauma)*factor(recemnas),
data=hg2, method="breslow")
summary(fit2)
cox.zph(fit2, transform="identity")
### Figura 6.6 ###
par(mfrow=c(2,3))
plot(cox.zph(fit2))Modelo de Cox estratificado
hg2a<-na.omit(hg2[,c(1:5,7)])
rendac<-ifelse(hg2a$renda<4,1,2)
alt<-ifelse(hg2a$ialtura<120,1,2)
fit4<-coxph(Surv(tempos,cens)~factor(raca) + factor(trauma) + factor(rendac) +
strata(alt), data=hg2a, method="breslow")
summary(fit4)
cox.zph(fit4, transform="identity")
### Figura 6.7 ###
par(mfrow=c(1,3))
plot(cox.zph(fit4)) ### Figura 6.8 ###
resm<-resid(fit4, type="martingale")
res<-hg2a$cens-resm
ekm <- survfit(Surv(res, hg2a$cens)~1)
summary(ekm)
par(mfrow=c(1,2))
plot(ekm, mark.time=F,conf.int=F,xlab="Resíduos de Cox-Snell",ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(1,0.8, lty=c(1,3), c("Kaplan Meier","Exponencial padrão"),
lwd=1, bty="n", cex=0.9)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): exponencial padrão", pch=16)
abline(a=0,b=1, lty=3) ### Figura 6.9 ###
H0<-basehaz(fit4, centered=F); H0
H01<-as.matrix(H0[1:21,1])
H02<-as.matrix(H0[22:39,1])
tempo1<-H0$time[1:21]
S01<-exp(-H01)
round(cbind(tempo1,S01,H01), digits=5)
tempo2<- H0$time[22:39]
S02<-exp(-H02)
round(cbind(tempo2,S02,H02), digits=5)
par(mfrow=c(1,2))
plot(tempo2[H02<8], H02[H02<8], lty=4, type="s", xlab="Tempo (em meses)",
ylim=c(0,6),xlim=range(c(10,50)), ylab=expression(Lambda[0][j]*(t)))
lines(tempo1[H01<6], H01[H01<6], type="s", lty=1)
legend(10,5.6, lty=c(1,4), c(expression("Altura inicial"<120~cm),
expression("Altura inicial">=120~cm)),bty="n",cex=0.9)
plot(c(0,tempo2),c(1,S02), lty=4, type="s", xlab="Tempo (em meses)",
ylim=range(c(0,1)), xlim=range(c(10,50)), ylab=expression(S[0][j]*(t)))
lines(c(0,tempo1),c(1,S01), lty=1, type="s")
legend(22,0.95, lty=c(1,4), c(expression("Altura inicial"<120~cm),
expression("Altura inicial">=120~cm)),bty="n",cex=0.9)Capítulo 7
Modelo aditivo de Aalen
Seção 7.8.1 - Estudo sobre câncer de laringe
library(survival)
library(timereg)
laringe<-read.table("https://docs.ufpr.br/~giolo/asa/dados/laringe.txt",h=T)
attach(laringe)
laringe$idadec<-idade-mean(idade)
i<-order(laringe$tempos)
dados<-laringe[i,]## Obs: avaliar, inicialmente, as saídas sem fixar max.time para, então, defini-lo.
fit1<- aalen(Surv(tempos,cens)~factor(estadio)+idadec+factor(estadio)*idadec,dados)
par(mfrow=c(2,4)); plot(fit1)
## Note nos gráficos que a partir de t=3.9 as estimativas se repetem. Isto indica
## que X(t) deixa de ter posto completo em t = 3.9 ==> tau = max.time = 3.9.
fit1<- aalen(Surv(tempos,cens)~factor(estadio)+idadec+factor(estadio)*idadec,
dados, max.time=3.9)
summary(fit1); par(mfrow=c(2,4)); plot(fit1)
## Obs: Como dito no texto, a abordagem de reamostragem é utilizada para obter a
## distribuição das estatísticas de teste. Assim, os valores p não são exatamente
## os mesmos para cada vez que os comandos acima são executados. ## Obs: procedendo do mesmo modo, tem-se tau = 4.3 para fit2 e fit3.
fit2<- aalen(Surv(tempos,cens)~factor(estadio)+idadec,dados,max.time=4.3)
summary(fit2); par(mfrow=c(2,3)); plot(fit2)
fit3<- aalen(Surv(tempos,cens)~factor(estadio),dados,max.time=4.3)
summary(fit3)
par(mfrow=c(2,2)); plot(fit3) # Figura 7.3
## Nota: alternativamente, pode-se utilizar a função Addreg.r para
## encontrar tau = max.time. Isto é ilustrado no próximo exemplo. ### Figura 7.4 ###
options(OutDec = ",")
par(mfrow=c(1,1))
plot(fit3$cum[,1],fit3$cum[,2], typ="s", ylim=c(0,2),lwd=1, xlab="Tempo (anos)",
ylab = expression(hat(Lambda)(t~"|"~x)), cex.lab=1.2)
lines(fit3$cum[,1],(fit3$cum[,2]+fit3$cum[,3]), typ="s", ylim=c(0,2),lwd=1, lty=2)
lines(fit3$cum[,1],(fit3$cum[,2]+fit3$cum[,4]), typ="s", ylim=c(0,2),lwd=1, lty=4)
lines(fit3$cum[,1],(fit3$cum[,2]+fit3$cum[,5]), typ="s", ylim=c(0,2),lwd=1, lty=4)
lines(fit3$cum[,1],(fit3$cum[,2]+fit3$cum[,5]), typ="s", ylim=c(0,2),lwd=1, lty=5)
legend(0,1.8, c("Estádio I", "Estádio II", "Estádio III", "Estádio IV"),
lty=c(1,2,4,5), bty="n", cex =1.1)Seção 7.8.2 - Estudo sobre sinusite em pacientes com HIV
aids<-read.table("https://docs.ufpr.br/~giolo/asa/dados/aids.txt",h=T)
aids1<-subset(aids,ti<tf)
aids1$idc<- aids1$id-mean(aids1$id[!is.na(aids1$id)])
attach(aids1)
aids2<-aids1[,c(1:7,13)]
aids3<-na.omit(aids2)
## Utilizando a função Addreg.r para encontrar o tempo máximo = tau.
library(survival)
source("https://docs.ufpr.br/~giolo/asa/arquivos/Addreg.r")
fit1<-addreg(Surv(ti,tf,cens) ~ idc + factor(grp), aids3)
summary(fit1)
## Note nas saídas: "(Last estimate at time 617)" ==> tau = 617.
library(timereg)
fit3<-aalen(Surv(ti,tf,cens) ~ idc + factor(grp), aids3, max.time=617)
par(mfrow=c(2,3))
plot(fit3, cex.lab=1.2, cex.axis=1.2)### Figura 7.6 ###
fit5<-aalen(Surv(ti,tf,cens) ~ idc + factor(grp), aids3,
max.time=617, residuals=1)
n<-nrow(aids3)
rm<-matrix(0,n,1) # resíduos martingal
for(i in 1:n){
rm[i]<-sum(fit5$residuals$dM[,i])
}
delta<-aids3$cens
ei<-delta-rm # resíduos Cox-Snell
r.surv <- survfit(Surv(ei,delta)~1)
e<-r.surv$time
He<- -log(r.surv$surv)
par(mfrow=c(1,2))
plot(e[e>=0], He[e>=0], xlab=expression("Resíduos de Cox-Snell"~"="~e[i]),
ylab= expression(hat(Lambda)(e[i])), ylim=c(0,2), pch=16, xlim=c(0,2))
abline(a=0,b=1, lty=3)
title(main="(a)")
plot(e[e>=0], He[e>=0],type="s", xlab=expression("Resíduos de Cox-Snell"~"="~e[i]),
ylab= expression(hat(Lambda)(e[i])), ylim=c(0,2), pch=16, xlim=c(0,2))
abline(a=0,b=1, lty=3)
title(main="(b)")### Figura 7.7 ###
par(mfrow=c(1,2))
plot(fit3$cum[,1], fit3$cum[,2], type="s", ylim=c(0,1.6), xlab="Tempo (dias)",
ylab= expression(hat(Lambda)("t|x(t)")))
lines(fit3$cum[,1], (fit3$cum[,2]+fit3$cum[,4]), type="s", lty=3)
lines(fit3$cum[,1], (fit3$cum[,2]+fit3$cum[,5]), type="s", lty=2)
lines(fit3$cum[,1], (fit3$cum[,2]+fit3$cum[,6]), type="s", lty=4)
legend(0,1.5, c("HIV soronegativo","HIV assintomático","Grupo ARC","Grupo AIDS"),
lty=c(1,3,2,4), bty="n", cex =1.1)
title("(a)" )
plot(fit3$cum[1:6,1], fit3$cum[1:6,2], type="s", ylim=c(0,1.6), xlim=c(0,617),
xaxt="n", xlab="Tempo (dias)", ylab=expression(hat(Lambda)("t|x(t)"))) #(0;t1]
lines(fit3$cum[6:15,1],(fit3$cum[6:15,2]+fit3$cum[6:15,4]),
type="s",lty=2) #(t1;t2]
lines(fit3$cum[15:24,1],(fit3$cum[15:24,2]+fit3$cum[15:24,5]),
type="s",lty=1) #(t2;t3]
lines(fit3$cum[24:27,1],(fit3$cum[24:27,2]+fit3$cum[24:27,6]),
type="s",lty=2) #(t3;tau]
lines(c(209.5, 209.5), c(0.024,0.34), lty=3, col=1)
lines(c(415,415), c(0.77,1.02), lty=3, col=1)
axis(1, c(59.5, 209.5, 415.0, 617), c(expression(t[1]),expression(t[2]),
expression(t[3]),expression(tau)))
axis(1, c(100,300, 500), c(100, 300, 500))
text(150,1.4, expression("HIV soronegativo em (0,"~t[1]~"]"))
text(160,1.3, expression("HIV assintomático em ("~t[1]~","~t[2]~"]"))
text(131,1.2, expression("Grupo ARC em ("~t[2]~","~t[3]~"]"))
text(131,1.1, expression("Grupo AIDS em ("~t[3]~","~tau~"]"))
title("(b)" )library(timereg)
# ajuste do modelo aditivo semiparamétrico
fit4<-aalen(Surv(ti,tf,cens)~const(idc)+factor(grp),max.time=617,aids3)
summary(fit4); par(mfrow=c(2,2)); plot(fit4)
# ajuste do modelo aditivo de Lin e Ying
fit5<-aalen(Surv(ti,tf,cens)~const(idc)+const(factor(grp)),max.time=617,aids3)
summary(fit5); par(mfrow=c(1,1)); plot(fit5)Capítulo 8
Censura intervalar e dados grupados
Seção 8.2 - Estimador não paramétrico de Turnbull
library(survival)
source("http://docs.ufpr.br/~giolo/asa/arquivos/Turnbull.R")
left<-c(0,1,4,5,5); right<-c(5,8,9,8,9)
dat<-as.data.frame(cbind(left,right))
attach(dat)
right[is.na(right)] <- Inf
tau <- cria.tau(dat)
p <- S.ini(tau=tau)
A <- cria.A(data=dat,tau=tau)
tb <- Turnbull(p,A,dat); tbSeção 8.2.1 - Estudo sobre Câncer de Mama
rm(list=ls())
library(survival)
source("https://docs.ufpr.br/~giolo/asa/arquivos/Turnbull.R")
dat <- read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt",h=T)
dat1 <- dat[dat$ther==1,]
dat1$right[is.na(dat1$right)] <- Inf
tau <- cria.tau(dat1)
p <- S.ini(tau=tau)
A <- cria.A(data=dat1,tau=tau)
tb1 <- Turnbull(p,A,dat1)
tb1
dat1 <- dat[dat$ther==0,]
dat1$right[is.na(dat1$right)] <- Inf
tau <- cria.tau(dat1)
p <- S.ini(tau=tau)
A <- cria.A(data=dat1,tau=tau)
tb2 <- Turnbull(p,A,dat1)
tb2
### Figura 8.1 ###
par(mfrow=c(1,1))
plot(tb1$time, tb1$surv, lty=1, type = "s", ylim=c(0,1), xlim=c(0,50),
xlab="Tempo (em meses)", ylab="S(t)")
lines(tb2$time, tb2$surv, lty=4, type="s")
legend(1,0.3,lty=c(1,4),c("Radioterapia","Radioterapia e quimioterapia"),
bty="n",cex=0.9) ### Figura 8.2 ###
p <-dat$left+((dat$right-dat$left)/2)
pm <-ifelse(is.finite(p),p,dat$left)
cens <- ifelse(is.finite(p),1,0)
ekm<-survfit(Surv(pm,cens)~ther, type=c("kaplan-meier"), data=dat)
plot(tb1$time, tb1$surv, lty=2, type="s", ylim=c(0,1), xlim=c(0,50),
xlab="Tempo (em meses)", ylab="S(t)", cex.lab=1.0)
lines(tb2$time, tb2$surv, lty=2, type="s")
lines(ekm[1]$time, ekm[1]$surv, type="s", lty=1, lwd=1)
lines(ekm[2]$time, ekm[2]$surv, type="s", lty=1, lwd=1)
legend(0,0.2, lty=c(2,1), c("Turnbull","Ponto médio dos intervalos"),
lwd=1, bty="n", cex=1.0)
legend(39,0.6,c("Radioterapia"), bty="n", cex=1.0)
legend(37,0.28,c("Radioterapia e \n quimioterapia"), bty="n", cex=1.0)Seção 8.3.1 - Análise paramétrica dos dados de câncer de mama
breast<-read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt", h=T)
breast$left<-ifelse(breast$left==0,0.01,breast$left)
attach(breast)
library(survival)
fit1<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="weibull")
summary(fit1)
fit2<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="lognorm")
summary(fit2) ### Figura 8.3 ###
par(mfrow=c(1,2))
plot(tb1$time,tb1$surv, lty=1,type="s", ylim=c(0,1), xlim=c(0,50),
xlab="Tempo (em meses)", ylab="S(t)")
lines(tb2$time,tb2$surv, lty=2, type="s")
t1<-0:50
b0<-fit1$coefficients[1]
b1<-fit1$coefficients[2]
s<- 1/fit1$scale
a1<- exp(b0+b1)
st1<-exp(-(t1/a1)^s)
t2<-0:50
a2<- exp(b0)
st2<-exp(-(t2/a2)^s)
lines(t1,st1, type="l", lty=4)
lines(t2,st2, type="l", lty=4)
legend(0,0.25, lty=c(1,2,4), c("Turnbull: radioterapia",
"Turnbull: radioterapia e quimioterapia",
"Modelo de Weibull"), lwd=1, bty="n", cex=0.9)
plot(tb1$time,tb1$surv, lty=1,type="s", ylim=c(0,1), xlim=c(0,50),
xlab="Tempo (em meses)", ylab="S(t)")
lines(tb2$time,tb2$surv, lty=2, type="s")
t1<-0:50
b0<-fit2$coefficients[1]
b1<-fit2$coefficients[2]
s<- fit2$scale
a1<- (b0+b1)
st1<-pnorm((-log(t1) + a1)/s)
t2<-0:50
a2<- b0
st2<-pnorm((-log(t2) + a2)/s)
lines(t1,st1, type="l", lty=4)
lines(t2,st2, type="l", lty=4)
legend(0,0.25, lty=c(1,2,4), c("Turnbull: radioterapia",
"Turnbull: radioterapia e quimioterapia",
"Modelo log-normal"), lwd=1, bty="n", cex=0.9)## Figura 8.4
breast<-read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt", h=T)
breast$left<-ifelse(breast$left==0,0.01,breast$left)
attach(breast)
library(survival)
fit1<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="weibull")
summary(fit1)
fit2<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="lognorm")
summary(fit2)
# Resíduos de Cox-Snell ==> modelo Weibull
xb1<-fit1$coefficients[1]+fit1$coefficients[2]*ther
gama<-1/fit1$scale; gama
resl<-((breast$left)*exp(-xb1))^gama # resíduos (limite inferior)
resu<-((breast$right)*exp(-xb1))^gama # resíduos (limite superior)
left<-resl; right<-resu
source("https://docs.ufpr.br/~giolo/asa/arquivos/Turnbull.R")
dat1<-as.data.frame(cbind(left,right))
dat1$right[is.na(dat1$right)] <- Inf
tau <- cria.tau(dat1)
p <- S.ini(tau=tau)
A <- cria.A(data=dat1,tau=tau)
tb1 <- Turnbull(p,A,dat1); tb1
tw<-tb1$time
rcsw<- -log(tb1$surv)
par(mar=c(4.5,5,2.5,2)) # c(bottom,left,top,right)
par(mfrow=c(1,2))
plot(rcsw,tw, ylab=expression(-log(hat(S)(tau[j]))), ylim=c(0,3), xlim=c(0,3),
xlab=expression("Resíduos"~tau[j]), pch=16, cex.lab=1.3)
lines(abline(0,1))
title("(a) Weibull")
# Resíduos de Cox-Snell ==> modelo log-normal
xb<-fit2$coefficients[1]+fit2$coefficients[2]*ther
sigma<-fit2$scale
resl<-(log(breast$left)-(xb))/sigma # resíduos (limite inferior)
resu<-(log(breast$right)-(xb))/sigma # resíduos (limite superior)
el<- -log(1-pnorm(resl)); eu<- -log(1-pnorm(resu))
left<-el; right<-eu
dat<-as.data.frame(cbind(left,right))
dat$right[is.na(dat$right)] <- Inf
tau <- cria.tau(dat)
p <- S.ini(tau=tau)
A <- cria.A(data=dat,tau=tau); dim(A)
tb <- Turnbull(p,A,dat); tb
tl<-tb$time
rcsl<--log(tb$surv)
plot(rcsl,tl, ylim=c(0,3), xlim=c(0,3), ylab=expression(-log(hat(S)(tau[j]))),
xlab=expression("Resíduos"~tau[j]), pch=16, cex.lab=1.3)
lines(abline(0,1))
title("(b) Log-normal")
# tempo mediano
sigma<-exp(-0.479)
exp(3.331+0.568*1)*(-log(1-0.5))^sigma # radioterapia (x=1)
exp(3.331+0.568*0)*(-log(1-0.5))^sigma # radioterapia e quimioterapia (x=0)Seção 8.4.1 - Ilustração do modelo de Cox
rm(list=ls())
library(icenReg)
breast<-read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt", h=T)
m<-max(summary(breast$right))
li<-breast$left
ui<-ifelse(is.na(breast$right), m+1000, breast$right)
attach(breast)
fit1<- ic_sp(cbind(li,ui) ~ ther, model='ph', bs_samples=100, data=breast)
summary(fit1) ### Figura 8.5 ###
rm(list=ls())
library(survival)
source("https://docs.ufpr.br/~giolo/asa/arquivos/Turnbull.R")
breast<-read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt",h=T)
dat1 <- breast[breast$ther==1,]
dat1$right[is.na(dat1$right)] <- Inf
tau <- cria.tau(dat1); p <- S.ini(tau=tau)
A <- cria.A(data=dat1,tau=tau)
tb1 <- Turnbull(p,A,dat1); tb1
lls1<- log(-log(tb1$surv)); lls1
dat2 <- breast[breast$ther==0,]
dat2$right[is.na(dat2$right)] <- Inf
tau <- cria.tau(dat2); p <- S.ini(tau=tau)
A <- cria.A(data=dat2,tau=tau)
tb2 <- Turnbull(p,A,dat1); tb2
lls2<- log(-log(tb2$surv)); lls2
par(mfrow=c(1,1))
plot(tb2$time,lls2, type="s", lty=2, ylim=c(-9,1.5), xlab="Tempo (em meses)",
ylab= expression(log(Lambda[0](t))))
lines(tb1$time,lls1, type="s", lty=1)
legend(20,-4, lty=c(1,2), c("Radioterapia","Radioterapia e quimioterapia"),
bty="n", cex=1.0) ### Figura 8.6 ###
breast <- read.table("https://docs.ufpr.br/~giolo/asa/dados/breast.txt", h=T)
m<-max(summary(breast$right))
li<-breast$left
ui<-ifelse(is.na(breast$right), m+1000, breast$right)
attach(breast)
source("https://docs.ufpr.br/~giolo/asa/arquivos/Turnbull.R")
dat<-as.data.frame(cbind(left,right))
dat$right[is.na(dat$right)] <- Inf
tau <- cria.tau(dat); p <- S.ini(tau=tau)
A <- cria.A(data=dat,tau=tau)
tb <- Turnbull(p,A,dat); tb # S_0(tau)
dat1 <- breast[breast$ther==1,]
dat1$right[is.na(dat1$right)] <- Inf
tau <- cria.tau(dat1); p <- S.ini(tau=tau)
A <- cria.A(data=dat1,tau=tau)
tb1 <- Turnbull(p,A,dat1); tb1 # S_R(tau) -> Turnbull radioterapia
dat2 <- breast[breast$ther==0,]
dat2$right[is.na(dat2$right)] <- Inf
tau <- cria.tau(dat2); p <- S.ini(tau=tau)
A <- cria.A(data=dat2,tau=tau)
tb2 <- Turnbull(p,A,dat2); tb2 # S_RQ(tau) -> Turnbull radioterapia e quimio
breast$trat<-ifelse(ther==1,1,-1)
fit2<- ic_sp(cbind(li,ui) ~ trat, model='ph', bs_samples=100, data=breast)
s1<-(tb$surv)^exp(-fit2$coef) # S_0(tau)^exp(-beta_1) -> modelo de Cox: radio e quimio
s2<-(tb$surv)^exp(fit2$coef) # S_0(tau)^exp(beta_1) -> modelo de Cox: radio
par(mfrow=c(1,1))
plot(tb1$time, tb1$surv, lty=2, type="s", ylim=c(0,1), xlim=c(0,50),
xlab="Tempo (em meses)", ylab="S(t)", cex.lab=1.0)
lines(tb2$time, tb2$surv, lty=2, type="s")
lines(tb$time, s1, type="s", lty=1)
lines(tb$time, s2, type="s", lty=1)
legend(0,0.2, lty=c(2,1), c("Turnbull","Modelo de Cox"), lwd=1, bty="n", cex=1.0)
legend(39,0.6,c("Radioterapia"), bty="n", cex=1.0)
legend(26.5,0.28,c("Radioterapia e \n quimioterapia"), bty="n", cex=1.0)Seção 8.6 - Exemplo: tempos de vida de mangueiras
mang<-read.table("https://docs.ufpr.br/~giolo/asa/dados/mang.txt", h=T)
attach(mang)
library(survival)
ekm<-survfit(Surv(ti,cens)~1, conf.type="none")
summary(ekm)
mang1<-read.table("https://docs.ufpr.br/~giolo/asa/dados/dadmang.txt", h=T)
attach(mang1)
library(survival)
fit1<-glm(y~ -1 + int1 + int2 + int3 + int4 + int5 + int6 + int7 + int8 +
int9 + int10 + int11 + int12 + factor(bloco,levels=5:1) + factor(copa) +
factor(cavalo) + factor(copa)*factor(cavalo),
family=binomial(link="cloglog"))
anova(fit1, test="Chisq")
fit2 <-glm(y~ -1 + int1 + int2 + int3 + int4 + int5 + int6 + int7 + int8 +
int9 + int10 + int11 + int12 + factor(bloco,levels=5:1) + factor(copa) +
factor(cavalo) + factor(copa)*factor(cavalo),
family=binomial(link="logit"))
anova(fit2, test="Chisq")
fit1<-glm(y~ -1 + int1 + int2 + int3 + int4 + int5 + int6 + int7 + int8 +
int9 + int10 + int11 + int12 + factor(bloco,levels=5:1) + factor(copa),
family=binomial(link="cloglog"))
summary(fit1)
fit2<-glm(y~ -1 + int1 + int2 + int3 + int4 + int5 + int6 + int7 + int8 +
int9 + int10 + int11 + int12 + factor(bloco,levels=5:1) + factor(copa),
family=binomial(link="logit"))
summary(fit2) ### Figura 8.7 ####
mang1<-read.table("https://docs.ufpr.br/~giolo/asa/dados/dadmang.txt", h=T)
attach(mang1)
fit1<-glm(y~-1+int1+int2+int3+int4+int5+int6+int7+int8+int9+int10+int11+int12+
factor(bloco,levels=5:1)+as.factor(copa), family=binomial(link="cloglog"))
cf<-as.vector(fit1$coefficients[1:12]) # (gama_i)*
gi<-exp(-exp(cf)) # (gama_i)
S0<-gi
for(i in 1:11){
S0[i+1]<-prod(gi[1:(i+1)])}
S0<-c(1,S0) # S0(t)
cf1<-fit1$coefficients[18:22]
Sc1<-S0; Sc2<-(S0)^exp(cf1[1])
Sc3<-(S0)^exp(cf1[2]); Sc4<-(S0)^exp(cf1[3])
Sc5<-(S0)^exp(cf1[4]); Sc6<-(S0)^exp(cf1[5])
t<-c(0,2,3,4,10,12,14,15,16,17,18,19,21)
cbind(t,Sc1,Sc2,Sc3,Sc4,Sc5,Sc6)
par(mfrow=c(1,1))
plot(t,Sc1, type="s", lty = 1, ylim=range(c(0,1)), xlab="Tempo (em anos)",
ylab="Sobrevivência estimada", cex.lab=1.2)
points(t,Sc1,pch=21);
lines(t,Sc2,type="s",lty=2); points(t,Sc2,pch=15)
lines(t,Sc3,type="s",lty=3); points(t,Sc3,pch=14)
lines(t,Sc4,type="s",lty=4); points(t,Sc4,pch=8)
lines(t,Sc5,type="s",lty=5); points(t,Sc5,pch=16)
lines(t,Sc6,type="s",lty=6); points(t,Sc6,pch=17)
legend(1,0.5,lty=c(1,2,3,4,5,6),pch=c(21,15,14,8,16,17),
c("Copa 1 - Extrema","Copa 2 - Oliveira","Copa 3 - Pahiri",
"Copa 4 - Imperial","Copa 5 - Carlota","Copa 6 - Bourbon"),
bty="n", cex=1.2)
title("Modelo taxas de falha proporcionais de Cox", cex.main = 1.2)### Figura 8.8 ###
mang1<-read.table("https://docs.ufpr.br/~giolo/asa/dados/dadmang.txt", h=T)
attach(mang1)
fit2<-glm(y~-1+int1+int2+int3+int4+int5+int6+int7+int8+int9+int10+int11+int12+
factor(bloco,levels=5:1)+as.factor(copa), family=binomial(link="logit"))
cf<-fit2$coefficients[1:12] #(gama_i)*
gi<-exp(cf) #(gama_i)
cf1<-fit2$coefficients[18:22]
qi1<-(1/(1+gi))
Slc1<-qi1
for(i in 1:11){
Slc1[i+1]<-prod(qi1[1:(i+1)])}
Slc1<-c(1,Slc1) #S(t) copa 1
qi2<-(1/(1+gi*exp(cf1[1])))
Slc2<-qi2
for(i in 1:11){
Slc2[i+1]<-prod(qi2[1:(i+1)])}
Slc2<-c(1,Slc2) #S(t) copa 2
qi3<-(1/(1+gi*exp(cf1[2])))
Slc3<-qi3
for(i in 1:11){
Slc3[i+1]<-prod(qi3[1:(i+1)])}
Slc3<-c(1,Slc3) #S(t) copa 3
qi4<-(1/(1+gi*exp(cf1[3])))
Slc4<-qi4
for(i in 1:11){
Slc4[i+1]<-prod(qi4[1:(i+1)])}
Slc4<-c(1,Slc4) #S(t) copa 4
qi5<-(1/(1+gi*exp(cf1[4])))
Slc5<-qi5
for(i in 1:11){
Slc5[i+1]<-prod(qi5[1:(i+1)])}
Slc5<-c(1,Slc5) #S(t) copa 5
qi6<-(1/(1+gi*exp(cf1[5])))
Slc6<-qi6
for(i in 1:11){
Slc6[i+1]<-prod(qi6[1:(i+1)])}
Slc6<-c(1,Slc6) #S(t) copa 6
t<-c(0,2,3,4,10,12,14,15,16,17,18,19,21)
cbind(t,Slc1,Slc2,Slc3,Slc4,Slc5,Slc6)
par(mfrow=c(1,1))
plot(t, Slc1, type="s", lty=1, ylim=range(c(0,1)), xlab="Tempo (em anos)",
ylab="Sobrevivência estimada", cex.lab=1.2)
points(t, Slc1, pch=21)
lines(t, Slc2, type="s", lty=2); points(t, Slc2, pch=15)
lines(t, Slc3, type="s", lty=3); points(t, Slc3, pch=14)
lines(t, Slc4, type="s", lty=4); points(t, Slc4, pch=8)
lines(t, Slc5, type="s", lty=5); points(t, Slc5, pch=16)
lines(t, Slc6, type="s", lty=6); points(t, Slc6, pch=17)
legend(1, 0.5, lty=c(1,2,3,4,5,6), pch=c(21,15,14,8,16,17),
c("Copa 1 - Extrema","Copa 2 - Oliveira","Copa 3 - Pahiri",
"Copa 4 - Imperial","Copa 5 - Carlota","Copa 6 - Bourbon"),
bty="n", cex=1.2)
title("Modelo logístico", cex.main=1.2)Capítulo 9
Riscos competitivos
Seção 9.3.3 - Exemplo de estimação da FIA
library(survival)
tempo<-c(1,3,4,5,7,7,8,9,12,12,13,15,17,25,34,40,52,63,66,71)
evento<-c(1,2,1,2,1,1,2,1,2,1,2,1,2,2,0,2,2,1,2,2) # 1=óbito e 2=alta
KMo<-survfit(Surv(tempo,evento==1)~1) # KM evento 1 = óbito
summary(KMo)
1-KMo$surv
cens<-ifelse(evento>=1,1,0)
KMg<-survfit(Surv(tempo,cens)~1) # Sobrevida global
summary(KMg)
fia<-survfit(Surv(tempo,evento,type="mstate")~1)
summary(fia)
### KM para o evento alta hospitalar ###
KMa<-survfit(Surv(tempo,evento==2)~1) # KM evento 2 = alta
summary(KMa)
1-KMa$surv### Figura 9.2 ###
rm(list=ls())
library(cmprsk)
tempo<-c(1,3,4,5,7,7,8,9,12,12,13,15,17,25,34,40,52,63,66,71)
evento<-c(1,2,1,2,1,1,2,1,2,1,2,1,2,2,0,2,2,1,2,2)
da<-data.frame(tempo, evento)
fia1<-with(da,cuminc(tempo,evento,rho=0,cencode=0))
par(mfrow=c(1,2))
plot(survfit(Surv(tempo,evento==1)~1), conf.int=F, fun=function(x){1-x},
ylim=c(0,1), col=1, lty=1, lwd=1, bty="n", xlab="Tempo (em dias)",
ylab=expression(Kaplan-Meier: 1-S[k](t)), main="(a)")
lines(survfit(Surv(tempo,evento==2)~1), conf.int=F, fun=function(x){1-x},
col=1, lty=2, lwd=1)
legend(0,1, c("Alta hospitalar","Óbito"), lty=c(1,2), col=c(1,1),
bty='n', cex=1.0)
plot(fia1, col=c(1,1), lty=c(1,2), lwd=c(1,1),
xlab="Tempo (em dias)", ylab=expression(FIA: F[k](t)),
curvlab=c("Alta hospitalar","Óbito"), main="(b)")### IC 95% ###
fia<-survfit(Surv(tempo,evento,type="mstate")~1,conf.type="log-log")
cbind(fia$time,fia$lower[,2],fia$upper[,2]) # evento 1### Figura 9.3 ###
rm(list=ls())
library(survival)
library(cmprsk)
covid<-read.table("https://docs.ufpr.br/~giolo/asa/dados/covid.txt",h=T)
attach(covid)
a.surv<-survfit(Surv(tempo,desfecho,type="mstate")~comorbidades,data=covid)
par(mfrow=c(1,1))
plot(a.surv,ylim=c(0,1),xlab="Tempo (em dias)",xlim=c(0,55),ylab='FIA',
cex=1,conf.int=F,lty=c(1,4, 1,4))
legend(0,1.0, lty=c(1,4), c("Sem comorbidades","Com comorbidades"),
bty='n',cex=1.0)
text(49,0.20, "Óbito", cex=0.8)
text(49,0.10, "Óbito", cex=0.8)
text(47.5,0.96, "Alta hospitalar", cex=0.8)
text(47.5,0.86, "Alta hospitalar", cex=0.8)
### Teste de Gray
com_test<-with(covid,cuminc(ftime=tempo,fstatus=desfecho,group=comorbidades))
com_test$Tests
### Teste logrank usual ###
survdiff(Surv(tempo,desfecho==1)~comorbidades,data=covid) # 1=óbito
survdiff(Surv(tempo,desfecho==2)~comorbidades,data=covid) # 2=altaSeção 9.6 - Exemplo Covid-19 em crianças e adolescentes
### Figura 9.5 ###
covid<-read.table("https://docs.ufpr.br/~giolo/asa/dados/covid.txt",h=T)
attach(covid)
fia1<-survfit(Surv(tempo,desfecho,type="mstate")~sexo,data=covid)
tg1<-with(covid, cuminc(ftime=tempo, fstatus=desfecho, group=sexo))
tg1$Tests
par(mfrow=c(1,3))
plot(fia1,ylim=c(0,1),xlab="Tempo (dias)",xlim=c(0,55),ylab='FIA estimada',
cex.lab=1.3,conf.int=F,lty=c(1,4, 1,4))
legend(0,1.0, lty=c(1,4), c("Masculino", "Feminino"), bty='n',cex=1.3)
text(45,0.15,"Óbito", cex=1.3)
text(45,0.94,"Alta hospitalar", cex=1.3)
text(45,0.83,"valor p = 0,68", cex=1.3)
text(45,0.05,"valor p = 0,54", cex=1.3)
title("Gênero", cex.main=1.3)
fia2<-survfit(Surv(tempo,desfecho,type="mstate")~idade,data=covid)
tg2<-with(covid, cuminc(ftime=tempo, fstatus=desfecho, group=idade))
tg2$Tests
plot(fia2,ylim=c(0,1),xlab="Tempo (dias)",xlim=c(0,55),ylab='FIA estimada',
cex.lab=1.3,conf.int=F,lty=c(1,2,4, 1,2,4))
legend(0,1.0, lty=c(1,2,4), c("[0, 2)", "[2, 11)", "[11, 20)"), bty='n',cex=1.3)
text(46,0.17,"Óbito", cex=1.3)
text(45,0.96,"Alta hospitalar", cex=1.3)
text(45,0.83,"valor p = 0,04", cex=1.3)
text(45,0.03,"valor p = 0,28", cex=1.3)
title("Idade (anos)", cex.main=1.3)
fia3<-survfit(Surv(tempo,desfecho,type="mstate")~comorbidades,data=covid)
tg3<-with(covid, cuminc(ftime=tempo, fstatus=desfecho, group=comorbidades))
tg3$Tests
plot(fia3,ylim=c(0,1),xlab="Tempo (dias)",xlim=c(0,55),ylab='FIA estimada',
cex.lab=1.3,conf.int=F,lty=c(1,4, 1,4))
legend(0,1.0, lty=c(1,4), c("Sem", "Com"), bty='n',cex=1.3)
text(46,0.21,"Óbito", cex=1.3)
text(45,0.96,"Alta hospitalar", cex=1.3)
text(45,0.76,"valor p = 0,0008", cex=1.3)
text(45,0.02,"valor p = 0,0003", cex=1.3)
title("Comorbidades", cex.main=1.3)Seção 9.6.1 - Modelos para o evento óbito
# Modelo taxas de falha causa-específica
library(survival)
covid<-read.table("https://docs.ufpr.br/~giolo/asa/dados/covid.txt",h=T)
attach(covid)
covid$idade <- relevel(factor(covid$idade), ref = 2)
mod<- coxph(Surv(tempo,desfecho==1)~factor(idade)+sexo+comorbidades,data=covid)
mod
mod1<- coxph(Surv(tempo,desfecho==1)~factor(idade)+comorbidades,data=covid)
mod1
summary(mod1)
exp(cbind(RTF=coef(mod1),confint(mod1, level=0.95)))
exp(cbind(RTF=coef(mod1),confint(mod1, level=0.90)))
cox.zph(mod1, transform="km", terms=TRUE)
### Figura 9.6 ###
par(mfrow=c(1,3))
xx<-cox.zph(mod1, transform="km", terms=FALSE)
plot(xx[1], xlab="Tempo (em dias)",ylab=expression(Resíduos~Schoenfeld~+hat(beta[1])),
cex.lab=1.2, main="Idade [0,2) anos")
abline(h=0.6545, lty=3)
plot(xx[2], xlab="Tempo (em dias)",ylab=expression(Resíduos~Schoenfeld~+hat(beta[2])),
cex.lab=1.2, main="Idade [11,20) anos")
abline(h=0.5682, lty=3)
plot(xx[3], xlab="Tempo (dias)",ylab=expression(Resíduos~Schoenfeld~+hat(beta[3])),
cex.lab=1.2, main="Comorbidades: sim")
abline(h=0.8085, lty=3)# Modelo de Fine-Gray
# xmod1 = matriz de delinamento com as 3 covariáveis na forma matricial
cens<-ifelse(desfecho>1,0,desfecho)
mod1<-coxph(Surv(tempo,cens)~factor(idade)+factor(sexo)+factor(comorbidades),data=covid,x=TRUE)
xmod1<-model.matrix(mod1)
fit1<-with(covid, crr(tempo,desfecho,xmod1,failcode=1)) # 1=óbito
summary(fit1)
# xmod2 = Matriz de delinamento com 2 covariáveis na forma matricial (excluído sexo)
mod2<-coxph(Surv(tempo,cens)~factor(idade)+factor(comorbidades),data=covid,x=T)
xmod2<-model.matrix(mod2)
fit2<-with(covid, crr(tempo,desfecho,xmod2,failcode=1)) # 1 = óbito
summary(fit2)### Figura 9.7 ###
par(mar=c(4.5,5,2.5,2))
par(mfrow=c(1,3))
names(fit2$coef)=c("Idade [0, 2) anos","Idade [11, 20) anos","Comorbidades: sim")
for(j in 1:ncol(fit2$res)){
scatter.smooth(fit2$uft, fit2$res[,j], cex.lab=1.3, cex.main=1.4,
lpars = list(col='black', lwd=1),
main=names(fit2$coef)[j],
xlab='Tempo (em dias)', ylab='Resíduos tipo-Schoenfeld')
abline(h=0,lty=3)
} ### Figura 9.8 ###
par(mfrow=c(1,3))
fit1<-with(covid, crr(tempo,desfecho,xmod2, failcode=1)) # 1 = óbito
fp1<-predict(fit1, rbind(c(1,0,0), # idade [0,2), comorbidades = não
c(1,0,1))) # idade [0,2), comorbidades = sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,0.5), xlab="Tempo (em dias)",
ylab=expression(F[1]~"(t | x)"), main="Óbito", cex.lab=1.2)
legend(0,0.45, lty=c(1,4), col=c(1,1),
c("[0,2) anos sem comorbidades","[0,2) anos com comorbidades"),
bty="n", cex=1.2)
fp1<-predict(fit1, rbind(c(0,0,0), # idade [2,11), comorbidades = não
c(0,0,1))) # idade [2,11), comorbidades = sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,0.5), xlab="Tempo (em dias)",
ylab=expression(F[1]~"(t | x)"), main="Óbito", cex.lab=1.2)
legend(0,0.45, lty=c(1,4), col=c(1,1),
c("[2,11) anos sem comorbidades","[2,11) anos com comorbidades"),
bty="n", cex=1.2)
fp1<-predict(fit1, rbind(c(0,1,0), # idade [11,20), comorbidades = não
c(0,1,1))) # idade [11,20), comorbidades = sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,0.5), xlab="Tempo (em dias)",
ylab=expression(F[1]~"(t | x)"), main="Óbito", cex.lab=1.2)
legend(0,0.45, lty=c(1,4), col=c(1,1),
c("[11,20) anos sem comorbidades","[11,20) anos com comorbidades"),
bty="n", cex=1.2)9.6.2 - Modelos para o evento alta hospitalar
# Modelo taxas de falha causa-específica
library(survival)
cc1<-coxph(Surv(tempo,desfecho==2)~factor(idade)+sexo+comorbidades,data=covid)
cc1
exp(cbind(RTF=coef(cc1),confint(cc1, level=0.90)))
cc2<-coxph(Surv(tempo,desfecho==2)~factor(idade)+comorbidades,data=covid)
cc2
exp(cbind(RTF=coef(cc2),confint(cc2, level=0.90)))
cox.zph(cc2, transform="km", terms=TRUE)
### Figura 9.9 ###
par(mar=c(4.5,5,2.5,2)) # c(bottom,left,top,right)
par(mfrow=c(1,3))
xx<-cox.zph(cc2, transform="km", terms=FALSE)
plot(xx[1], xlab="Tempo (em dias) [g(t)=km]",
ylab=expression(Resíduos~Schoenfeld~+hat(beta[1])),
cex.lab=1.3, main="Idade [2,11) anos")
abline(h=cc2$coef[1], lty=3)
plot(xx[2], xlab="Tempo (em dias) [g(t)=km]",
ylab=expression(Resíduos~Schoenfeld~+hat(beta[2])),
cex.lab=1.3, main="Idade [11,20) anos")
abline(h=cc2$coef[2], lty=3)
plot(xx[3], xlab="Tempo (em dias) [g(t)=km]",
ylab=expression(Resíduos~Schoenfeld~+hat(beta[3])),
cex.lab=1.3, main="Comorbidades: sim")
abline(h=cc2$coef[3], lty=3) ### Figura 9.10 ###
resm<-resid(cc2,type="martingale")
cens<-ifelse(desfecho==2,1,0)
res<-cens - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, cens)~1)
par(mfrow=c(1,2))
plot(ekm, mark.time=F,conf.int=F,xlab="Resíduos de Cox-Snell",ylab="S(e) estimada")
res<-sort(res)
exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(1,0.8, lty=c(1,3),c("Kaplan Meier","Exponencial padrão"),
lwd=1, bty="n", cex=1.0)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st,sexp1,xlab="S(e): Kaplan-Meier",ylab="S(e): exponencial padrão",pch=16)# Modelo de Fine-Gray
library(cmprsk)
## xmod1 = matriz de delinamento com as 3 covariáveis na forma matricial
cens<-ifelse(desfecho==2,1,0)
mod1<-coxph(Surv(tempo,cens)~factor(sexo)+factor(idade)+
factor(comorbidades),data=covid,x=TRUE)
xmod1<-model.matrix(mod1)
fit1 <- with(covid, crr(tempo, desfecho, xmod1, failcode=2)) # 2=alta
summary(fit1)
## xmod2 = matriz com 2 covariáveis na forma matricial (excluído gênero)
cens<-ifelse(desfecho>1,0,desfecho)
mod2<-coxph(Surv(tempo,cens)~factor(idade)+factor(comorbidades),data=covid,x=TRUE)
xmod2<-model.matrix(mod2)
fit2 <- with(covid, crr(tempo, desfecho, xmod2, failcode=2)) # 2=alta
summary(fit2)### Figura 9.11 ###
par(mar=c(4.5,5,2.5,2)) # c(bottom,left,top,right)
par(mfrow=c(1,3))
names(fit2$coef) = c("Idade [0, 2) anos","Idade [11, 20) anos","Comorbidades: sim")
for(j in 1:ncol(fit2$res)){
scatter.smooth(fit2$uft, fit2$res[,j], cex.lab=1.3, cex.main=1.4,
lpars = list(col='black', lwd=1),
main = names(fit2$coef)[j],
xlab = 'Tempo (em dias)', ylab = 'Resíduos tipo-Schoenfeld')
abline(h=0,lty=3)
}### Figura 9.12 -> log[-log(1-F2(t))] versus log(t)
par(mfrow=c(1,2))
fit<-cuminc(covid$tempo,covid$desfecho,covid$idade)
a<-timepoints(fit,times=covid$tempo)
fia<-t(a$est[4:6,])
llfia<-log(-log(1-fia))
matplot(log(unique(sort(covid$tempo))),llfia, pch=c(22,12,10),col=c(1,1,1),
xlab="log(t)", ylab="log[-log(1-F(t))]")
legend('topleft', pch=c(22,12,10), c("Grupo [2, 11) anos","Grupo [0, 2) anos","Grupo [11, 20) anos"),
col=c(1,1,1), bty='n',cex=1.0)
title("Grupos de idade")
fit<-cuminc(covid$tempo,covid$desfecho,covid$comorbidades)
a<-timepoints(fit,times=covid$tempo)
fia<-t(a$est[3:4,])
llfia<-log(-log(1-fia))
matplot(log(unique(sort(covid$tempo))),llfia,pch=c(22,12),col=c(1,1),
xlab="log(t)", ylab="log[-log(1-F(t))]")
legend('topleft', pch=c(22,12), c("Comorbidades: não","Comorbidades: sim"),
col=c(1,1), bty='n',cex=1.0)
title("Comorbidades")### Figura 9.13 ###
### empírica ###
comb1<-subset(covid,idade==1 & comorbidades=="nao")
fia1 <- with(comb1, cuminc(tempo, desfecho, rho=0, cencode = 0))
fia12 <- cbind(fia1$`1 2`$ time,fia1$`1 2`$ est)
comb2<-subset(covid,idade==2 & comorbidades=="nao")
fia2 <- with(comb2, cuminc(tempo, desfecho, rho=0, cencode = 0))
fia22 <- cbind(fia2$`1 2`$ time,fia2$`1 2`$ est)
comb3<-subset(covid,idade==3 & comorbidades=="sim")
fia3 <- with(comb3, cuminc(tempo, desfecho, rho=0, cencode = 0))
fia32 <- cbind(fia3$`1 2`$ time,fia3$`1 2`$ est)
### modelo ###
cens<-ifelse(desfecho==2,0,desfecho)
mod2<-coxph(Surv(tempo,cens)~factor(idade)+factor(comorbidades),data=covid,x=TRUE)
xmod2<- model.matrix(mod2)
fit2 <- with(covid, crr(tempo, desfecho, xmod2, failcode=2)) # 2=alta
fp1 <- predict(fit2, rbind(c(1,0,0))) # idade [0,2), comorbidades = não
fp2 <- predict(fit2, rbind(c(0,0,0))) # idade [2,11), comorbidades = não
fp3 <- predict(fit2, rbind(c(0,1,1))) # idade [11,20), comorbidades = sim
par(mfrow=c(1,3))
plot(fia12, type ="s", xlab="Tempo (em dias)", ylab="F(t|x)", cex.lab=1.4,
ylim=c(0,1), xlim=c(0,55), main="[0,2) anos sem comorbidades")
lines(fp1, type="s", lty=2)
legend(20,0.3, lty=c(1,2), c("Empírica","Modelo"), bty="n", cex=1.1)
plot(fia22, type ="s", xlab="Tempo (em dias)", ylab="F(t|x)", cex.lab=1.4,
ylim=c(0,1), xlim=c(0,55), main="[2,11) anos sem comorbidades")
lines(fp2, type="s", lty=2)
legend(20,0.3, lty=c(1,2), c("Empírica","Modelo"), bty="n", cex=1.1)
plot(fia32, type ="s", xlab="Tempo (em dias)", ylab="F(t|x)", cex.lab=1.4,
ylim=c(0,1), xlim=c(0,55), main="[11,20) anos com comorbidades")
lines(fp3, type="s", lty=2)
legend(20,0.3, lty=c(1,2), c("Empírica","Modelo"), bty="n", cex=1.1)### Figura 9.14 ###
fit2<- with(covid, crr(tempo, desfecho, xmod2, failcode=2)) # 2=alta
fp1<-predict(fit2, rbind(c(0,0,0), # idade: [0,2), comorbidades: não
c(0,0,1))) # idade: [0,2), comorbidades: sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,1), xlab="Tempo (em dias)",
ylab=expression(F[2]~"(t | x)"), main="Alta hospitalar", cex.lab=1.2)
legend(10,0.3, lty=c(1,4), col=c(1,1), c("[0,2) anos sem comorbidades",
"[0,2) anos com comorbidades"), bty="n", cex=1.2)
fp1<-predict(fit2, rbind(c(1,0,0), # idade: [2,11), comorbidades: não
c(1,0,1))) # idade: [2,11), comorbidades: sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,1), xlab="Tempo (em dias)",
ylab=expression(F[2]~"(t | x)"), main="Alta hospitalar", cex.lab=1.2)
legend(10,0.3, lty=c(1,4), col=c(1,1), c("[2,11) anos sem comorbidades",
"[2,11) anos com comorbidades"), bty="n", cex=1.2)
fp1<-predict(fit2, rbind(c(0,1,0), # idade: [11,19), comorbidades: não
c(0,1,1))) # idade: [11,19), comorbidades: sim
plot(fp1, lty=c(1,4), col=c(1,1), ylim=c(0,1), xlab="Tempo (em dias)",
ylab=expression(F[2]~"(t | x)"), main="Alta hospitalar", cex.lab=1.2)
legend(10,0.3, lty=c(1,4), col=c(1,1), c("[11,20) anos sem comorbidades",
"[11,20) anos com comorbidades"), bty="n", cex=1.2)Capítulo 10
Estudos com fração de imunes
Seção 10.4 - Estudo sobre câncer de pele
pele<-read.table("https://docs.ufpr.br/~giolo/asa/dados/pele.txt",h=T)
i<-order(pele$survtime)
pele<-pele[i,]
attach(pele)Figura 10.2 - Investigando a presença de imunes ao evento
library(survival)
par(mfrow=c(1,2))
ekm<-survfit(Surv(survtime,status)~1)
plot(ekm,lty=c(1), col=c(1), mark.time=F, conf.int=F, xlab="Tempo (em meses)",
ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
lines(c(98,98),c(0.0,0.86), lty=3)
points(98,0.0,pch=25, col=1)
plot(ekm,lty=c(1), col=c(1), mark.time=F, conf.int=F, xlab="Tempo (em meses)",
xlim=c(0,125),ylim=c(0.8,1),ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
lines(c(98,98),c(0.8,0.86), lty=3)
points(98,0.8,pch=25, col=1)Figura 10.3 - Análise exploratória das covariáveis
par(mfrow=c(2,2))
ekm1<-survfit(Surv(survtime,status)~Sexo)
plot(ekm1,lty=c(1,2), col=c(1,1), mark.time=F, xlab="Tempo (em meses)",
ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
legend(1,0.4,lty=c(1,2),col=c(1,1), c("Feminino","Masculino"),
lwd=1,bty="n",cex=1.0)
title("(a) Sexo", cex.main=1.0)
legend(65,0.4, "logrank = 8,8, g.l.=1", cex=1.0, bty="n")
legend(65,0.31, "valor p = 0,003", cex=1.0, bty="n")
# Definindo a idade em faixas etárias #
pele$age<-ifelse(pele$Idade<=39,"18-39",0)
pele$age<-ifelse(pele$Idade>39 & pele$Idade<60,"40-59",pele$age)
pele$age<-ifelse(pele$age==0,"60-80",pele$age)
table(pele$age)
ekm2<-survfit(Surv(survtime,status)~pele$age)
plot(ekm2,lty=c(1,2,4),col=c(1,1,1),mark.time=F, xlab="Tempo (em meses)",
ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
legend(1,0.4,lty=c(1,2,4),col=c(1,1,1),c("18 a 39 anos","40 a 59 anos","60 a 80 anos"),
lwd=1,bty="n",cex=1.0)
title("(b) Faixa etária", cex.main=1.0)
legend(65,0.4, "logrank = 5,6, g.l.=2", cex=1.0, bty="n")
legend(65,0.31, "valor p = 0,06", cex=1.0, bty="n")
ekm3<-survfit(Surv(survtime,status)~Estadio)
plot(ekm3,lty=c(1,2,3,4),col=c(1,1,1,1), mark.time=F, xlab="Tempo (em meses)",
ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
legend(1,0.4,lty=c(1,2,3,4),col=c(1,1,1,1),
c("Estádio 0/1","Estádio 2","Estádio 3","Estádio 4"),
lwd=1, bty="n",cex=1.0)
title("(c) Estadiamento", cex.main=1.0)
legend(63,0.4, "logrank = 139, g.l.=3", cex=1.0, bty="n")
legend(63,0.31, "valor p < 0,001", cex=1.0, bty="n")
ekm4<-survfit(Surv(survtime,status)~Trat)
plot(ekm4,lty=c(1,2),col=c(1,1), mark.time=F, xlab="Tempo (em meses)",
ylab="S(t) estimada", cex.lab=1.2, cex.axis=1.2)
legend(1,0.4,lty=c(1,2),col=c(1,1),c("Cirurgia","Quimioterapia"),
lwd=1, bty="n",cex=1.0)
title("(d) Tratamento", cex.main=1.0)
legend(65,0.4, "logrank = 4,9, g.l.=1", cex=1.0, bty="n")
legend(65,0.31, "valor p = 0,03", cex=1.0, bty="n")10.4.2 - Ajuste do modelo de mistura semiparamétrico
x11<-ifelse(pele$Sexo=="M",1,0)
x21<-ifelse(pele$age=="40-59",1,0)
x22<-ifelse(pele$age=="60-80",1,0)
x31<-ifelse(Estadio==2,1,0)
x32<-ifelse(Estadio==3,1,0)
x33<-ifelse(Estadio==4,1,0)
x41<-ifelse(Trat=="Quimio",1,0)
dados<-as.data.frame(cbind(survtime,status,Idade,x11,x21,x22,x31,x32,x33,x41))
# iniciando com z = x = (sexo, faixa etária, estadiamento) #
library(smcure)
msm<-smcure(Surv(survtime,status) ~ x11 + x21 + x22 + x31 + x32 + x33,
cureform= ~ x11 + x21 + x22 + x31 + x32 + x33,
link="logit", data=dados, model="ph", Var=TRUE,
emmax=300, eps=1e-3, nboot=100)
msm<-smcure(Surv(survtime,status) ~ x11 + x31 + x32 + x33,
cureform= ~ x11 + x21 + x22 + x31 + x32 + x33,
link="logit", data=dados, model="ph", Var=TRUE,
emmax=300, eps=1e-3, nboot=100)
msm<-smcure(Surv(survtime,status) ~ x11 + x21 + x22 + x31 + x32 + x33,
cureform= ~ x11 + x31 + x32 + x33,
link="logit", data=dados, model="ph", Var=TRUE,
emmax=300, eps=1e-3, nboot=100)
msm<-smcure(Surv(survtime,status) ~ x11 + x31 + x32 + x33,
cureform= ~ x11 + x31 + x32 + x33,
link="logit", data=dados, model="ph", Var=TRUE,
emmax=300, eps=1e-3, nboot=100)
# Modelo final #
msm<-smcure(Surv(survtime,status) ~ x31 + x32 + x33,
cureform= ~ x11 + x31 + x32 + x33,
link="logit", data=dados, model="ph", Var=TRUE,
emmax=300, eps=1e-3, nboot=100)# Função de sobrevivência de base S0(t) => objeto s das saídas
# Se os tempos estiverem ordenados #
names(msm); par(mar=c(5,5,2,2)); par(mfrow=c(1,2))
plot(msm$Time, msm$s, type="s", xlab="Tempos", ylab="S0(t)")
# Se tempos não estiverem ordenados #
su0<-as.data.frame(cbind(msm$Time, msm$s)) # (tempos,S0(t))
names(su0)<-c("Time","S0")
i<-order(su0$Time)
sbase<-su0[i,]
plot(sbase$Time,sbase$S0, type="s",xlab="Tempos",ylab="S0(t)")Figura 10.4 - Qualidade de ajuste Sp(t|x,z)
# Combinação 1: feminino e estádio 0/1 #
comb1<-subset(pele, Sexo=="F" & Estadio==1)
ekm1<-survfit(Surv(survtime,status)~1, data=comb1)
pred1<-predictsmcure(msm,newX = cbind(c(0),c(0),c(0)), # estádio 0/1
newZ = cbind(c(0),c(0),c(0),c(0)), # feminino e estádio 0/1
model="ph")
# Combinação 5: masculino e estádio 0/1 #
comb5<-subset(pele,Sexo=="M" & Estadio==1)
ekm5<-survfit(Surv(survtime,status)~1, data=comb5)
pred5<-predictsmcure(msm,newX = cbind(c(0),c(0),c(0)), # estádio 0/1
newZ = cbind(c(1),c(0),c(0),c(0)), # masculino e estádio 0/1
model="ph")
# Correlação combinação 1 #
x1<-data.frame(cbind(pred1$prediction[,2],pred1$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm1$time, ekm1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Correlação combinação 5 #
x1<-data.frame(cbind(pred5$prediction[,2],pred5$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm5$time, ekm5$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
### Figura 10.4 ###
par(mfrow=c(1,2))
plot(c(0,ekm1$time),ylim=c(0,1), c(1,ekm1$surv), type="s", col=1,
lty=4, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(pred1$prediction[,2],pred1$prediction[,1], type="s", lty=1)
legend(70,0.5, lty=c(1,4), col=c(1, 1), c("Modelo","Kaplan-Meier"),bty="n")
text(100,0.1, expression(paste(rho)==0.863))
title("Feminino e estádio 0 ou 1", cex.main=0.9)
plot(c(0,ekm5$time),ylim=c(0,1), c(1,ekm5$surv), type="s", col=1,
lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(pred5$prediction[,2],pred5$prediction[,1], type="s", lty=1)
legend(70,0.5, lty=c(1,4), col=c(1, 1), c("Modelo","Kaplan-Meier"),bty="n")
text(100,0.1, expression(paste(rho)==0.951))
title("Masculino e estádio 0 ou 1", cex.main=0.9)Tabela 10.2 - Correlações
# Combinação 2: feminino e estádio 2 #
comb2<-subset(pele,Sexo=="F" & Estadio==2)
ekm2<-survfit(Surv(survtime,status)~1, data=comb2)
pred2<-predictsmcure(msm,newX = cbind(c(1),c(0),c(0)), # estádio 2
newZ = cbind(c(0),c(1),c(0),c(0)), # feminino e estádio 2
model="ph")
# Correlação #
x1<-data.frame(cbind(pred2$prediction[,2],pred2$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm2$time, ekm2$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Combinação 3: feminino e estádio 3 #
comb3<-subset(pele,Sexo=="F" & Estadio==3)
ekm3<-survfit(Surv(survtime,status)~1, data=comb3)
pred3<-predictsmcure(msm,newX = cbind(c(0),c(1),c(0)), # estádio 3
newZ = cbind(c(0),c(0),c(1),c(0)), # feminino e estádio 3
model="ph")
# Correlação #
x1<-data.frame(cbind(pred3$prediction[,2],pred3$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm3$time, ekm3$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Combinação 4: feminino e estádio 4 #
comb4<-subset(pele,Sexo=="F" & Estadio==4)
ekm4<-survfit(Surv(survtime,status)~1, data=comb4)
pred4<-predictsmcure(msm,newX = cbind(c(0),c(0),c(1)), # estádio 4
newZ = cbind(c(0),c(0),c(0),c(1)), # feminino e estádio 4
model="ph")
# Correlação #
x1<-data.frame(cbind(pred4$prediction[,2],pred4$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm4$time, ekm4$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Combinação 6: Masculino e estádio 2 #
comb6<-subset(pele,Sexo=="M" & Estadio==2)
ekm6<-survfit(Surv(survtime,status)~1, data=comb6)
pred6<-predictsmcure(msm,newX = cbind(c(1),c(0),c(0)), # estádio 2
newZ = cbind(c(1),c(1),c(0),c(0)), # masculino e estádio 2
model="ph")
# Correlação #
x1<-data.frame(cbind(pred6$prediction[,2],pred6$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm6$time, ekm6$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Combinação 7: masculino e estádio 3 #
comb7<-subset(pele,Sexo=="M" & Estadio==3)
ekm7<-survfit(Surv(survtime,status)~1, data=comb7)
pred7<-predictsmcure(msm,newX = cbind(c(0),c(1),c(0)), # estádio 3
newZ = cbind(c(1),c(0),c(1),c(0)), # masculino e estádio 3
model="ph")
# Correlação #
x1<-data.frame(cbind(pred7$prediction[,2],pred7$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm7$time, ekm7$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
# Combinação 8: masculino e estádio 4 #
comb8<-subset(pele,Sexo=="M" & Estadio==4)
ekm8<-survfit(Surv(survtime,status)~1, data=comb8)
pred8<-predictsmcure(msm,newX = cbind(c(0),c(0),c(1)), # estádio 4
newZ = cbind(c(1),c(0),c(0),c(1)), # masculino e estádio 4
model="ph")
# Correlação #
x1<-data.frame(cbind(pred8$prediction[,2],pred8$prediction[,1]))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm8$time, ekm8$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)Figura 10.5 - Qualidade de ajuste Su(t|x)
# x = Estádio 0/1 #
par(mfrow=c(2,2))
comb1<-subset(pele, Estadio==1 & status==1)
ekmc1<-survfit(Surv(survtime,status)~1, data=comb1)
plot(msm$Time, msm$s, type="s", xlab="Tempo (em meses)", ylab="S(t|x)")
lines(c(0,ekmc1$time),c(1,ekmc1$surv), type="s", lty=4, col=1)
text(95,0.6, expression(paste(rho)==0.997))
legend(70,0.9, lty=c(1,4), col=c(1,1), c("Modelo", "Kaplan-Meier"), bty="n")
title("Estádio 0 ou 1", cex.main=1.0)
# correlação #
x1<-data.frame(cbind(msm$Time, msm$s))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekmc1$time, ekmc1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
x1[x1$tempo==30,]
x1[x1$tempo==60,]
# x = Estádio 2 #
pred2<-(msm$s)^(exp(msm$beta[1]))
plot(msm$Time, pred2, type="s", xlab="Tempo (em meses)", ylab="S(t|x)" )
comb2<-subset(pele, Estadio==2 & status==1)
ekmc2<-survfit(Surv(survtime,status)~1, data=comb2)
lines(c(0,ekmc2$time),c(1,ekmc2$surv), type="s", lty=4, col=1)
legend(70,0.9, lty=c(1,4), col=c(1,1), c("Modelo", "Kaplan-Meier"), bty="n")
text(95,0.6, expression(paste(rho)==0.992))
title("Estádio 2", cex.main=1.0)
# correlação #
x1<-data.frame(cbind(msm$Time, pred2))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekmc2$time, ekmc2$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
x1[x1$tempo==30,]
x1[x1$tempo==60,]
# x= Estádio 3 #
pred3<-(msm$s)^(exp(msm$beta[2]))
plot(msm$Time, pred3, type="s", xlab="Tempo (em meses)", ylab="S(t|x)" )
comb3<-subset(pele, Estadio==3 & status==1)
ekmc3<-survfit(Surv(survtime,status)~1, data=comb3)
summary(ekmc3)
lines(c(0,ekmc3$time),c(1,ekmc3$surv), type="s", lty=4, col=1)
legend(70,0.9, lty=c(1,4), col=c(1,1), c("Modelo", "Kaplan-Meier"), bty="n")
text(95,0.6, expression(paste(rho)==0.996))
title("Estádio 3", cex.main=1.0)
# correlação #
x1<-data.frame(cbind(msm$Time, pred3))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekmc3$time, ekmc3$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
x1[x1$tempo==30,]
x1[x1$tempo==60,]
# x = Estádio 4 #
pred4<-(msm$s)^(exp(msm$beta[3]))
plot(msm$Time, pred4, type="s", ylim=c(0,1), xlab="Tempo (em meses)", ylab="S(t|x)")
comb4<-subset(pele, Estadio==4 & status==1)
ekmc4<-survfit(Surv(survtime,status)~1, data=comb4)
summary(ekmc4)
lines(c(0,ekmc4$time),c(1,ekmc4$surv),type="s", lty=4, col=1)
legend(70,0.9, lty=c(1,4), col=c(1, 1), c("Modelo", "Kaplan-Meier"), bty="n")
text(95,0.6, expression(paste(rho)==0.961))
title("Estádio 4", cex.main=0.9)
# correlação #
x1<-data.frame(cbind(msm$Time, pred4))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekmc4$time, ekmc4$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
x1[x1$tempo==30,]
x1[x1$tempo==60,]Figura 10.6 - Curvas preditas a partir do modelo de mistura
par(mar=c(5,5,2,2)); par(mfrow=c(1,2))
predm<-predictsmcure(msm, newX=cbind(c(1,1),c(0,0),c(0,0)), # Estádio 2
newZ=cbind(c(0,1),c(1,1),c(0,0),c(0,0)), # Fem Est2/Masc Est2
model="ph")
1-predm$newuncureprob
time<-predm$prediction[,3]
perfil1<-predm$prediction[,1]
perfil2<-predm$prediction[,2]
plot(c(0,time), c(1,perfil1), type="s", xlab="Tempo (em meses)",
ylab="Sp(t| z, x))", ylim=c(0,1))
lines(c(0,time), c(1,perfil2), type="s", lty=4)
text(113,0.79,"Estádio 2", cex=0.8)
legend(10,0.2,lty=c(1,4), c("Feminino","Masculino"), bty="n")
predm<-predictsmcure(msm, newX=cbind(c(0,0),c(0,0),c(1,1)), # Estádio 4
newZ=cbind(c(0,1),c(0,0),c(0,0),c(1,1)), # Fem Est4/Masc Est4
model="ph")
1-predm$newuncureprob
time<-predm$prediction[,3]
perfil1<-predm$prediction[,1]
perfil2<-predm$prediction[,2]
lines(c(0,time), c(1,perfil1), type="s", lty=1)
lines(c(0,time), c(1,perfil2), type="s", lty=4)
text(113,0.35,"Estádio 4", cex=0.8)
pred2<-(msm$s)^(exp(msm$beta[1]))
pred4<-(msm$s)^(exp(msm$beta[3]))
plot(c(0,msm$Time), c(1,pred2), type="s", xlab="Tempo (em meses)",
ylab="Su(t|x))", ylim=c(0,1))
lines(c(0,msm$Time), c(1,pred4), type="s", lty=4)
legend(65,0.8,lty=c(1,4), c("Estádio 2","Estádio 4"), bty="n")Seção 10.4.3 - Ajuste do modelo tempo de promoção
library(flexcure)
pele<-read.table("https://docs.ufpr.br/~giolo/asa/dados/pele.txt",h=T)
i<-order(pele$survtime)
pele<-pele[i,]
pele$age<-ifelse(pele$Idade<=39,"18-39",0)
pele$age<-ifelse(pele$Idade>39 & pele$Idade<60,"40-59",pele$age)
pele$age<-ifelse(pele$age==0,"60-80",pele$age)
table(pele$age)
attach(pele)
m1<-curereg(Surv(survtime, status) ~ factor(Estadio) + factor(Sexo) + factor(age),
cureformula = ~ factor(Estadio) + factor(Sexo) + factor(age),
data = pele, ncausedist = "poisson", timedist="weibull",
method="L-BFGS-B") # help(curereg) => ver métodos disponíveis
summary(m1)
m2<-curereg(Surv(survtime, status) ~ factor(Estadio),
cureformula = ~ factor(Sexo) + factor(Estadio),
data = pele, ncausedist = "poisson", timedist="weibull",
method="L-BFGS-B")
summary(m2)
coef(m2)
coef(m2, terms = "time")
coef(m2, terms = "cure")
AIC(m2) # Akaike information criterion
BIC(m2) # Bayesian information criterion# M ~ Poisson(theta(z)) e log(theta(z)) => theta(z) = exp(beta_0 + beta_1*z1 + ...)
# Proporção curados = exp(-theta(z)) = exp(-exp(exp(beta_0 + beta_1*z1 + ...))
exp(-exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[5])) # b0+b1*1+b4*1 (b1=1=masc, b4=1=est4)
exp(-exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[5])) # b0+b1*0+b4*1 (b1=0=fem, b4=1=est4)
cure <- curefraction(m2) # proporção de cura para todas as combinações
cure
mean(cure[,"curefraction"])
# Sp(t|x,z) para todas as combinações #
library(ggplot2)
par(mfrow=c(1,1))
plot(m2, mark.time=F, ylab="Sp(t)", xlab="Tempos")Figura 10.7 - Qualidade de ajuste do modelo
# Combinação: Masculino e Estádio 0 ou 1 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2]) # exp(b0+b1*1) (masculino e estádio 0/1)
b<- exp(m2$coefficients$coef.time[1]) # exp(g0) (estadio 0/1)
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a); Stx # Stx = exp(-(x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)); Spt # exp(-thetaz)*[1-S(t|x)])
comb1<-subset(pele,Sexo=="M" & Estadio==1)
ekm1<-survfit(Surv(survtime,status)~1, data=comb1)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm1$time, ekm1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
par(mfrow=c(1,2))
par(mar=c(5,5,2,2))
plot(c(0,ekm1$time),ylim=c(0,1), c(1,ekm1$surv), type="s", col=1,
lty=2, ylab="Sp(t|z,x)", xlab="Tempo (em meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1, 1), c("Modelo","Kaplan-Meier"), bty="n")
text(90,0.4, expression(paste(rho)==0.968))
title("Masculino e estádio 0 ou 1", cex.main=0.9)
# Combinação: Feminino e Estádio 0 ou 1 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1]) # exp(b0+b1*0): feminino e estádio 0/1
b<- exp(m2$coefficients$coef.time[1]) # exp(g0): estadio 0/1
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a); Stx # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)); Spt # exp(-thetaz)*[1 - S(t|x)])
comb2<-subset(pele,Sexo=="F" & Estadio==1)
ekm2<-survfit(Surv(survtime,status)~1, data=comb2)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm2$time, ekm2$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
par(mar=c(5,5,2,2))
plot(c(0,ekm2$time),ylim=c(0,1), c(1,ekm2$surv), type="s", col=1,
lty=2, ylab="Sp(t|z,x)", xlab="Tempo (em meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1, 1), c("Modelo","Kaplan-Meier"), bty="n")
text(99,0.4, expression(paste(rho)==0.943))
title("Feminino e estádio 0 ou 1", cex.main=0.9)Tabela 10.6 - Correlações
# Feminino e Estádio 2 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[3]) # feminino e estádio 2
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[2] ) # exp(g0+g1) estadio 2
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb<-subset(pele,Sexo=="F" & Estadio==2)
ekm<-survfit(Surv(survtime,status)~1, data=comb)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm$time, ekm$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
par(mar=c(5,5,2,2))
plot(c(0,ekm$time),ylim=c(0,1), c(1,ekm$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo", "Kaplan-Meier"), bty="n")
# Feminino e Estádio 3 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[4]) # feminino e estádio 3
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[3] ) # estádio 3
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb<-subset(pele,Sexo=="F" & Estadio==3)
ekm<-survfit(Surv(survtime,status)~1, data=comb)
summary(ekm)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm$time, ekm$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
plot(c(0,ekm$time),ylim=c(0,1), c(1,ekm$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo","Kaplan-Meier"), bty="n")
# Feminino e Estádio 4 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[5]) # feminino e estádio 4
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[4] ) # estádio 4
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb<-subset(pele,Sexo=="F" & Estadio==4)
ekm<-survfit(Surv(survtime,status)~1, data=comb)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm$time, ekm$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
plot(c(0,ekm$time),ylim=c(0,1), c(1,ekm$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo","Kaplan-Meier"), bty="n")
# Masculino e Estádio 2 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[3]) # masculino e estádio 2
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[2]) # estádio 2
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb1<-subset(pele,Sexo=="M" & Estadio==2)
ekm1<-survfit(Surv(survtime,status)~1, data=comb1)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm1$time, ekm1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
plot(c(0,ekm1$time),ylim=c(0,1), c(1,ekm1$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo","Kaplan-Meier"), bty="n")
# Masculino e Estádio 3
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[4]) # masculino e estádio 3
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[3]) # estádio 3
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb1<-subset(pele,Sexo=="M" & Estadio==3)
ekm1<-survfit(Surv(survtime,status)~1, data=comb1)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm1$time, ekm1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
plot(c(0,ekm1$time),ylim=c(0,1), c(1,ekm1$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo","Kaplan-Meier"), bty="n")
# Masculino e Estádio 4 #
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[5]) # masculino e estádio 4
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[4]) # estádio 4
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
comb1<-subset(pele,Sexo=="M" & Estadio==4)
ekm1<-survfit(Surv(survtime,status)~1, data=comb1)
# Correlação #
x1<-data.frame(cbind(t,Spt))
names(x1)<-c("tempo","surv")
x2<-data.frame(cbind(ekm1$time, ekm1$surv))
names(x2)<-c("tempo","surv")
xx<-merge(x1,x2, by="tempo")
cor(xx$surv.x,xx$surv.y)
plot(c(0,ekm1$time),ylim=c(0,1), c(1,ekm1$surv), type="s",
col=1, lty=2, ylab="Sp(t|z,x)", xlab="Tempo (meses)")
lines(t,Spt, type="s", lty=1)
legend(70,0.6, lty=c(1,4), col=c(1,1), c("Modelo","Kaplan-Meier"), bty="n")Figura 10.8 - Modelo de mistura vs. modelo tempo de promoção
### Feminino ###
t<-c(1:122)
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[3]) # feminino e estádio 2
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[2] ) # estádio 2
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
par(mar=c(5,5,2,2)); par(mfrow=c(1,2))
plot(t,Spt,ylim=c(0,1), type="s", col=1, lty=2,
ylab="Sp(t|z,x)", xlab="Tempo (meses)")
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[5]) # feminino e estádio 4
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[4] ) # estadio 4
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
lines(t, Spt, lty=2, type="l")
pred2<-predictsmcure(msm, newX = cbind(c(1),c(0),c(0)), # estádio 2
newZ = cbind(c(0),c(1),c(0),c(0)), # feminino e estádio 2
model="ph")
lines(pred2$prediction[,2],pred2$prediction[,1], type="s", lty=1)
pred4<-predictsmcure(msm,newX = cbind(c(0),c(0),c(1)), # estádio 4
newZ = cbind(c(0),c(0),c(0),c(1)), # feminino e estádio 4
model="ph")
lines(pred4$prediction[,2],pred4$prediction[,1], type="s", lty=1)
text(113,0.82,"Estádio 2", cex=0.8)
text(113,0.44,"Estádio 4", cex=0.8)
legend(10,0.2,lty=c(1,2), c("Modelo de mistura", "Modelo tempo de promoção"), bty="n")
title("Feminino", cex.main=1.1)
### Masculino ###
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[3]) # masculino e estádio 2
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[2]) # estádio 2
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
plot(t,Spt,ylim=c(0,1), type="s", col=1, lty=2,
ylab="Sp(t|z,x)", xlab="Tempo (meses)")
thetaz<- exp(m2$coefficients$coef.cure[1] +
m2$coefficients$coef.cure[2] +
m2$coefficients$coef.cure[5]) # masculino e estádio 4
b<- exp(m2$coefficients$coef.time[1] +
m2$coefficients$coef.time[4]) # estádio 4
a<- 1/exp(m2$coefficients$coef.time[5]) # gama = 1/sigma
Stx<- exp(-(t/b)^a) # Stx = exp(- (x/b)^a)
Spt<- exp((-thetaz)*(1 - Stx)) # exp(-thetaz)*[1 - S(t|x)])
lines(t, Spt, lty=2, type="l")
pred6<-predictsmcure(msm,newX = cbind(c(1),c(0),c(0)), # estádio 2
newZ = cbind(c(1),c(1),c(0),c(0)), # masculino e estádio 2
model="ph")
lines(pred6$prediction[,2],pred6$prediction[,1], type="s", lty=1)
pred8<-predictsmcure(msm,newX = cbind(c(0),c(0),c(1)), # estádio 4
newZ = cbind(c(1),c(0),c(0),c(1)), # masculino e estádio 4
model="ph")
lines(pred8$prediction[,2],pred8$prediction[,1], type="s", lty=1)
text(113,0.77,"Estádio 2", cex=0.79)
text(113,0.33,"Estádio 4", cex=0.8)
legend(10,0.2,lty=c(1,2), c("Modelo de mistura","Modelo tempo de promoção"), bty="n")
title("Masculino", cex.main=1.1)Capítulo 11
Análise de sobrevivência multivariada
Seção 11.6.1 - Estudo sobre Leucemia Pediátrica
library(survival)
leucc<-read.table("https://docs.ufpr.br/~giolo/asa/dados/leucc.txt",h=T)
attach(leucc)
id<-1:103
fit1<-coxph(Surv(tempos,cens)~leuinic+idadec+zpesoc+pasc+vacc+
frailty(id,dist="gamma"), data=leucc, method="breslow")
summary(fit1)
wi<-fit1$frail; zi<-exp(wi)
plot(id,zi, xlab="Crianças 1 a 103", ylab=expression(z[i]~estimados), pch=16)
abline(h=1,lty=2)Seção 11.6.2 - Estudo com animais da raça Nelore
cattle<-read.table("https://docs.ufpr.br/~giolo/asa/dados/cattle.txt",h=T)
attach(cattle)
library(survival)
fit1<-coxph(Surv(tempo,censura)~factor(sex) + agedam +
frailty(sire,dist="gamma"), data=cattle)
summary(fit1)
fit2<-coxph(Surv(tempo,censura)~factor(sex) +
frailty(sire,dist="gamma"), data=cattle)
summary(fit2)
## Figura 11.5
H0<-basehaz(fit2, centered=F)
S0<-exp(-H0$hazard)
S3m<-S0^(1.798*exp(0.797)) # machos touro 3
S3f<-S0^(1.798) # fêmeas touro 3
S1m<-S0^(0.767*exp(0.797)) # machos touro 1
S1f<-S0^(0.767) # fêmeas touro 1
par(mfrow=c(1,2))
t<-H0$time
plot(t,S1m, type="s", ylim=range(c(0,1)), xlab="Tempo (em dias)",
ylab="Sobrevivência estimada")
lines(t, S1f, type="s",lty=4)
legend(142, 0.25, lty=c(1,4),c("Machos","Fêmeas"), bty="n", cex=0.8)
title("Touro 1")
plot(t,S3m, type="s", ylim=range(c(0,1)), xlab="Tempo (em dias)",
ylab="Sobrevivência estimada")
lines(t, S3f, type="s", lty=4)
legend(142, 0.25, lty=c(1,4), c("Machos","Fêmeas"), bty="n", cex=0.8)
title("Touro 3")Seção 11.6.3 - Estudo sobre infecções recorrentes
## Modelo AG (escala TC)
library(survival)
data(cgd)
fitag <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age +
inherit + steroids + cluster(id), data=cgd); fitag
fitag <- coxph(Surv(tstart, tstop, status) ~ treat + age + inherit +
steroids + cluster(id), data=cgd); fitag
fitag <- coxph(Surv(tstart, tstop, status) ~ treat + age + steroids +
cluster(id), data=cgd); fitag
# Análise dos resíduos
cox.zph(fitag, transform="identity")
resm<-resid(fitag,type="martingale")
res<-cgd$status - resm # resíduos de Cox-Snell
ekm <- survfit(Surv(res, cgd$status)~1)
par(mfrow=c(1,2))
plot(ekm, mark.time=F, conf.int=F, xlab="Resíduos de Cox-Snell",
ylab="S(e) estimada")
res<-sort(res); exp1<-exp(-res)
lines(res, exp1, lty=3)
legend(0.5, 0.9, lty=c(1,3), c("Kaplan Meier","Exponencial padrão"),
lwd=1, bty="n", cex=0.9)
st<-ekm$surv; t<-ekm$time; sexp1<-exp(-t)
plot(st,sexp1, ylim=c(0,1), xlim=c(0,1),
xlab="S(e): Kaplan-Meier", ylab="S(e): exponencial padrão",
pch=16, cex=1.0)
abline(0,1)Comandos para ajuste dos modelos PWP e WLW
## Modelo PWP (escala TC) ##
## Considerando beta comum ao estratos ##
library(survival)
data(cgd)
pwpfit<- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit +
steroids + cluster(id) + strata(enum), data = cgd)
summary(pwpfit)
## limitando a N = 2
cgd1<-subset(cgd, enum<=2); dim(cgd1) # variável enum = estratos
pwpfit1<- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit +
steroids + cluster(id) + strata(enum), data = cgd1)
summary(pwpfit1)
## Considerando beta_k e N = 2 ##
pwpfit2<- coxph(Surv(tstart, tstop, status) ~ treat*strata(enum) + age*strata(enum) +
steroids*strata(enum) + cluster(id), data = cgd1)
summary(pwpfit2)## Modelo WLW (escala TT)
## Nota: Necessário criar o arquivo cgd2 => dados no layout da Tabela 11.2
## Considerando beta comum ##
wlwfit<-coxph(Surv(time, status) ~ treat + age + steroids +
cluster(id), data=cgd2)
## Considerando beta_k ##
wlwfit<-coxph(Surv(time, status) ~ treat*strata(enum) + age*strata(enumm) +
steroids*strata(enum) + cluster(id), data=cgd2)
## Obs: Para limitar a N = 2, criar cgd3 => subset de cgd2
cgd3<-subset(cgd2, enum<=2)Copyright © 2024 Suely Ruiz Giolo