TP 2-Corrigé
EXO1
#1)
library(survival)
head(aml)## time status x
## 1 9 1 Maintained
## 2 13 1 Maintained
## 3 13 0 Maintained
## 4 18 1 Maintained
## 5 23 1 Maintained
## 6 28 0 Maintaineddim(aml)## [1] 23 3attach(aml)
#2)
base=Surv(time,status==1)
base## [1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+ 5 5 8
## [15] 8 12 16+ 23 27 30 33 43 451-sum(status)/23## [1] 0.2173913#21.74% d'obs censurées
#3)
#Survie globale
fit=survfit(base~1,conf.type="plain")
summary(fit)## Call: survfit(formula = base ~ 1, conf.type = "plain")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 23 2 0.9130 0.0588 0.7979 1.000
## 8 21 2 0.8261 0.0790 0.6712 0.981
## 9 19 1 0.7826 0.0860 0.6140 0.951
## 12 18 1 0.7391 0.0916 0.5597 0.919
## 13 17 1 0.6957 0.0959 0.5076 0.884
## 18 14 1 0.6460 0.1011 0.4477 0.844
## 23 13 2 0.5466 0.1073 0.3364 0.757
## 27 11 1 0.4969 0.1084 0.2844 0.709
## 30 9 1 0.4417 0.1095 0.2270 0.656
## 31 8 1 0.3865 0.1089 0.1731 0.600
## 33 7 1 0.3313 0.1064 0.1227 0.540
## 34 6 1 0.2761 0.1020 0.0762 0.476
## 43 5 1 0.2208 0.0954 0.0339 0.408
## 45 4 1 0.1656 0.0860 0.0000 0.334
## 48 2 1 0.0828 0.0727 0.0000 0.225plot(fit,xlim=c(0,100))
#Survie par groupe
fitgr=survfit(base~x)
summary(fitgr)## Call: survfit(formula = base ~ x)
##
## x=Maintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 9 11 1 0.909 0.0867 0.7541 1.000
## 13 10 1 0.818 0.1163 0.6192 1.000
## 18 8 1 0.716 0.1397 0.4884 1.000
## 23 7 1 0.614 0.1526 0.3769 0.999
## 31 5 1 0.491 0.1642 0.2549 0.946
## 34 4 1 0.368 0.1627 0.1549 0.875
## 48 2 1 0.184 0.1535 0.0359 0.944
##
## x=Nonmaintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 12 2 0.8333 0.1076 0.6470 1.000
## 8 10 2 0.6667 0.1361 0.4468 0.995
## 12 8 1 0.5833 0.1423 0.3616 0.941
## 23 6 1 0.4861 0.1481 0.2675 0.883
## 27 5 1 0.3889 0.1470 0.1854 0.816
## 30 4 1 0.2917 0.1387 0.1148 0.741
## 33 3 1 0.1944 0.1219 0.0569 0.664
## 43 2 1 0.0972 0.0919 0.0153 0.620
## 45 1 1 0.0000 NaN NA NA#4)
plot(fitgr, col = c(1,2))
title(main="Durées de survie de patients atteints d'aml", col.main="black", font.main=4)
title(xlab="Semaines")
title(ylab="Estimation de la fonction de survie")
legend(50, 1, c("Maintained","Nonmaintened"), cex=0.8,
col=c("black","red"),lty=1)
#5)
print(fit, print.rmean=TRUE)## Call: survfit(formula = base ~ 1, conf.type = "plain")
##
## records n.max n.start events *rmean *se(rmean)
## 23.00 23.00 23.00 18.00 36.36 9.85
## median 0.95LCL 0.95UCL
## 27.00 18.00 34.00
## * restricted mean with upper limit = 161#mean=36.36, se(mean)=9.85 mais dernière valeur prise=161 !!
print(fit, print.rmean=TRUE,rmean = 48)## Call: survfit(formula = base ~ 1, conf.type = "plain")
##
## records n.max n.start events *rmean *se(rmean)
## 23.00 23.00 23.00 18.00 27.01 3.21
## median 0.95LCL 0.95UCL
## 27.00 18.00 34.00
## * restricted mean with upper limit = 48#mean=27.01, se(mean)=3.21
print(fitgr, print.rmean=TRUE)## Call: survfit(formula = base ~ x)
##
## records n.max n.start events *rmean *se(rmean) median
## x=Maintained 11 11 11 7 42.0 11.26 31
## x=Nonmaintained 12 12 12 11 22.7 4.18 23
## 0.95LCL 0.95UCL
## x=Maintained 18 NA
## x=Nonmaintained 8 NA
## * restricted mean with upper limit = 103print(fitgr, print.rmean=TRUE,rmean=48)## Call: survfit(formula = base ~ x)
##
## records n.max n.start events *rmean *se(rmean) median
## x=Maintained 11 11 11 7 31.8 4.53 31
## x=Nonmaintained 12 12 12 11 22.7 4.18 23
## 0.95LCL 0.95UCL
## x=Maintained 18 NA
## x=Nonmaintained 8 NA
## * restricted mean with upper limit = 48Global - 1st quartile=12 - Median=27 - 3rd quartile=43 - Moy=27.01 - Ecart-type=3.21
Maintened - 1st quartile=18 - Mediane=31 - 3rd quartile=48 - Moy=31.8 - Ecart-type=4.53
Non-maintened - 1st quartile=8 - Mediane=23 - 3rd quartile=33 - Moy=22.7 - Ecart-type=4.18
plot(fitgr, col = c(1,2))
title(main="Durées de survie de patients atteints d'aml", col.main="black", font.main=4)
title(xlab="Semaines")
title(ylab="Estimation de la fonction de survie")
legend(50, 1, c("Maintained","Nonmaintened"), cex=0.8,
col=c("black","red"),lty=1)
xM=aml[aml$x=="Maintained",]
xNM=aml[aml$x=="Nonmaintained",]
fitM=survfit(Surv(xM$time,xM$status==1)~1)
fitNM=survfit(Surv(xNM$time,xNM$status==1)~1)
#Courbes survie par groupe+intervalle de confiance à 90%
lines(stepfun(fitM$time,c(1,fitM$lower)),do.points=FALSE,col=1,lty=2)
lines(stepfun(fitM$time,c(1,fitM$upper)),do.points=FALSE,col=1,lty=2)
lines(stepfun(fitNM$time,c(1,fitNM$lower)),do.points=FALSE,col=2,lty=2)
lines(stepfun(fitNM$time,c(1,fitNM$upper)),do.points=FALSE,col=2,lty=2)
#7)
#Tests non-paramétriques de comparaison des survies
survdiff(base~x,rho=0)## Call:
## survdiff(formula = base ~ x, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained 11 7 10.69 1.27 3.4
## x=Nonmaintained 12 11 7.31 1.86 3.4
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.0653# Chisq= 3.4 on 1 degrees of freedom, p= 0.0653
survdiff(base~x,rho=1)## Call:
## survdiff(formula = base ~ x, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained 11 3.85 6.14 0.859 2.78
## x=Nonmaintained 12 7.18 4.88 1.081 2.78
##
## Chisq= 2.8 on 1 degrees of freedom, p= 0.0955# Chisq= 2.8 on 1 degrees of freedom, p= 0.0955
#8)
#Risque cumulé pour tous les patients#
str(fit)## List of 13
## $ n : int 23
## $ time : int [1:18] 5 8 9 12 13 16 18 23 27 28 ...
## $ n.risk : num [1:18] 23 21 19 18 17 15 14 13 11 10 ...
## $ n.event : num [1:18] 2 2 1 1 1 0 1 2 1 0 ...
## $ n.censor : num [1:18] 0 0 0 0 1 1 0 0 0 1 ...
## $ surv : num [1:18] 0.913 0.826 0.783 0.739 0.696 ...
## $ type : chr "right"
## $ std.err : num [1:18] 0.0643 0.0957 0.1099 0.1239 0.1379 ...
## $ upper : num [1:18] 1 0.981 0.951 0.919 0.884 ...
## $ lower : num [1:18] 0.798 0.671 0.614 0.56 0.508 ...
## $ conf.type: chr "plain"
## $ conf.int : num 0.95
## $ call : language survfit(formula = base ~ 1, conf.type = "plain")
## - attr(*, "class")= chr "survfit"risqueCum=cumsum(fit$n.event/fit$n.risk)
plot.stepfun(stepfun(fit$time,c(0,risqueCum)),do.points=FALSE,main="Estimations du risque cumulé",xlab="Temps",ylab="Risque cumulé estimé")
#, col.lab=rgb(0,0.5,0)
H=-log(fit$surv)
lines(stepfun(fit$time,c(0,H)),do.points=FALSE,col="red")
legend(100, 0.2, c("Nelson-Aalen","Breslow"), cex=0.8,
col=c("black","red"),lty=1)
#9)
#Risque cumulé par x#
risqueCumM=cumsum(fitM$n.event/fitM$n.risk)
risqueCumNM=cumsum(fitNM$n.event/fitNM$n.risk)
plot.stepfun(stepfun(fitNM$time,c(0,risqueCumNM)),main="Estimations du risque cumulé par maintien",xlab="Semaines",ylab="Risque cumulé estimé",do.points=FALSE,col="red")
lines(stepfun(fitM$time,c(0,risqueCumM)),do.points=FALSE)
legend(33, 0.2, c("Maintened","Non-maintened"), cex=0.8,
col=c("black","red"),lty=1)
detach(aml)Exo 2
library(ISwR)##
## Attaching package: 'ISwR'## The following object is masked from 'package:survival':
##
## lunghead(melanom)## no status days ulc thick sex
## 1 789 3 10 1 676 2
## 2 13 3 30 2 65 2
## 3 97 2 35 2 134 2
## 4 16 3 99 2 290 1
## 5 21 1 185 1 1208 2
## 6 469 1 204 1 484 2attach(melanom)status=1 décès, status=2 ou 3 censure
base=Surv(days,status==1)
#3)
fit=survfit(base~1)
summary(fit)## Call: survfit(formula = base ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 185 201 1 0.995 0.00496 0.985 1.000
## 204 200 1 0.990 0.00700 0.976 1.000
## 210 199 1 0.985 0.00855 0.968 1.000
## 232 198 1 0.980 0.00985 0.961 1.000
## 279 196 1 0.975 0.01100 0.954 0.997
## 295 195 1 0.970 0.01202 0.947 0.994
## 386 193 1 0.965 0.01297 0.940 0.991
## 426 192 1 0.960 0.01384 0.933 0.988
## 469 191 1 0.955 0.01465 0.927 0.984
## 529 189 1 0.950 0.01542 0.920 0.981
## 621 188 1 0.945 0.01615 0.914 0.977
## 629 187 1 0.940 0.01683 0.907 0.973
## 659 186 1 0.935 0.01748 0.901 0.970
## 667 185 1 0.930 0.01811 0.895 0.966
## 718 184 1 0.925 0.01870 0.889 0.962
## 752 183 1 0.920 0.01927 0.883 0.958
## 779 182 1 0.915 0.01981 0.877 0.954
## 793 181 1 0.910 0.02034 0.871 0.950
## 817 180 1 0.904 0.02084 0.865 0.946
## 833 178 1 0.899 0.02134 0.859 0.942
## 858 177 1 0.894 0.02181 0.853 0.938
## 869 176 1 0.889 0.02227 0.847 0.934
## 872 175 1 0.884 0.02272 0.841 0.930
## 967 174 1 0.879 0.02315 0.835 0.926
## 977 173 1 0.874 0.02357 0.829 0.921
## 982 172 1 0.869 0.02397 0.823 0.917
## 1041 171 1 0.864 0.02436 0.817 0.913
## 1055 170 1 0.859 0.02474 0.812 0.909
## 1062 169 1 0.854 0.02511 0.806 0.904
## 1075 168 1 0.849 0.02547 0.800 0.900
## 1156 167 1 0.844 0.02582 0.794 0.896
## 1228 166 1 0.838 0.02616 0.789 0.891
## 1252 165 1 0.833 0.02649 0.783 0.887
## 1271 164 1 0.828 0.02681 0.777 0.883
## 1312 163 1 0.823 0.02713 0.772 0.878
## 1435 161 1 0.818 0.02744 0.766 0.874
## 1506 159 1 0.813 0.02774 0.760 0.869
## 1516 155 1 0.808 0.02805 0.755 0.865
## 1548 152 1 0.802 0.02837 0.749 0.860
## 1560 150 1 0.797 0.02868 0.743 0.855
## 1584 148 1 0.792 0.02899 0.737 0.851
## 1621 146 1 0.786 0.02929 0.731 0.846
## 1667 137 1 0.780 0.02963 0.725 0.841
## 1690 134 1 0.775 0.02998 0.718 0.836
## 1726 131 1 0.769 0.03033 0.712 0.831
## 1933 110 1 0.762 0.03085 0.704 0.825
## 2061 95 1 0.754 0.03155 0.694 0.818
## 2062 94 1 0.746 0.03221 0.685 0.812
## 2103 90 1 0.737 0.03290 0.676 0.805
## 2108 88 1 0.729 0.03358 0.666 0.798
## 2256 80 1 0.720 0.03438 0.656 0.791
## 2388 75 1 0.710 0.03523 0.645 0.783
## 2467 69 1 0.700 0.03619 0.633 0.775
## 2565 63 1 0.689 0.03729 0.620 0.766
## 2782 57 1 0.677 0.03854 0.605 0.757
## 3042 52 1 0.664 0.03994 0.590 0.747
## 3338 35 1 0.645 0.04307 0.566 0.735plot(fit)
#4)
fitSex=survfit(base~sex)
summary(fitSex)## Call: survfit(formula = base ~ sex)
##
## sex=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 279 124 1 0.992 0.00803 0.976 1.000
## 295 123 1 0.984 0.01131 0.962 1.000
## 386 121 1 0.976 0.01384 0.949 1.000
## 469 120 1 0.968 0.01593 0.937 0.999
## 667 119 1 0.959 0.01775 0.925 0.995
## 817 118 1 0.951 0.01937 0.914 0.990
## 833 116 1 0.943 0.02087 0.903 0.985
## 858 115 1 0.935 0.02224 0.892 0.980
## 869 114 1 0.927 0.02351 0.882 0.974
## 872 113 1 0.919 0.02469 0.871 0.968
## 982 112 1 0.910 0.02580 0.861 0.962
## 1055 111 1 0.902 0.02684 0.851 0.956
## 1156 110 1 0.894 0.02782 0.841 0.950
## 1252 109 1 0.886 0.02875 0.831 0.944
## 1271 108 1 0.878 0.02963 0.821 0.938
## 1312 107 1 0.869 0.03046 0.812 0.931
## 1548 102 1 0.861 0.03134 0.802 0.924
## 1560 100 1 0.852 0.03218 0.791 0.918
## 1621 97 1 0.843 0.03303 0.781 0.911
## 1667 90 1 0.834 0.03396 0.770 0.903
## 1726 86 1 0.824 0.03493 0.759 0.896
## 1933 74 1 0.813 0.03619 0.745 0.887
## 2062 63 1 0.800 0.03785 0.729 0.878
## 2108 60 1 0.787 0.03950 0.713 0.868
## 2256 54 1 0.772 0.04137 0.695 0.858
## 2467 45 1 0.755 0.04386 0.674 0.846
## 3042 34 1 0.733 0.04787 0.645 0.833
## 3338 25 1 0.704 0.05419 0.605 0.818
##
## sex=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 185 76 1 0.987 0.0131 0.962 1.000
## 204 75 1 0.974 0.0184 0.938 1.000
## 210 74 1 0.961 0.0223 0.918 1.000
## 232 73 1 0.947 0.0256 0.898 0.999
## 426 72 1 0.934 0.0284 0.880 0.992
## 529 70 1 0.921 0.0310 0.862 0.984
## 621 69 1 0.908 0.0333 0.845 0.975
## 629 68 1 0.894 0.0354 0.827 0.966
## 659 67 1 0.881 0.0373 0.811 0.957
## 718 66 1 0.867 0.0390 0.794 0.947
## 752 65 1 0.854 0.0407 0.778 0.938
## 779 64 1 0.841 0.0422 0.762 0.928
## 793 63 1 0.827 0.0435 0.746 0.917
## 967 62 1 0.814 0.0448 0.731 0.907
## 977 61 1 0.801 0.0461 0.715 0.896
## 1041 60 1 0.787 0.0472 0.700 0.886
## 1062 59 1 0.774 0.0482 0.685 0.875
## 1075 58 1 0.761 0.0492 0.670 0.864
## 1228 57 1 0.747 0.0501 0.655 0.852
## 1435 55 1 0.734 0.0510 0.640 0.841
## 1506 53 1 0.720 0.0519 0.625 0.829
## 1516 51 1 0.706 0.0528 0.610 0.817
## 1584 50 1 0.692 0.0536 0.594 0.805
## 1690 47 1 0.677 0.0544 0.578 0.792
## 2061 32 1 0.656 0.0567 0.554 0.777
## 2103 29 1 0.633 0.0591 0.527 0.760
## 2388 25 1 0.608 0.0619 0.498 0.742
## 2565 22 1 0.580 0.0650 0.466 0.723
## 2782 21 1 0.553 0.0675 0.435 0.702plot(fitSex, col = c(1,2))
title(main="Durées de survie de patients après opération d'un mélanome", col.main="black", font.main=4)
title(xlab="Semaines")
title(ylab="Estimation de la fonction de survie")
legend(4000, 1, c("Femme","Homme"), cex=0.8,
col=c("black","red"),lty=1)
#Les femmes semblent vivrent plus longtemps que les hommes.
#5)
#Tests non-paramétriques de comparaison des survies
survdiff(base~sex,data=melanom,rho=0)## Call:
## survdiff(formula = base ~ sex, data = melanom, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 126 28 37.1 2.25 6.47
## sex=2 79 29 19.9 4.21 6.47
##
## Chisq= 6.5 on 1 degrees of freedom, p= 0.011# Chisq= 6.5 on 1 degrees of freedom, p= 0.011
survdiff(base~sex,data=melanom,rho=1)## Call:
## survdiff(formula = base ~ sex, data = melanom, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 126 23.4 31.6 2.14 7.09
## sex=2 79 25.2 17.0 3.98 7.09
##
## Chisq= 7.1 on 1 degrees of freedom, p= 0.00776# Chisq= 7.1 on 1 degrees of freedom, p= 0.00776
#Les tests sont significatifs : il y a une différence significative entre la durée de vie des hommes et des femmes.
#(avec un risque d'erreur de 1% environ pour chaque test)
#6)
fitUlc=survfit(base~ulc)
summary(fitUlc)## Call: survfit(formula = base ~ ulc)
##
## ulc=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 185 89 1 0.989 0.0112 0.967 1.000
## 204 88 1 0.978 0.0157 0.947 1.000
## 210 87 1 0.966 0.0191 0.930 1.000
## 232 86 1 0.955 0.0220 0.913 0.999
## 279 84 1 0.944 0.0245 0.897 0.993
## 295 83 1 0.932 0.0267 0.881 0.986
## 386 81 1 0.921 0.0287 0.866 0.979
## 426 80 1 0.909 0.0306 0.851 0.971
## 469 79 1 0.898 0.0323 0.837 0.963
## 529 77 1 0.886 0.0339 0.822 0.955
## 621 76 1 0.874 0.0354 0.808 0.947
## 629 75 1 0.863 0.0368 0.794 0.938
## 659 74 1 0.851 0.0381 0.780 0.929
## 667 73 1 0.839 0.0393 0.766 0.920
## 718 72 1 0.828 0.0405 0.752 0.911
## 752 71 1 0.816 0.0416 0.739 0.902
## 779 70 1 0.805 0.0426 0.725 0.892
## 793 69 1 0.793 0.0435 0.712 0.883
## 833 67 1 0.781 0.0444 0.699 0.873
## 967 66 1 0.769 0.0453 0.685 0.863
## 977 65 1 0.757 0.0461 0.672 0.853
## 982 64 1 0.746 0.0469 0.659 0.843
## 1055 63 1 0.734 0.0476 0.646 0.833
## 1062 62 1 0.722 0.0483 0.633 0.823
## 1075 61 1 0.710 0.0490 0.620 0.813
## 1156 60 1 0.698 0.0495 0.608 0.802
## 1228 59 1 0.686 0.0501 0.595 0.792
## 1252 58 1 0.675 0.0506 0.582 0.781
## 1271 57 1 0.663 0.0511 0.570 0.771
## 1312 56 1 0.651 0.0515 0.557 0.760
## 1506 55 1 0.639 0.0519 0.545 0.749
## 1516 53 1 0.627 0.0523 0.532 0.738
## 1621 51 1 0.615 0.0527 0.520 0.727
## 1667 50 1 0.602 0.0531 0.507 0.716
## 1726 49 1 0.590 0.0534 0.494 0.705
## 2108 32 1 0.572 0.0548 0.474 0.690
## 2256 29 1 0.552 0.0564 0.452 0.674
## 2467 24 1 0.529 0.0585 0.426 0.657
## 2565 19 1 0.501 0.0617 0.394 0.638
## 3042 16 1 0.470 0.0653 0.358 0.617
## 3338 12 1 0.431 0.0706 0.312 0.594
##
## ulc=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 817 112 1 0.991 0.00889 0.974 1.000
## 858 111 1 0.982 0.01251 0.958 1.000
## 869 110 1 0.973 0.01526 0.944 1.000
## 872 109 1 0.964 0.01754 0.931 0.999
## 1041 108 1 0.955 0.01951 0.918 0.994
## 1435 106 1 0.946 0.02131 0.905 0.989
## 1548 101 1 0.937 0.02307 0.893 0.983
## 1560 99 1 0.928 0.02470 0.880 0.977
## 1584 97 1 0.918 0.02623 0.868 0.971
## 1690 85 1 0.907 0.02806 0.854 0.964
## 1933 71 1 0.894 0.03043 0.837 0.956
## 2061 62 1 0.880 0.03318 0.817 0.947
## 2062 61 1 0.866 0.03564 0.798 0.938
## 2103 58 1 0.851 0.03802 0.779 0.928
## 2388 49 1 0.833 0.04102 0.757 0.918
## 2782 41 1 0.813 0.04477 0.730 0.906plot(fitUlc, col = c(1,2))
title(main="Durées de survie de patients après opération d'un mélanome", col.main="black", font.main=4)
title(xlab="Semaines")
title(ylab="Estimation de la fonction de survie")
legend(4000, 1, c("Présent","Absent"), cex=0.8,
col=c("black","red"),lty=1)
#Tests non-paramétriques de comparaison des survies
survdiff(base~ulc,data=melanom,rho=0)## Call:
## survdiff(formula = base ~ ulc, data = melanom, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ulc=1 90 41 21.2 18.5 29.6
## ulc=2 115 16 35.8 10.9 29.6
##
## Chisq= 29.6 on 1 degrees of freedom, p= 5.41e-08# Chisq= 29.6 on 1 degrees of freedom, p= 5.41e-08
survdiff(base~ulc,data=melanom,rho=1)## Call:
## survdiff(formula = base ~ ulc, data = melanom, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ulc=1 90 35.7 18.3 16.7 31.1
## ulc=2 115 12.9 30.4 10.1 31.1
##
## Chisq= 31.1 on 1 degrees of freedom, p= 2.45e-08# Chisq= 31.1 on 1 degrees of freedom, p= 2.45e-08
#Différence hautement significative entre la survie de patients ayant ulcération et ceux n'ayant pas ulcération.
#7)
#On se demande si l'effet sexe n'est pas amplifié par l'effet ulcération.
#Ce serait le cas par exemple si les hommes avaient plus tendance à avoir de l'ulcération que les femmes.
TabProp=prop.table(table(sex,ulc),margin=2)
rownames(TabProp)=c("Femme","Homme")
barplot(t(TabProp),beside=TRUE,col=c(1,8))
legend(4, 0.6, c("ulc present","ulc absent"), cex=0.8, col=c(1,8),lty=1)
#Parmi ceux dont ulc est présent il y a 52,22% de femmes, alors que parmi ceux
#dont ulc est absent il y a 68,7% de femmes etc...
#Test du chi-deux
chisq.test(sex,ulc)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sex and ulc
## X-squared = 5.1099, df = 1, p-value = 0.02379#X-squared = 5.1099, df = 1, p-value = 0.02379
#Il y a bien un lien entre sexe et ulcération (avec un risque de 2.4% d'erreur)!
fitSexUlc=survfit(base~sex+strata(ulc),data=melanom)
survdiff(base~sex+strata(ulc),data=melanom,rho=0)## Call:
## survdiff(formula = base ~ sex + strata(ulc), data = melanom,
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 126 28 34.7 1.28 3.31
## sex=2 79 29 22.3 1.99 3.31
##
## Chisq= 3.3 on 1 degrees of freedom, p= 0.0687#Chisq= 3.3 on 1 degrees of freedom, p= 0.0687
#Test beaucoup moins significatif que pour juste sex.
#Les différences entre sexes sont réduites après avoir pris en
#compte l'ulcération!!
#8)
summary(thick)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10 97 194 292 356 1742# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 10 97 194 292 356 1742
thick1 = (thick<97)
thick2 = (thick>=97)&(thick<194)
thick3 = (thick>=194)&(thick<356)
thick4 = (thick>=356)
fitTh=survfit(base~thick1+thick2+thick3+thick4)
summary(fitTh)## Call: survfit(formula = base ~ thick1 + thick2 + thick3 + thick4)
##
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 185 51 1 0.980 0.0194 0.943 1.000
## 204 50 1 0.961 0.0272 0.909 1.000
## 210 49 1 0.941 0.0329 0.879 1.000
## 232 48 1 0.922 0.0376 0.851 0.998
## 279 47 1 0.902 0.0416 0.824 0.987
## 295 46 1 0.882 0.0451 0.798 0.975
## 386 45 1 0.863 0.0482 0.773 0.963
## 426 44 1 0.843 0.0509 0.749 0.949
## 529 42 1 0.823 0.0535 0.725 0.935
## 621 41 1 0.803 0.0559 0.701 0.920
## 629 40 1 0.783 0.0580 0.677 0.905
## 659 39 1 0.763 0.0598 0.654 0.890
## 667 38 1 0.743 0.0615 0.631 0.874
## 752 37 1 0.723 0.0631 0.609 0.858
## 779 36 1 0.703 0.0644 0.587 0.841
## 793 35 1 0.683 0.0656 0.565 0.824
## 858 33 1 0.662 0.0668 0.543 0.807
## 967 32 1 0.641 0.0679 0.521 0.789
## 982 31 1 0.620 0.0688 0.499 0.771
## 1041 30 1 0.600 0.0695 0.478 0.753
## 1062 29 1 0.579 0.0701 0.457 0.734
## 1252 28 1 0.558 0.0706 0.436 0.715
## 1312 27 1 0.538 0.0710 0.415 0.696
## 1506 26 1 0.517 0.0712 0.395 0.677
## 1621 24 1 0.496 0.0714 0.374 0.657
## 1667 23 1 0.474 0.0715 0.353 0.637
##
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 469 52 1 0.981 0.0190 0.944 1.000
## 718 51 1 0.962 0.0267 0.911 1.000
## 833 50 1 0.942 0.0323 0.881 1.000
## 869 49 1 0.923 0.0370 0.853 0.998
## 1055 48 1 0.904 0.0409 0.827 0.988
## 1075 47 1 0.885 0.0443 0.802 0.976
## 1228 46 1 0.865 0.0473 0.777 0.963
## 1271 45 1 0.846 0.0500 0.754 0.950
## 1435 44 1 0.827 0.0525 0.730 0.936
## 1516 42 1 0.807 0.0548 0.707 0.922
## 1560 41 1 0.788 0.0569 0.684 0.907
## 1933 32 1 0.763 0.0602 0.654 0.891
## 2061 28 1 0.736 0.0639 0.621 0.872
## 2062 27 1 0.708 0.0671 0.588 0.853
## 2256 24 1 0.679 0.0705 0.554 0.832
## 2467 21 1 0.647 0.0742 0.516 0.810
## 2565 20 1 0.614 0.0772 0.480 0.786
## 2782 18 1 0.580 0.0801 0.443 0.760
## 3042 16 1 0.544 0.0829 0.403 0.733
## 3338 11 1 0.494 0.0889 0.348 0.703
##
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 872 53 1 0.981 0.0187 0.945 1.000
## 977 52 1 0.962 0.0262 0.912 1.000
## 1156 51 1 0.943 0.0317 0.883 1.000
## 1548 47 1 0.923 0.0369 0.854 0.998
## 1690 37 1 0.898 0.0435 0.817 0.988
## 1726 36 1 0.873 0.0489 0.783 0.975
## 2108 20 1 0.830 0.0630 0.715 0.963
##
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 817 43 1 0.977 0.0230 0.933 1
## 1584 39 1 0.952 0.0334 0.889 1
## 2103 27 1 0.916 0.0472 0.828 1
## 2388 24 1 0.878 0.0587 0.770 1plot(fitTh, col = c(1,2,3,4))
title(main="Durées de survie de patients après opération d'un mélanome", col.main="black", font.main=4)
title(xlab="Semaines")
title(ylab="Estimation de la fonction de survie")
legend(4000, 0.3, c("thick>356","194<thick<356","97<thick<194","thick<97"), cex=0.8,
col=c("black","red","green","blue"),lty=1)
survdiff(base~thick1+thick2+thick3+thick4,data=melanom,rho=0)## Call:
## survdiff(formula = base ~ thick1 + thick2 + thick3 + thick4,
## data = melanom, rho = 0)
##
## N Observed
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 52 26
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 54 20
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 54 7
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 45 4
## Expected (O-E)^2/E
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 10.7 21.71
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 16.0 1.01
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 15.8 4.89
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 14.5 7.61
## (O-E)^2/V
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 26.89
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 1.41
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 6.79
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 10.25
##
## Chisq= 35.4 on 3 degrees of freedom, p= 1.01e-07# Chisq= 35.4 on 3 degrees of freedom, p= 1.01e-07
survdiff(base~thick1+thick2+thick3+thick4,data=melanom,rho=1)## Call:
## survdiff(formula = base ~ thick1 + thick2 + thick3 + thick4,
## data = melanom, rho = 1)
##
## N Observed
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 52 23.79
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 54 15.97
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 54 5.72
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 45 3.17
## Expected (O-E)^2/E
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 9.27 22.732
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 13.58 0.422
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 13.55 4.530
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 12.25 6.727
## (O-E)^2/V
## thick1=FALSE, thick2=FALSE, thick3=FALSE, thick4=TRUE 32.406
## thick1=FALSE, thick2=FALSE, thick3=TRUE , thick4=FALSE 0.681
## thick1=FALSE, thick2=TRUE , thick3=FALSE, thick4=FALSE 7.277
## thick1=TRUE, thick2=FALSE, thick3=FALSE, thick4=FALSE 10.523
##
## Chisq= 39.8 on 3 degrees of freedom, p= 1.17e-08# Chisq= 35.4 on 3 degrees of freedom, p= 1.01e-07
detach(melanom)