require(survival) source("Rfun.r") dados <- read.csv('data.csv', header=T, dec=",") dados$ano <- as.factor(dados$ano) dados$arvore <- as.factor(dados$arvore) (s1 <- survfit(Surv(tempo, status)~ cult, data=subset(dados, ano=="1"))) (s2 <- survfit(Surv(tempo, status)~ cult, data=subset(dados, ano=="2"))) (s <- survfit(Surv(tempo, status)~ cult, data=dados)) survdiff(Surv(tempo, status)~ cult + strata(ano), data=subset(ensac, ano==1)) survdiff(Surv(tempo, status)~ cult + strata(ano), data=subset(ensac, ano==2)) survdiff(Surv(tempo, status)~ cult + strata(ano), data=ensac) fit <- coxph(Surv(tempo, status) ~ strata(ano)+ cult + frailty(arvore), data=dados) summary(fit) (residuo.sch <- cox.zph(fit)) ### Plots par(mfrow=c(2,2),mar=c(5,4,1,1), oma=c(1,1,1,1)) # kaplan meier par(mar=c(4,4,1,1)) plot(s, conf.int=F, lty=c(2:4), ylim=c(0.5, 1), xlab="Survival time (days)", font.lab=2, font=2, cex=1.5, ylab="Probability of survival", main="", lwd=c(2, 3, 4)) legend("bottomright", bty="n", cex=2, "a") legend(0.5, 0.7, c('Cultivar A','Cultivar B',' Cultivar C'), cex=1.2, bty="n",x.intersp=0.2, y.intersp=1, xjust=0, lty=2:4, lwd=c(2, 3, 4)) # frailty par(mar=c(4,4,1,1)) plot.frail(ensac$arvore, fit, xlab="Peach tree", font=2, font.lab=2, ylim=c(0.3,2.3), main="") legend("bottomright", bty="n", cex=2, "b") legend("topleft", bty="n", cex=1.2, "Frailty (p = 0,17)") # residuals par(mar=c(4,4,1,1)) plot(residuo.sch, var="cultB",font=2, font.lab=2, main="", ylim=c(-21,21)) legend("bottomright", bty="n", cex=2, "c") legend("topleft", bty="n", cex=1.2, "r = 0,14 (p = 0,0925)") par(mar=c(4,4,1,1)) plot(residuo.sch, var="cultC",font=2, font.lab=2, main="", ylim=c(-21,21)) legend("bottomright", bty="n", cex=2, "d") legend("topleft", bty="n", cex=1.2, "r = 0,12 (p = 0,1442)") dev.off()