Exo 3
library(survival)
head(kidney)
## id time status age sex disease frail
## 1 1 8 1 28 1 Other 2.3
## 2 1 16 1 28 1 Other 2.3
## 3 2 23 1 48 2 GN 1.9
## 4 2 13 0 48 2 GN 1.9
## 5 3 22 1 32 1 Other 1.2
## 6 3 28 1 32 1 Other 1.2
dim(kidney)
## [1] 76 7
? kidney
attach(kidney)
#(1) Nombre et prop. de données censurées sur l'échantillon total
sum(1-status)
## [1] 18
sum(1-status)/76
## [1] 0.2368421
#0.2368421
#Nombre et prop. de données censurées par groupe
sum((1-status)[sex==1])
## [1] 2
sum((1-status)[sex==1])/(sum(sex==1))
## [1] 0.1
#0.1
sum((1-status)[sex==2])
## [1] 16
sum((1-status)[sex==2])/(sum(sex==2))
## [1] 0.2857143
#0.2857143
#(1bis) Stat univariées
table(sex)
## sex
## 1 2
## 20 56
boxplot(age)
summary(age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.0 34.0 45.5 43.7 54.0 69.0
table(disease)
## disease
## Other GN AN PKD
## 26 18 24 8
#(2) Comparaison de l'âge par sexe
boxplot(age~sex,data=kidney)
title(main="Boite a moustache de l'age en fonction du sexe")
#ageH : age des hommes
#ageF : age des femmes
ageH=kidney$age[sex==1]
ageF=kidney$age[sex==2]
#(3) Age/variance moyen des hommes et femmes
mean(ageH)
## [1] 43.5
#43.5
mean(ageF)
## [1] 43.76786
#43.76786
var(ageH)
## [1] 237.8421
sd(ageH)
## [1] 15.42213
#15.42213
var(ageF)
## [1] 214.1088
sd(ageF)
## [1] 14.63246
#14.63246
summary(ageH)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.0 32.0 49.0 43.5 53.0 60.0
summary(ageF)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.00 34.75 44.00 43.77 54.25 69.00
#(4) Comparaison moyenne âge H/F
length(ageH)
## [1] 20
#20
length(ageF)
## [1] 56
#56
shapiro.test(ageH)
##
## Shapiro-Wilk normality test
##
## data: ageH
## W = 0.86996, p-value = 0.01173
#p-value = 0.01173
shapiro.test(ageF)
##
## Shapiro-Wilk normality test
##
## data: ageF
## W = 0.94174, p-value = 0.009228
#p-value = 0.009228
#normalité des âges n'est pas vérifié. Le t-test effectué est alors asymptotique. Attention car seulement 20 individus pour les hommes!!
#Attention test de Fisher (var.test(ageH,ageF)) non valide car variables non gaussiennes...!
bartlett.test(age~sex)#test des variances moins sensible à la gaussianité des variables.
##
## Bartlett test of homogeneity of variances
##
## data: age by sex
## Bartlett's K-squared = 0.077859, df = 1, p-value = 0.7802
#p-value = 0.7802. Variances semblent égales
t.test(ageH,ageF,paired=FALSE)
##
## Welch Two Sample t-test
##
## data: ageH and ageF
## t = -0.067568, df = 32.037, p-value = 0.9465
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -8.342454 7.806740
## sample estimates:
## mean of x mean of y
## 43.50000 43.76786
#p-value = 0.9465
#Espérances semblent égales.
###Tests non-paramétriques (plus approprié ici)###
ks.test(ageH,ageF)
## Warning in ks.test(ageH, ageF): cannot compute exact p-value with ties
##
## Two-sample Kolmogorov-Smirnov test
##
## data: ageH and ageF
## D = 0.13571, p-value = 0.9489
## alternative hypothesis: two-sided
#p-value = 0.949 pb d'ex-aequo, la p-valeur n'est pas "exacte"
wilcox.test(ageH,ageF)
##
## Wilcoxon rank sum test with continuity correction
##
## data: ageH and ageF
## W = 563.5, p-value = 0.9717
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.9717
#Dans tous les cas on ne rejette pas (H0), égalité des âges H/F !
#(5) Représentation graphique de disease/sexe
par(mfrow=c(2,1))
barplot(prop.table(table(disease[sex==1])),main="Proportions des types de maladies parmi les hommes")
barplot(prop.table(table(disease[sex==2])),main="Proportions des types de maladies parmi les femmes")
#On veut faire un test du chi-deux pour comparer disease/sexe
#Tableau effectif
table(disease,sex)
## sex
## disease 1 2
## Other 6 20
## GN 6 12
## AN 4 20
## PKD 4 4
#Vérification eff. théoriques suffisament grands
(4+4)*(6+6+4+4)/sum(table(disease,sex))
## [1] 2.105263
#=2.105263<<5 !!! Les effectifs PKD sont trop petits !
(6+12)*(6+6+4+4)/sum(table(disease,sex))
## [1] 4.736842
#=4.736842<5 ! Les effectifs GN sont trop petits !
#On regroupe les PKD et GN d'eff. trop petits avec Other
tab=matrix(c(16,4,36,20),ncol=2)
rownames(tab)=c("Other","AN")
colnames(tab)=c("H","F")
chisq.test(tab)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tab
## X-squared = 1.0355, df = 1, p-value = 0.3089
#p-value = 0.3089 à priori on ne rejette pas l'hypothèse d'indépendance mais il faut être prudent quand même, le test est moyennement puissant car effectif total=76
#Représentation graphique de disease/sexe avec nouvelles classes
par(mfrow=c(1,1))
barplot(t(prop.table(tab,margin=2)),beside=TRUE)
#Chez les hommes, 80% d'other, 20%d'AN, chez les femmes 64.3% d'Other, 35.7% d'AN
#Test Fisher exact (pour éviter de regrouper des classes entres elles)
fisher.test(table(disease,sex))
##
## Fisher's Exact Test for Count Data
##
## data: table(disease, sex)
## p-value = 0.2628
## alternative hypothesis: two.sided
#p-value=0.262 on ne rejette toujours l'hypothèse d'indépendance !
#(6) Représentation graphique de age/disease
boxplot(age~disease)
bartlett.test(age~disease)
##
## Bartlett test of homogeneity of variances
##
## data: age by disease
## Bartlett's K-squared = 7.962, df = 3, p-value = 0.0468
#p-value = 0.0468 variances inégales à priori mais attention effectifs très petits...
Anov=aov(age~disease)
#disease Residuals
#Sum of Squares 6151.874 10144.166
#Deg. of Freedom 3 72
summary(Anov)
## Df Sum Sq Mean Sq F value Pr(>F)
## disease 3 6152 2050.6 14.55 1.67e-07 ***
## Residuals 72 10144 140.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Df Sum Sq Mean Sq F value Pr(>F)
#disease 3 6151.9 2050.62 14.555 1.666e-07 ***
#Residuals 72 10144.2 140.89
#p-valeur très très petite, au moins un groupe très différent des autres
#Test comparaison deux à deux
? pairwise.t.test
pairwise.t.test(age,disease,p.adjust.method="bonferroni")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: age and disease
##
## Other GN AN
## GN 0.00052 - -
## AN 8.7e-07 1.00000 -
## PKD 8.9e-05 0.96273 1.00000
##
## P value adjustment method: bonferroni
#Avec correction de Bonferroni car multi-test ici
# Other GN AN
# GN 0.00052 - -
# AN 8.7e-07 1.00000 -
# PKD 8.9e-05 0.96273 1.00000
#Avec correction de Holm, même conclusion
pairwise.t.test(age,disease)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: age and disease
##
## Other GN AN
## GN 0.00035 - -
## AN 8.7e-07 0.48137 -
## PKD 7.4e-05 0.48137 0.57796
##
## P value adjustment method: holm
# Other GN AN
# GN 0.00035 - -
# AN 8.7e-07 0.48137 -
# PKD 7.4e-05 0.48137 0.57796
#P value adjustment method: holm
#Si on regroupe PKD et GN dans Other
ageOther=c(age[disease=="Other"],age[disease=="PKD"],age[disease=="GN"])
ageAN=age[disease=="AN"]
boxplot(ageOther,ageAN)
#Tests comparaison age/disease
shapiro.test(ageAN)
##
## Shapiro-Wilk normality test
##
## data: ageAN
## W = 0.94645, p-value = 0.2266
shapiro.test(ageOther)
##
## Shapiro-Wilk normality test
##
## data: ageOther
## W = 0.90574, p-value = 0.0005787
#normalité des âges n'est pas vérifié. Le t-test effectué est alors asymptotique. Attention car seulement 24 individus pour les AN!!
sd(ageAN)
## [1] 9.906812
sd(ageOther)
## [1] 15.39146
#On teste l'égalité des variances
bartlett.test(c(ageAN,ageOther)~c(rep(1,length(ageAN)),rep(0,length(ageOther))))
##
## Bartlett test of homogeneity of variances
##
## data: c(ageAN, ageOther) by c(rep(1, length(ageAN)), rep(0, length(ageOther)))
## Bartlett's K-squared = 5.3094, df = 1, p-value = 0.02121
#p-value = 0.02121 on rejette l'hypothèse d'ég. des variances avec un risque de 2.12% de chances de se tromper.
mean(ageAN)
## [1] 51.16667
mean(ageOther)
## [1] 40.25
t.test(ageAN,ageOther)
##
## Welch Two Sample t-test
##
## data: ageAN and ageOther
## t = 3.7128, df = 65.904, p-value = 0.0004238
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 5.046103 16.787230
## sample estimates:
## mean of x mean of y
## 51.16667 40.25000
#p-value = 0.0004238 Test très significatif, on rejette (H0) : l'âge moyen des patients AN est différent de celui de l'âge moyen des autres patients.
#Test non-paramétrique
ks.test(ageAN,ageOther)
## Warning in ks.test(ageAN, ageOther): cannot compute exact p-value with ties
##
## Two-sample Kolmogorov-Smirnov test
##
## data: ageAN and ageOther
## D = 0.38462, p-value = 0.01553
## alternative hypothesis: two-sided
#p-value = 0.01553 Test également significatif, on rejette l'hypothèse d'galité des lois entre ageAN et ageOther avec un risque de 1.553% de chances de se tromper
wilcox.test(ageAN,ageOther)
##
## Wilcoxon rank sum test with continuity correction
##
## data: ageAN and ageOther
## W = 853, p-value = 0.01057
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.01057 Idem !