Curso R - UFMG - Resolução do exercício de sobrevivência
Ricardo Solar - UFMG (ricardosolar@ufmg.br)
December, 2023
© Ricardo Solar/UFMG - compartilhamento e utilização não-comercial livres. Não reproduzir sem autorização
1 Procedimentos iniciais
rm(list=ls()) ##remove todos os elementos existentes no Environment
dev.off() ##apaga todo o device grafico## null device
## 11.1 Carregando os pacotes que usaremos
library(survival)
library(ggplot2)
library(survminer)## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myeloma1.2 Carregando os dados
1.3 Utilizaremos o arquivo ataque.txt
dados<- read.table("ataque.txt",header=T, stringsAsFactors = T)
#attach(dados)
summary(dados)## Interacao Distancia tempo Ataque
## hom:62 longe:61 Min. : 5.0 Min. :0.0000
## nef:60 perto:61 1st Qu.:166.2 1st Qu.:0.0000
## Median :600.0 Median :0.0000
## Mean :406.1 Mean :0.4836
## 3rd Qu.:600.0 3rd Qu.:1.0000
## Max. :600.0 Max. :1.0000plot(dados)
1.4 Nesta tabela de dados temos:
- tempo: eh o tempo para ataque que cada individuo (cupim) sobre por formigas
- ataque: se o individuo foi atacado (1) ou nao (0)
- distancia: perto ou longe (ja pensou que legal se essa distancia fosse contiua?)
- Tipo de interacao ecologica: NEF ou hom
1.4.1 Pergunta: Tempo de ataque e dependente do tipo de interacao ecologica e da distancia?
2 Resolvendo o exercicio
m.completo <- survreg(Surv(tempo,Ataque)~Interacao*Distancia, data=dados)
m.n <- survreg(Surv(tempo,Ataque)~1, data=dados)
anova(m.n,m.completo,test="Chisq")anova(m.completo,test="Chisq")# simplificando:
#retirando a interacao:
m.completo2 <- update(m.completo, .~.-Interacao:Distancia)
anova(m.completo, m.completo2, test="Chi")anova(m.completo2,test="Chisq")summary(m.completo2)##
## Call:
## survreg(formula = Surv(tempo, Ataque) ~ Interacao + Distancia,
## data = dados)
## Value Std. Error z p
## (Intercept) 7.411 0.362 20.47 < 2e-16
## Interacaonef 0.816 0.364 2.25 0.025
## Distanciaperto -1.632 0.404 -4.04 5.4e-05
## Log(scale) 0.294 0.116 2.55 0.011
##
## Scale= 1.34
##
## Weibull distribution
## Loglik(model)= -439 Loglik(intercept only)= -451.1
## Chisq= 24.33 on 2 degrees of freedom, p= 5.2e-06
## Number of Newton-Raphson Iterations: 5
## n= 122- o tempo de ataque qd tem NEF e superior (p=0.049)
- quanto mais perto mais rapido e encontrado nos dois tratamentos (p<0.001)
3 Fazendo o grafico, como no exemplo anterior
3.1 Vamos fazer uma variavel por vez, pra entender o que esta acontecendo
3.1.1 Primeiro o tipo de interacao ecologica
m.completo2.1 = survfit(Surv(tempo,Ataque)~Interacao, data=dados)
summary(m.completo2.1)## Call: survfit(formula = Surv(tempo, Ataque) ~ Interacao, data = dados)
##
## Interacao=hom
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 62 1 0.984 0.0160 0.953 1.000
## 10 61 3 0.935 0.0312 0.876 0.999
## 15 58 1 0.919 0.0346 0.854 0.990
## 20 57 1 0.903 0.0375 0.833 0.980
## 33 56 1 0.887 0.0402 0.812 0.969
## 35 55 1 0.871 0.0426 0.791 0.959
## 40 54 2 0.839 0.0467 0.752 0.935
## 46 52 1 0.823 0.0485 0.733 0.923
## 56 51 1 0.806 0.0502 0.714 0.911
## 70 50 1 0.790 0.0517 0.695 0.898
## 82 49 1 0.774 0.0531 0.677 0.886
## 95 48 1 0.758 0.0544 0.659 0.873
## 102 47 1 0.742 0.0556 0.641 0.859
## 112 46 1 0.726 0.0567 0.623 0.846
## 115 45 1 0.710 0.0576 0.605 0.832
## 132 44 1 0.694 0.0585 0.588 0.818
## 145 43 1 0.677 0.0594 0.571 0.804
## 170 42 2 0.645 0.0608 0.536 0.776
## 202 40 1 0.629 0.0613 0.520 0.762
## 206 39 1 0.613 0.0619 0.503 0.747
## 265 38 1 0.597 0.0623 0.486 0.732
## 266 37 1 0.581 0.0627 0.470 0.717
## 275 36 1 0.565 0.0630 0.454 0.702
## 311 35 1 0.548 0.0632 0.438 0.687
## 353 34 1 0.532 0.0634 0.421 0.672
## 470 33 1 0.516 0.0635 0.406 0.657
## 510 32 2 0.484 0.0635 0.374 0.626
## 537 30 1 0.468 0.0634 0.359 0.610
## 540 29 1 0.452 0.0632 0.343 0.594
## 545 28 1 0.435 0.0630 0.328 0.578
##
## Interacao=nef
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 10 60 2 0.967 0.0232 0.922 1.000
## 30 58 2 0.933 0.0322 0.872 0.999
## 40 56 1 0.917 0.0357 0.849 0.989
## 52 55 1 0.900 0.0387 0.827 0.979
## 60 54 1 0.883 0.0414 0.806 0.968
## 67 53 1 0.867 0.0439 0.785 0.957
## 120 52 1 0.850 0.0461 0.764 0.945
## 125 51 1 0.833 0.0481 0.744 0.933
## 165 50 1 0.817 0.0500 0.724 0.921
## 193 49 1 0.800 0.0516 0.705 0.908
## 200 48 1 0.783 0.0532 0.686 0.895
## 231 47 1 0.767 0.0546 0.667 0.882
## 245 46 1 0.750 0.0559 0.648 0.868
## 300 45 1 0.733 0.0571 0.630 0.854
## 319 44 1 0.717 0.0582 0.611 0.840
## 360 43 1 0.700 0.0592 0.593 0.826
## 367 42 1 0.683 0.0601 0.575 0.812
## 410 41 1 0.667 0.0609 0.557 0.797
## 450 40 1 0.650 0.0616 0.540 0.783
## 470 39 1 0.633 0.0622 0.522 0.768
## 480 38 1 0.617 0.0628 0.505 0.753
## 510 37 1 0.600 0.0632 0.488 0.738plot(m.completo2.1,mark=c(16,1),lty=c(1,2),las=1,ylab= "Probabilidade de sobrevivencia", xlab="Tempo (s)",bty="l")
# legenda
legend("bottomleft",legend=levels(dados$Interacao),lty=c(1,2),pch=c(16,1),bty="n")
3.1.2 Agora a distancia
m.completo2.2 = survfit(Surv(tempo,Ataque)~Distancia, data=dados)
summary(m.completo2.2)## Call: survfit(formula = Surv(tempo, Ataque) ~ Distancia, data = dados)
##
## Distancia=longe
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 30 61 1 0.984 0.0163 0.952 1.000
## 70 60 1 0.967 0.0228 0.924 1.000
## 82 59 1 0.951 0.0277 0.898 1.000
## 112 58 1 0.934 0.0317 0.874 0.999
## 145 57 1 0.918 0.0351 0.852 0.990
## 165 56 1 0.902 0.0381 0.830 0.980
## 170 55 1 0.885 0.0408 0.809 0.969
## 200 54 1 0.869 0.0432 0.788 0.958
## 266 53 1 0.852 0.0454 0.768 0.946
## 311 52 1 0.836 0.0474 0.748 0.934
## 319 51 1 0.820 0.0492 0.729 0.922
## 353 50 1 0.803 0.0509 0.709 0.909
## 360 49 1 0.787 0.0524 0.691 0.897
## 367 48 1 0.770 0.0538 0.672 0.884
## 470 47 2 0.738 0.0563 0.635 0.857
## 510 45 2 0.705 0.0584 0.599 0.829
## 540 43 1 0.689 0.0593 0.582 0.815
##
## Distancia=perto
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 61 1 0.984 0.0163 0.952 1.000
## 10 60 5 0.902 0.0381 0.830 0.980
## 15 55 1 0.885 0.0408 0.809 0.969
## 20 54 1 0.869 0.0432 0.788 0.958
## 30 53 1 0.852 0.0454 0.768 0.946
## 33 52 1 0.836 0.0474 0.748 0.934
## 35 51 1 0.820 0.0492 0.729 0.922
## 40 50 3 0.770 0.0538 0.672 0.884
## 46 47 1 0.754 0.0551 0.653 0.870
## 52 46 1 0.738 0.0563 0.635 0.857
## 56 45 1 0.721 0.0574 0.617 0.843
## 60 44 1 0.705 0.0584 0.599 0.829
## 67 43 1 0.689 0.0593 0.582 0.815
## 95 42 1 0.672 0.0601 0.564 0.801
## 102 41 1 0.656 0.0608 0.547 0.786
## 115 40 1 0.639 0.0615 0.530 0.772
## 120 39 1 0.623 0.0621 0.512 0.757
## 125 38 1 0.607 0.0625 0.496 0.742
## 132 37 1 0.590 0.0630 0.479 0.727
## 170 36 1 0.574 0.0633 0.462 0.712
## 193 35 1 0.557 0.0636 0.446 0.697
## 202 34 1 0.541 0.0638 0.429 0.682
## 206 33 1 0.525 0.0639 0.413 0.666
## 231 32 1 0.508 0.0640 0.397 0.650
## 245 31 1 0.492 0.0640 0.381 0.635
## 265 30 1 0.475 0.0639 0.365 0.619
## 275 29 1 0.459 0.0638 0.350 0.603
## 300 28 1 0.443 0.0636 0.334 0.587
## 410 27 1 0.426 0.0633 0.319 0.570
## 450 26 1 0.410 0.0630 0.303 0.554
## 480 25 1 0.393 0.0625 0.288 0.537
## 510 24 1 0.377 0.0621 0.273 0.521
## 537 23 1 0.361 0.0615 0.258 0.504
## 545 22 1 0.344 0.0608 0.243 0.487plot(m.completo2.2,mark=c(16,1),lty=c(1,2),las=1,ylab= "Probabilidade de sobrevivencia", xlab="Tempo (s)",bty="l")
# legenda
legend("bottomleft",legend=levels(dados$Distancia),lty=c(1,2),pch=c(16,1),bty="n")
################3.1.3 Agora vamos juntar tudo
m.completo2.3 = survfit(Surv(tempo,Ataque)~Distancia+Interacao, data=dados)
summary(m.completo2.3)## Call: survfit(formula = Surv(tempo, Ataque) ~ Distancia + Interacao,
## data = dados)
##
## Distancia=longe, Interacao=hom
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 70 31 1 0.968 0.0317 0.908 1.000
## 82 30 1 0.935 0.0441 0.853 1.000
## 112 29 1 0.903 0.0531 0.805 1.000
## 145 28 1 0.871 0.0602 0.761 0.997
## 170 27 1 0.839 0.0661 0.719 0.979
## 266 26 1 0.806 0.0710 0.679 0.958
## 311 25 1 0.774 0.0751 0.640 0.936
## 353 24 1 0.742 0.0786 0.603 0.913
## 470 23 1 0.710 0.0815 0.567 0.889
## 510 22 1 0.677 0.0840 0.531 0.864
## 540 21 1 0.645 0.0859 0.497 0.838
##
## Distancia=longe, Interacao=nef
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 30 30 1 0.967 0.0328 0.905 1.000
## 165 29 1 0.933 0.0455 0.848 1.000
## 200 28 1 0.900 0.0548 0.799 1.000
## 319 27 1 0.867 0.0621 0.753 0.997
## 360 26 1 0.833 0.0680 0.710 0.978
## 367 25 1 0.800 0.0730 0.669 0.957
## 470 24 1 0.767 0.0772 0.629 0.934
## 510 23 1 0.733 0.0807 0.591 0.910
##
## Distancia=perto, Interacao=hom
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 31 1 0.968 0.0317 0.908 1.000
## 10 30 3 0.871 0.0602 0.761 0.997
## 15 27 1 0.839 0.0661 0.719 0.979
## 20 26 1 0.806 0.0710 0.679 0.958
## 33 25 1 0.774 0.0751 0.640 0.936
## 35 24 1 0.742 0.0786 0.603 0.913
## 40 23 2 0.677 0.0840 0.531 0.864
## 46 21 1 0.645 0.0859 0.497 0.838
## 56 20 1 0.613 0.0875 0.463 0.811
## 95 19 1 0.581 0.0886 0.431 0.783
## 102 18 1 0.548 0.0894 0.398 0.755
## 115 17 1 0.516 0.0898 0.367 0.726
## 132 16 1 0.484 0.0898 0.336 0.696
## 170 15 1 0.452 0.0894 0.306 0.666
## 202 14 1 0.419 0.0886 0.277 0.635
## 206 13 1 0.387 0.0875 0.249 0.603
## 265 12 1 0.355 0.0859 0.221 0.570
## 275 11 1 0.323 0.0840 0.194 0.537
## 510 10 1 0.290 0.0815 0.167 0.503
## 537 9 1 0.258 0.0786 0.142 0.469
## 545 8 1 0.226 0.0751 0.118 0.433
##
## Distancia=perto, Interacao=nef
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 10 30 2 0.933 0.0455 0.848 1.000
## 30 28 1 0.900 0.0548 0.799 1.000
## 40 27 1 0.867 0.0621 0.753 0.997
## 52 26 1 0.833 0.0680 0.710 0.978
## 60 25 1 0.800 0.0730 0.669 0.957
## 67 24 1 0.767 0.0772 0.629 0.934
## 120 23 1 0.733 0.0807 0.591 0.910
## 125 22 1 0.700 0.0837 0.554 0.885
## 193 21 1 0.667 0.0861 0.518 0.859
## 231 20 1 0.633 0.0880 0.482 0.832
## 245 19 1 0.600 0.0894 0.448 0.804
## 300 18 1 0.567 0.0905 0.414 0.775
## 410 17 1 0.533 0.0911 0.382 0.745
## 450 16 1 0.500 0.0913 0.350 0.715
## 480 15 1 0.467 0.0911 0.318 0.684plot(m.completo2.3,mark=c(16,1,16,1), col=1:4, lty=c(1,2,1,2),las=1,ylab= "Probabilidade de ataque", xlab="Tempo (s)",bty="l", lwd=2)
# legenda
legend("bottomleft",legend=c("Hom-Longe", "Nef-Longe", "Hom-Perto", "Nef-Perto"),lty=c(1,2,1,2),pch=c(16,1, 16, 1),bty="n", col=1:4)
3.1.4 Com o ggplot2
ggsurvplot(
m.completo2.3, # survfit object with calculated statistics.
xlab = "Time in seconds", # customize X axis label.
break.time.by = 100, # break X axis in time intervals by 100.
ncensor.plot = TRUE, # plot the number of censored subjects at time t
legend.title="Tratamento",
legend.labs = c("Longe - Hom", "Longe - NEF", "Perto - Hom", "Perto - NEF") # change legend labels.
)