Códigos em Linguagem R Apresentados no Texto
Códigos em Linguagem R Apresentados no Texto
Capítulo 2 - Técnicas não-paramétricas
- Seção 2.3 - Estimador de Kaplan-Meier
require(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))
ekm<- survfit(Surv(tempos,cens)~grupos)
summary(ekm)
plot(ekm, lty=c(2,1), xlab="Tempo (semanas)",ylab="S(t) estimada")
legend(1, 0.3, lty=c(2,1), c("Controle","Esteróide"),lwd=1, bty="n")
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)
summary(ekm)- Seção 2.4.1 - Estimador de Nelson-Aalen
require(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))
ss<- survfit(coxph(Surv(tempos[grupos==2],cens[grupos==2])~1, method = "breslow"))
summary(ss)
racum<- -log(ss$surv)
cbind(ss$time, racum)- Seção 2.5.1 - Exemplo: Reincidência de Tumor Sólido
require(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", bty="n")
t<- tempos[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ekm$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
tm<-sum(prod)
tm
Função para obtenção da variância do tempo médio. Autora: Vanessa F. Sehaber. Clique para ver a função var_tm
** Utilizando a função var_tm aos dados de reincidência de tumor sólido
source("https://docs.ufpr.br/~giolo/Livro/var_tm.txt") # carregando var_tm
ekm<- survfit(Surv(tempos,cens)~1, conf.type="plain")
var.tmedio <- var_tm(ekm)
var.tmedio- Seção 2.6 - Comparação de Curvas de Sobrevivência
require(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 - Análise dos Dados da Malária
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))
require(survival)
ekm<-survfit(Surv(tempos,cens)~grupo)
summary(ekm)
plot(ekm, lty=c(1,4,2), xlab="Tempo",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)
survdiff(Surv(tempos[1:31],cens[1:31])~grupo[1:31], rho=0)
survdiff(Surv(tempos[17:44],cens[17:44])~grupo[17:44], rho=0)
survdiff(Surv(c(tempos[1:16],tempos[32:44]),c(cens[1:16],cens[32:44]))~c(grupo[1:16],grupo[32:44]), rho=0)
# pode-se, alternativamente, usar a opção "subset"
data<-as.data.frame(cbind(tempos, cens, grupo))
dat1<-subset(data, subset = grupo %in% c(1,2)) # grupos 1 e 2
survdiff(Surv(tempos,cens) ~ grupo, rho=0, dat1)
dat2<-subset(data, subset = grupo %in% c(2,3)) # grupos 2 e 3
survdiff(Surv(tempos,cens) ~ grupo, rho=0, dat2)
dat3<-subset(data, subset = grupo %in% c(1,3)) # grupos 1 e 3
survdiff(Surv(tempos,cens) ~ grupo, rho=0, dat3)Capítulo 3 - Modelos Probabilísticos
- Seção 3.6.1 - Exemplo 1: Pacientes com Câncer de Bexiga
require(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=range(c(0.0,1)), xlim=range(c(0,1)), xlab = "S(t): Kaplan-Meier",
ylab="S(t): exponencial")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(st, stw, pch=16, ylim=range(c(0.0,1)), xlim=range(c(0,1)), xlab = "S(t): Kaplan-Meier",
ylab="S(t): Weibull")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(st, stln, pch=16, ylim=range(c(0.0,1)), xlim=range(c(0,1)), xlab = "S(t): Kaplan-Meier",
ylab="S(t): log-normal")
lines(c(0,1), c(0,1), type="l", lty=1)
par(mfrow=c(1,3))
invst<-qnorm(st)
plot(time, -log(st), pch=16, xlab="tempos", ylab="-log(S(t))")
plot(log(time), log(-log(st)), pch=16, xlab="log(tempos)", ylab="log(-log(S(t)))")
plot(log(time), invst, pch=16, xlab="log(tempos)", ylab=expression(Phi^-1*(S(t))))
ajust1$loglik[2]
ajust2$loglik[2]
ajust3$loglik[2]
par(mfrow=c(1,2))
plot(ekm, conf.int=F, xlab="Tempos", 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.8)
plot(ekm, conf.int=F, xlab="Tempos", 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.8) Ajuste do modelo gama generalizado no R, via o pacote “flexsurv”. Veja a seguir.
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)
require(flexsurv) # instalar previamente o pacote
ajust4<-flexsurvreg(Surv(tempos,cens)~1, dist='gengamma')
ajust4$loglik
ajust4
summary(ajust4)
ekm<-survfit(Surv(tempos,cens)~1)
time<-ekm$time
st<-ekm$surv
gg<-as.data.frame(summary(ajust4))
stgg<-gg$est
par(mfrow=c(1,2))
plot(st, stgg, pch=16, ylim=range(c(0.0,1)), xlim=range(c(0,1)), xlab = "S(t): Kaplan-Meier",
ylab="S(t): gama generalizada")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(ekm, conf.int=F, xlab="Tempos", ylab="S(t)")
lines(c(0,gg$time), c(1,stgg), lty=2)
legend(25, 0.8, lty=c(1,2), c("Kaplan-Meier", "Gama generalizada"), bty="n", cex=0.8) Capítulo 4 - Modelos de Regressão Paramétricos
- Seção 4.5.1 - Sobrevida de Pacientes com Leucemia Aguda
temp<-c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
cens<-rep(1,17)
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,lwbc)
require(survival)
dados<-as.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)
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))))
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- Seção 4.5.3 - Análise dos Dados de Aleitamento Materno
rm(list=ls(all=TRUE))
desmame<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/desmame.txt",h=T)
attach(desmame)
require(survival)
ekm<- survfit(Surv(tempo,cens)~V4)
summary(ekm)
survdiff(Surv(tempo,cens)~V4, rho=0)
plot(ekm, lty=c(1,4), mark.time=F, xlab="Tempo até o desmame (meses)", ylab="S(t)")
text(18.5, 0.93, c("Dificuldades para Amamentar"), bty="n", cex=0.85)
legend(15.5, 0.9, lty=c(4), c("Sim"), bty="n", cex=0.8)
legend(18.5, 0.9, lty=c(1), c("Não"), bty="n", cex=0.8)
ajust1<-survreg(Surv(tempo,cens)~V1+V3+V4+V6, dist='lognorm')
ajust1
summary(ajust1)
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")
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 - Análise de um Estudo sobre Câncer de Laringe
laringe<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/laringe.txt", h=T)
attach(laringe)
require(survival)
fit2<-coxph(Surv(tempos,cens)~factor(estagio), data=laringe, x = T, method="breslow")
summary(fit2)
fit2$loglik
fit3<- coxph(Surv(tempos,cens)~factor(estagio) + idade, data=laringe, x = T, method="breslow")
summary(fit3)
fit3$loglik
fit4<-coxph(Surv(tempos,cens) ~ factor(estagio) + idade + factor(estagio)*idade, data=laringe,
x = T, method="breslow")
summary(fit4)
fit4$loglik
resid(fit4, type="scaledsch")
cox.zph(fit4, transform="identity") # g(t) = t
par(mfrow=c(2,4))
plot(cox.zph(fit4))
resid(fit3, type="scaledsch")
cox.zph(fit3, transform="identity") # g(t) = t
par(mfrow=c(1,4))
plot(cox.zph(fit3))
Ht<-basehaz(fit4, centered=F)
tempos<-Ht$time
H0<-Ht$hazard
S0<- exp(-H0)
round(cbind(tempos, S0, H0), digits=5) Análise gráfica dos resíduos de Cox-Snell do modelo ajustado aos dados de laringe.
resm<-resid(fit4,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", 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(1)"), lwd=1, bty="n", cex=0.7)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): Exponencial(1)", pch=16)- Seção 5.7.2 - Análise dos Dados de Aleitamento Materno
require(survival)
desmame<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/desmame.txt",h=T)
attach(desmame)
fit<-coxph(Surv(tempo,cens)~V1+V3+V4+V6, data=desmame, x = T, method="breslow")
summary(fit)
fit$loglik
resid(fit, type="scaledsch")
cox.zph(fit, transform="identity") # g(t) = t
par(mfrow=c(2,2))
plot(cox.zph(fit)) Análise gráfica dos resíduos de Cox-Snell do modelo ajustado aos dados de aleitamento materno.
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", 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(1)"), lwd=1, bty="n", cex=0.7)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): Exponencial(1)", pch=16)- Seção 5.7.3 - Análise dos Dados de Leucemia Pediátrica
rm(list=ls())
leuc<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/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)
require(survival)
fit<-coxph(Surv(tempos,cens)~leuinic + idadec + zpesoc + zestc + pasc + vacc + riskc + r6c, data=leucc,
x = T, method="breslow")
summary(fit)
fit3<-coxph(Surv(tempos,cens)~leuinic + idadec + zpesoc + pasc + vacc, data=leucc, x = T, method="breslow")
summary(fit3)
-2*fit3$loglik[2]
resid(fit3, type="scaledsch")
cox.zph(fit3, transform="identity") # g(t) = t
par(mfrow=c(2,3))
plot(cox.zph(fit3))
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íduo martingal", pch=16)
plot(pl, rd, xlab="Preditor linear", ylab="Resíduo deviance" , pch=16)
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") Análise gráfica dos resíduos de Cox-Snell do modelo ajustado aos dados de leucemia.
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(1)"), lwd=1, bty="n", cex=0.7)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): Exponencial(1)", pch=16)Capítulo 6 - Extensões do Modelo de Cox
- Seção 6.4.1 - Análise dos Dados de Pacientes HIV
aids<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/aids.txt", h=T)
attach(aids)
require(survival)
fit1<-coxph(Surv(ti[ti<tf], tf[ti<tf], cens[ti<tf])~id[ti<tf] + factor(grp)[ti<tf], method="breslow")
summary(fit1) NOTA: Alternativamente, pode-se usar subset() aids1<-subset(aids, ti < tf)
fit2<-coxph(Surv(ti, tf, cens)~id + factor(grp), method="breslow", data=aids1)
summary(fit2) Análise gráfica dos resíduos de Cox-Snell do modelo ajustado aos dados de HIV.
aids2<-na.omit(aids1[, 1:7])
fit2<- coxph(Surv(ti, tf, cens)~id + factor(grp), method="breslow", data=aids2)
resm<-resid(fit2, 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(1, 0.8, lty=c(1,3), c("Kaplan-Meier","Exponencial(1)"), lwd=1, bty="n", cex=0.7)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): Exponencial(1)", pch=16)- Seção 6.5 - Modelo de Cox Estratificado ajustado aos Dados de Leucemia
leuc<-read.table("http://www.ufpr.br/~giolo/Livro/ApendiceA/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)
require(survival)
fit1<-coxph(Surv(tempos,cens)~idadec + zpesoc + pasc + vacc + strata(leuinic), data=leucc,
x = T,method="breslow")
summary(fit1)
leucc1<-as.data.frame(cbind(tempos[leuinic==0],cens[leuinic==0],idadec[leuinic==0],
zpesoc[leuinic==0],pasc[leuinic==0],vacc[leuinic==0]))
leucc2<-as.data.frame(cbind(tempos[leuinic==1],cens[leuinic==1],idadec[leuinic==1],
zpesoc[leuinic==1],pasc[leuinic==1],vacc[leuinic==1]))
fit2<-coxph(Surv(V1,V2)~V3+V4+V5+V6, data=leucc1, x = T, method="breslow")
summary(fit2)
fit3<-coxph(Surv(V1,V2)~V3+V4+V5+V6, data=leucc2, x = T, method="breslow")
summary(fit3)
trv<-2*(-fit1$loglik[2] + fit2$loglik[2] + fit3$loglik[2])
trv
1-pchisq(trv,4)
cox.zph(fit1, transform="identity") # g(t) = t
par(mfrow=c(1,4))
plot(cox.zph(fit1))
H0<-basehaz(fit1, centered=F) # risco acumulado de base
H0
H01<-as.matrix(H0[1:81,1]) # risco acumulado de base do estrato 1 (leuinic=0)
H02<-as.matrix(H0[82:101,1]) # risco acumulado de base do estrato 2 (leuinic=1)
tempo1<- H0$time[1:81] # tempos do 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 do 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, H02, lty=4, type="s", xlab="Tempos", xlim=range(c(0,4)), ylab=expression(Lambda[0]* (t)))
lines(tempo1, H01, type="s", lty=1)
legend(0.0, 9, 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="Tempos", xlim=range(c(0,4)), ylab="So(t)")
lines(c(0, tempo2), c(1,S02), lty=4, type="s")
legend(2.2, 0.9, lty=c(1,4), c("Leuini < 75000","Leuini > 75000"), lwd=1,bty="n", cex=0.8)- Seção 6.5 - Estudo sobre Hormônio de Crescimento
hg2<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/hg2.txt", h=T)
attach(hg2)
require(survival)
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")
par(mfrow=c(2,3))
plot(cox.zph(fit2)) Análise gráfica dos resíduos de Cox-Snell do modelo ajustado aos dados de hormônio de crescimento.
hg2a<-na.omit(hg2)
rendac<-ifelse(hg2a$renda < 4, 1, 2)
fit2<-coxph(Surv(tempos,cens)~factor(raca) + factor(trauma) + factor(recemnas) + factor(rendac) + ialtura +
factor(trauma)*factor(recemnas), data=hg2a, method="breslow")
resm<-resid(fit2,type="martingale")
res<-hg2a$cens - resm # resíduos de Cox-Snell
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", 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(1)"), lwd=1, bty="n", cex=0.7)
st<-ekm$surv
t<-ekm$time
sexp1<-exp(-t)
plot(st, sexp1, xlab="S(e): Kaplan-Meier", ylab= "S(e): Exponencial(1)", pch=16)Capítulo 7 - Modelo Aditivo de Aalen
- Seção 7.6 - Análise dos Dados de Câncer de Laringe
laringe<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/laringe.txt", h=T)
attach(laringe)
require(survival)
source("https://docs.ufpr.br/~giolo/Livro/ApendiceA/Addreg.r")
idadec<-idade - mean(idade)
fit1<- addreg(Surv(tempos,cens)~factor(estagio) + idadec, laringe)
summary(fit1)
fit2<- addreg(Surv(tempos,cens)~factor(estagio), laringe)
summary(fit2)
i<-order(tempos)
laringe<-laringe[i,] # dados ordenados pelos tempos
laringe1<-laringe[1:51,] # obs: como tau = 4.3 deve-se usar somente as linhas em que t <= 4.3
xo<-rep(1,51)
x1<-ifelse(laringe1$estagio==2,1,0)
x2<-ifelse(laringe1$estagio==3,1,0)
x3<-ifelse(laringe1$estagio==4,1,0)
x <-as.matrix(cbind(xo,x1,x2,x3))
t<-fit2$times
coef<-fit2$increments
xt<-t(x)
Bt<-coef%*%xt
riscoacum<-diag(Bt)
for(i in 1:50){
riscoacum[i+1]<-riscoacum[i+1]+riscoacum[i]
}
riscoacum
par(mfrow=c(1,1))
plot(t,riscoacum,xlab="Tempos", ylab = expression(Lambda*(t)), pch=16)
plot(fit2,xlab="Tempo",ylab="FRA",labelofvariable=c("Estadio I","Estadio II em relação ao I",
"Estadio III em relação ao I","Estadio IV em relação ao I"))- Seção 7.7 - Análise dos Dados de Pacientes HIV
rm(list=ls())
source("https://docs.ufpr.br/~giolo/Livro/ApendiceA/Addreg.r") # lendo função Addreg.r
aids<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/aids.txt", h=T) # lendo os dados
attach(aids)
require(survival)
idade<-id - mean(id[!is.na(id)])
fit1<-addreg( Surv(ti[ti<tf],tf[ti<tf],cens[ti<tf])~idade[ti<tf]+sex[ti<tf]+factor(grp)[ti<tf], data=aids)
summary(fit1)
fit2<-addreg( Surv(ti[ti<tf],tf[ti<tf],cens[ti<tf])~idade[ti<tf]+factor(grp)[ti<tf], data=aids)
summary(fit2)
aids<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/aids.txt", h=T)
attach(aids)
aids1<-as.data.frame(cbind(tf,id,grp))
aids1<-na.omit(aids1) # eliminando valores missing = NA
attach(aids1)
i<-order(aids1$tf)
aids1<-aids1[i,] # dados ordenados por tf e sem NA nas covariaveis
aids2<-aids1[10:121,] # tau = 617 -> foram mantidas as linhas em que 0 < tf <= 617
n<-nrow(aids2)
xo<-rep(1,n)
x1<-(aids2$id) - mean(aids2$id)
x2<-ifelse(aids2$grp==2,1,0)
x3<-ifelse(aids2$grp==3,1,0)
x4<-ifelse(aids2$grp==4,1,0)
x <-as.matrix(cbind(xo,x1,x2,x3,x4))
t<-fit2$times
coef<-fit2$increments
xt<-t(x)
Bt<-coef%*%xt
riscoacum<-diag(Bt)
for(i in 1:(n-1)){
riscoacum[i+1]<-riscoacum[i+1]+riscoacum[i]
}
riscoacum
plot(t, riscoacum, xlab="Tempos", ylab = expression(Lambda*(t)), pch=16)
plot(fit2, xlab="Tempo", ylab="FRA",label=c("(a) Constante","(b) Idade","(c) Assintomático em relação HIV-",
"(d) ARC em relação HIV-", "(e) AIDS em relação HIV-"))Capítulo 8 - Censura Intervalar e Dados Grupados
- Seção 8.2 - Técnicas não-paramétricas
require(survival)
source("http://docs.ufpr.br/~giolo/Livro/ApendiceE/Turnbull.R") # lendo a função 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)
tb- Seção 8.2.1 - Exemplo de Câncer de Mama
rm(list=ls())
require(survival)
source("https://docs.ufpr.br/~giolo/Livro/ApendiceE/Turnbull.R") # lendo Turnbull.R
dat <- read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/breast.txt", h=T) # lendo os dados
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
plot(tb1$time, tb1$surv, lty=1, type = "s", ylim=c(0,1), xlim=c(0,50), xlab="Tempos (meses)", ylab="S(t)")
lines(tb2$time, tb2$surv, lty=4, type="s")
legend(1, 0.3, lty=c(1,4), c("Radioterapia","Radioterapia + Quimioterapia"), bty="n", cex=0.8)
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=1, type="s", ylim=c(0,1), xlim=c(0,50), xlab="Tempos (meses)", ylab="S(t)")
lines(tb2$time, tb2$surv, lty=2, type="s")
lines(ekm[1]$time, ekm[1]$surv, type="s", lty=3)
lines(ekm[2]$time, ekm[2]$surv, type="s", lty=3)
legend(1, 0.3, lty=c(1,2), c("Radioterapia","Radioterapia + Quimioterapia"), bty="n", cex=0.8)
legend(1, 0.21, lty=3, "Usando Ponto Médio dos Intervalos", bty="n", cex=0.8)- Seção 8.3.1 - Análise dos Dados de Câncer de Mama (Modelo paramétrico)
breast<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/breast.txt", h=T)
attach(breast)
cens1<-ifelse(cens==1,3,0)
require(survival)
fit1<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="logistic")
summary(fit1)
fit2<-survreg(Surv(left,right,type="interval2")~ther, breast, dist="gaussian")
summary(fit2)
t1<-0:50
b0<-fit1$coefficients[1]
b1<-fit1$coefficients[2]
s<- fit1$scale
a1<- t1-(b0+b1)
e1<- exp(a1/s)
st1<-1/(1+e1)
t2<-0:50
a2<- t2-(b0)
e2<- exp(a2/s)
st2<-1/(1+e2)
plot(t1, st1, type="l", lty=3, ylim=range(c(0,1)), xlab="Tempos", ylab="S(t) estimada")
lines(t2, st2, type="l", lty=3)
t1<-0:45
b0<-fit2$coefficients[1]
b1<-fit2$coefficients[2]
s<- fit2$scale
a1<- t1-(b0+b1)
st11<- 1-pnorm(a1/s)
t2<-0:60
a2<-t2-(b0)
st22<- 1 -pnorm(a2/s)
lines(t2, st22, type="l", lty=2)
lines(t1, st11, type="l", lty=2)
legend(1,0.2, lty=c(3,2), c("Logística","Gaussiana"), lwd=1, bty="n", cex=0.8) A análise a seguir usa o pacote “icenReg”, tendo em vista o pacote “intcox” não estar mais disponível no R.
- Seção 8.4.1 - Análise dos Dados de Câncer de Mama via o pacote icenReg
breast<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/breast.txt", h=T)
attach(breast)
require(icenReg) # instalar previamente o pacote
fit1<- ic_sp(cbind(left, right) ~ ther, model = 'ph', bs_samples = 100, data = breast)
summary(fit1)- Obs: as estimativas não são idênticas às do intcox => procedimentos de estimação diferentes são considerados nesses pacotes.
- Para visualizar os comandos os quais eram utilizados no intcox, clique no link a seguir -> cmd_intcox
- Seção 8.8 - Ensaio de Vida de Mangueiras (Dados Grupados)
mang<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/mang.txt", h=T)
attach(mang)
require(survival)
ekm<-survfit(Surv(ti,cens)~1, conf.type="none")
summary(ekm)
mang1<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/dadmang.txt", h=T) # ver Apêndice A5
attach(mang1)
require(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)
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)
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)Capítulo 9 - Análise de Sobrevivência Multivariada
- Seção 9.11.1 - Fragilidade no Estudo de Leucemina Pediátrica
leucc<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/leucc.txt", h=T) # dados dicotomizados
attach(leucc)
require(survival)
id<-1:103
fit3a<-coxph(Surv(tempos,cens)~leuinic+idadec+zpesoc+pasc+vacc+frailty(id,dist="gamma"), data=leucc,
x = T, method="breslow")
summary(fit3a)
wi<-fit3a$frail
zi<-exp(wi)
plot(id,zi, xlab="Crianças (1 a 103)", ylab="zi estimados", pch=16)
abline(h=1,lty=2)- Seção 9.11.2 - Estudo com Animais da Raça Nelore
cattle<-read.table("https://docs.ufpr.br/~giolo/Livro/ApendiceA/cattle.txt",h=T)
attach(cattle)
require(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)
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 (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 (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")- Revisado em 15/02/2018.