Jour 2 : Analyses statistiques de base
Packages utilisés :
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)
library(effsize)
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(ggplot2)
library(esquisse)1 Variables continues
1.1 Analyse descriptive
- On ouvre dans R studio la base de données obtenue hier en fin de séance (BD2)
BD2 <- read.csv("~/Nextcloud/Formation doctorale/BD2.csv")- On va principalement utiliser le package dplyr qui permet de simplifier les appels de commande et on va construire des tables descriptives à partir de ces commandes.
1.2 Fonction summarise_at
Les arguments de cette fonction sont : - La base de données - Un vecteur contenant les noms des variables dont on veut calculer les statistiques descriptives. - Une liste de fonctions (funs) contenant les fonctions statistiques à utiliser.
Premier exemple : statistiques descriptives des 3 temps de mesure de SGE :
stats_SGE1<-BD2 %>% summarise_at("SGE1",
lst(M=mean,
SD=sd,
MED=median),na.rm=TRUE)
stats_SGE2<-BD2%>% summarise_at(
"SGE2",
lst(M=mean,
SD=sd,
MED=median),na.rm=TRUE)
stats_SGE3<-BD2 %>% summarise_at(
"SGE3",
lst(M=mean,
SD=sd,
MED=median),na.rm=TRUE)
stats_SGE<-rbind(stats_SGE1,stats_SGE2,stats_SGE3)
row.names(stats_SGE)<-c("SGE1","SGE2","SGE3")
kable(stats_SGE,digits=3)| M | SD | MED | |
|---|---|---|---|
| SGE1 | 4.892 | 0.780 | 5.0 |
| SGE2 | 4.638 | 0.868 | 4.8 |
| SGE3 | 4.529 | 0.774 | 4.6 |
Attention j’ai le droit de mettre na.rm=TRUE à l’extérieur des fonctions car cet argument correspond à toutes les fonctions
1.3 Ajout d’une fonction personnalisée
On va ajouter l’intervalle de confiance de la moyenne par bootstrap paramétrique
bmeanCI<-function(x){
x<-na.omit(x)
bmean<-NULL
for(b in 1:1000){
bmean<-c(bmean,mean(sample(x,length(x),replace=TRUE)))
}
return(quantile(bmean,probs=c(.025,.975)))
}On ajoute cette fonction aux statistiques descriptives précédentes :
stats_SGE1<-BD2%>% summarise_at(
"SGE1",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)) %>% setNames(c("n","M","SD","ciLower","ciUpper"))
stats_SGE1## n M SD ciLower ciUpper
## 1 295 4.891636 0.7795066 4.797036 4.9876551.4 Fonction group_by
Si on veut obtenir comme précédemment n,M,SD,ciLower,ciUpper de SGE1 regroupé par sexe :
SGE1_by_Sexe<-BD2%>%group_by(Sexe)%>%
summarise_at("SGE1",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Sexe","n","M","SD","ciLower","ciUpper"))
kable(SGE1_by_Sexe,digits=3) | Sexe | n | M | SD | ciLower | ciUpper |
|---|---|---|---|---|---|
| F | 163 | 4.940 | 0.736 | 4.831 | 5.055 |
| M | 132 | 4.829 | 0.832 | 4.684 | 4.971 |
Si on veut obtenir comme précédemment n,M,SD,ciLower,ciUpper de SGE1 regroupé par sexe puis par Strate :
SGE1_by_Sexe_strate<-BD2%>%group_by(Sexe,T1Strate)%>% summarise_at("SGE1",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Sexe","Strate","n","M","SD","ciLower","ciUpper"))
kable(SGE1_by_Sexe_strate,digits=3) | Sexe | Strate | n | M | SD | ciLower | ciUpper |
|---|---|---|---|---|---|---|
| F | 1 | 91 | 4.932 | 0.730 | 4.768 | 5.073 |
| F | 2 | 36 | 4.976 | 0.761 | 4.723 | 5.212 |
| F | 3 | 36 | 4.924 | 0.748 | 4.665 | 5.171 |
| M | 1 | 81 | 4.693 | 0.957 | 4.468 | 4.904 |
| M | 2 | 24 | 5.064 | 0.491 | 4.864 | 5.264 |
| M | 3 | 27 | 5.025 | 0.560 | 4.817 | 5.233 |
1.5 Exercice
Pour SGE1,2,3 en fonction de la variable T1StrateF et AgePass1, faire la table ci-dessous :
SGE1_by_Strate_AgePass<-BD2%>%group_by(T1StrateF,AgePass1)%>% summarise_at("SGE1",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Strate","Age","n","M","SD","ciLower","ciUpper"))
SGE2_by_Strate_AgePass<-BD2%>%group_by(T1StrateF,AgePass1)%>% summarise_at("SGE2",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Strate","Age","n","M","SD","ciLower","ciUpper"))
SGE3_by_Strate_AgePass<-BD2%>%group_by(T1StrateF,AgePass1)%>% summarise_at("SGE3",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Strate","Age","n","M","SD","ciLower","ciUpper"))
tab<-as.data.frame(rbind(SGE1_by_Strate_AgePass,
SGE2_by_Strate_AgePass,
SGE3_by_Strate_AgePass))
names<-c(rep("SGE1",8),rep("SGE2",8),rep("SGE3",8))
tab_f<-cbind(names,tab)
kable(tab_f,digits=3) | names | Strate | Age | n | M | SD | ciLower | ciUpper |
|---|---|---|---|---|---|---|---|
| SGE1 | PRIVE | De 10 à 12 ans | 34 | 5.042 | 0.610 | 4.836 | 5.237 |
| SGE1 | PRIVE | Moins de 10 ans | 17 | 5.100 | 0.800 | 4.650 | 5.463 |
| SGE1 | PRIVE | Plus de 12 ans | 9 | 4.657 | 0.538 | 4.314 | 5.057 |
| SGE1 | PUBLIC | De 10 à 12 ans | 55 | 4.764 | 0.864 | 4.542 | 4.964 |
| SGE1 | PUBLIC | Moins de 10 ans | 84 | 5.042 | 0.745 | 4.894 | 5.205 |
| SGE1 | PUBLIC | Plus de 12 ans | 33 | 4.359 | 0.884 | 4.069 | 4.669 |
| SGE1 | ZEP | De 10 à 12 ans | 41 | 4.780 | 0.619 | 4.600 | 4.956 |
| SGE1 | ZEP | Moins de 10 ans | 22 | 5.412 | 0.598 | 5.106 | 5.647 |
| SGE2 | PRIVE | De 10 à 12 ans | 34 | 4.641 | 0.649 | 4.429 | 4.847 |
| SGE2 | PRIVE | Moins de 10 ans | 17 | 5.024 | 0.681 | 4.694 | 5.341 |
| SGE2 | PRIVE | Plus de 12 ans | 9 | 3.978 | 0.724 | 3.533 | 4.422 |
| SGE2 | PUBLIC | De 10 à 12 ans | 55 | 4.447 | 0.922 | 4.189 | 4.669 |
| SGE2 | PUBLIC | Moins de 10 ans | 84 | 4.952 | 0.793 | 4.774 | 5.117 |
| SGE2 | PUBLIC | Plus de 12 ans | 33 | 4.133 | 0.729 | 3.879 | 4.370 |
| SGE2 | ZEP | De 10 à 12 ans | 41 | 4.556 | 0.966 | 4.273 | 4.810 |
| SGE2 | ZEP | Moins de 10 ans | 22 | 4.791 | 0.910 | 4.409 | 5.155 |
| SGE3 | PRIVE | De 10 à 12 ans | 34 | 4.553 | 0.623 | 4.329 | 4.753 |
| SGE3 | PRIVE | Moins de 10 ans | 17 | 4.718 | 0.889 | 4.282 | 5.094 |
| SGE3 | PRIVE | Plus de 12 ans | 9 | 4.133 | 0.819 | 3.599 | 4.600 |
| SGE3 | PUBLIC | De 10 à 12 ans | 55 | 4.367 | 0.784 | 4.171 | 4.578 |
| SGE3 | PUBLIC | Moins de 10 ans | 84 | 4.786 | 0.717 | 4.631 | 4.933 |
| SGE3 | PUBLIC | Plus de 12 ans | 33 | 4.291 | 0.781 | 4.024 | 4.552 |
| SGE3 | ZEP | De 10 à 12 ans | 41 | 4.488 | 0.662 | 4.297 | 4.678 |
| SGE3 | ZEP | Moins de 10 ans | 22 | 4.373 | 0.996 | 3.964 | 4.800 |
1.6 Solution
SGE1_by_Strate_AgePass<-BD2%>%group_by(T1StrateF,AgePass1)%>% summarise_at("SGE1",
lst(length,
~mean(.,na.rm=TRUE),
~sd(.,na.rm=TRUE),
~bmeanCI(.)[1],
~bmeanCI(.)[2]
)
) %>% setNames(c("Strate","Age","n","M","SD","ciLower","ciUpper"))
## Idem pour SGE2_by_Strate_AgePass et SGE3_by_Strate_AgePass
tab<-as.data.frame(rbind(SGE1_by_Strate_AgePass,
SGE2_by_Strate_AgePass,
SGE3_by_Strate_AgePass))
names<-c(rep("SGE1",8),rep("SGE2",8),rep("SGE3",8))
tab_f<-cbind(names,tab)
kable(tab_f,digits=3) | names | Strate | Age | n | M | SD | ciLower | ciUpper |
|---|---|---|---|---|---|---|---|
| SGE1 | PRIVE | De 10 à 12 ans | 34 | 5.042 | 0.610 | 4.830 | 5.255 |
| SGE1 | PRIVE | Moins de 10 ans | 17 | 5.100 | 0.800 | 4.700 | 5.463 |
| SGE1 | PRIVE | Plus de 12 ans | 9 | 4.657 | 0.538 | 4.286 | 5.057 |
| SGE1 | PUBLIC | De 10 à 12 ans | 55 | 4.764 | 0.864 | 4.513 | 4.982 |
| SGE1 | PUBLIC | Moins de 10 ans | 84 | 5.042 | 0.745 | 4.868 | 5.205 |
| SGE1 | PUBLIC | Plus de 12 ans | 33 | 4.359 | 0.884 | 4.062 | 4.669 |
| SGE1 | ZEP | De 10 à 12 ans | 41 | 4.780 | 0.619 | 4.590 | 4.966 |
| SGE1 | ZEP | Moins de 10 ans | 22 | 5.412 | 0.598 | 5.106 | 5.659 |
| SGE2 | PRIVE | De 10 à 12 ans | 34 | 4.641 | 0.649 | 4.429 | 4.847 |
| SGE2 | PRIVE | Moins de 10 ans | 17 | 5.024 | 0.681 | 4.694 | 5.341 |
| SGE2 | PRIVE | Plus de 12 ans | 9 | 3.978 | 0.724 | 3.533 | 4.422 |
| SGE2 | PUBLIC | De 10 à 12 ans | 55 | 4.447 | 0.922 | 4.189 | 4.669 |
| SGE2 | PUBLIC | Moins de 10 ans | 84 | 4.952 | 0.793 | 4.774 | 5.117 |
| SGE2 | PUBLIC | Plus de 12 ans | 33 | 4.133 | 0.729 | 3.879 | 4.370 |
| SGE2 | ZEP | De 10 à 12 ans | 41 | 4.556 | 0.966 | 4.273 | 4.810 |
| SGE2 | ZEP | Moins de 10 ans | 22 | 4.791 | 0.910 | 4.409 | 5.155 |
| SGE3 | PRIVE | De 10 à 12 ans | 34 | 4.553 | 0.623 | 4.329 | 4.753 |
| SGE3 | PRIVE | Moins de 10 ans | 17 | 4.718 | 0.889 | 4.282 | 5.094 |
| SGE3 | PRIVE | Plus de 12 ans | 9 | 4.133 | 0.819 | 3.599 | 4.600 |
| SGE3 | PUBLIC | De 10 à 12 ans | 55 | 4.367 | 0.784 | 4.171 | 4.578 |
| SGE3 | PUBLIC | Moins de 10 ans | 84 | 4.786 | 0.717 | 4.631 | 4.933 |
| SGE3 | PUBLIC | Plus de 12 ans | 33 | 4.291 | 0.781 | 4.024 | 4.552 |
| SGE3 | ZEP | De 10 à 12 ans | 41 | 4.488 | 0.662 | 4.297 | 4.678 |
| SGE3 | ZEP | Moins de 10 ans | 22 | 4.373 | 0.996 | 3.964 | 4.800 |
2 Variables catégorielles
2.1 Une seule variable
On peut calculer les effectifs de chaque modalité d’une variable catégorielle
BD2 %>%
group_by(AgePass1) %>%
summarise(effectif=n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 3 x 2
## AgePass1 effectif
## <chr> <int>
## 1 De 10 à 12 ans 130
## 2 Moins de 10 ans 123
## 3 Plus de 12 ans 42Mais aussi ajouter une colonne % :
BD2 %>%
group_by(AgePass1) %>%
summarise(effectif=n()) %>%
mutate(percent=100*effectif/sum(effectif))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 3 x 3
## AgePass1 effectif percent
## <chr> <int> <dbl>
## 1 De 10 à 12 ans 130 44.1
## 2 Moins de 10 ans 123 41.7
## 3 Plus de 12 ans 42 14.22.2 Deux variables catégorielles
On compte les effectifs du croisement des deux variables :
BD2 %>%
group_by(Sexe,AgePass1) %>%
summarise(effectif=n())## `summarise()` regrouping output by 'Sexe' (override with `.groups` argument)## # A tibble: 6 x 3
## # Groups: Sexe [2]
## Sexe AgePass1 effectif
## <chr> <chr> <int>
## 1 F De 10 à 12 ans 68
## 2 F Moins de 10 ans 76
## 3 F Plus de 12 ans 19
## 4 M De 10 à 12 ans 62
## 5 M Moins de 10 ans 47
## 6 M Plus de 12 ans 23On peut comme précédemment calculer une colonne % en fonction du sexe:
BD2 %>%
group_by(Sexe,AgePass1) %>%
summarise(effectif=n())%>%
mutate(percent=100*effectif/sum(effectif))## `summarise()` regrouping output by 'Sexe' (override with `.groups` argument)## # A tibble: 6 x 4
## # Groups: Sexe [2]
## Sexe AgePass1 effectif percent
## <chr> <chr> <int> <dbl>
## 1 F De 10 à 12 ans 68 41.7
## 2 F Moins de 10 ans 76 46.6
## 3 F Plus de 12 ans 19 11.7
## 4 M De 10 à 12 ans 62 47.0
## 5 M Moins de 10 ans 47 35.6
## 6 M Plus de 12 ans 23 17.4Ou bien avec la fonction ungroup() pour calculer les % sur l’effectif total :
BD2 %>%
group_by(Sexe,AgePass1) %>%
summarise(effectif=n())%>%
ungroup() %>%
mutate(percent=100*effectif/sum(effectif))## `summarise()` regrouping output by 'Sexe' (override with `.groups` argument)## # A tibble: 6 x 4
## Sexe AgePass1 effectif percent
## <chr> <chr> <int> <dbl>
## 1 F De 10 à 12 ans 68 23.1
## 2 F Moins de 10 ans 76 25.8
## 3 F Plus de 12 ans 19 6.44
## 4 M De 10 à 12 ans 62 21.0
## 5 M Moins de 10 ans 47 15.9
## 6 M Plus de 12 ans 23 7.803 Test d’hypothèse
3.1 Exemple
On considère \(X\mathcal N (\mu,1)\) et on veut tester \(\mathcal H_0 : \mu=0\) contre \(\mathcal H_1 : \mu \not =0\).
L’approche fréquentiste est basée sur le calcul de la \(p-\)value
3.2 Distribution de la p-value
On crée la fonction:
boot.pval<-function(mu,n){
B<-10000
pval<-NULL
for(b in 1:B){
X<-rnorm(n,mu,1)
pval<-c(pval,t.test(X)$p.value)
}
return(pval)
}3.3 Distribution de la p-value sous \(\mathcal H_0\)
proba_rejet<-matrix(0,nrow=4,ncol=4)
mu<-0
essai1<-boot.pval(mu,30)
essai2<-boot.pval(mu,50)
essai3<-boot.pval(mu,100)
essai4<-boot.pval(mu,1000)
proba_rejet[1,]<-apply(cbind(essai1,essai2,essai3,essai4),2,function(x) mean(x<.05))
par(mfrow=c(1,4))
hist(essai1,breaks = seq(0,1,by=.1))
hist(essai2,breaks = seq(0,1,by=.1))
hist(essai3,breaks = seq(0,1,by=.1))
hist(essai4,breaks = seq(0,1,by=.1))
3.4 Distribution de la p-value sous \(\mathcal H_1\)
mu<-0.2
essai1<-boot.pval(mu,30)
essai2<-boot.pval(mu,50)
essai3<-boot.pval(mu,100)
essai4<-boot.pval(mu,1000)
proba_rejet[2,]<-apply(cbind(essai1,essai2,essai3,essai4),2,function(x) mean(x<.05))
par(mfrow=c(3,4))
hist(essai1,breaks = seq(0,1,by=.1))
hist(essai2,breaks = seq(0,1,by=.1))
hist(essai3,breaks = seq(0,1,by=.1))
hist(essai4,breaks = seq(0,1,by=.1))
mu<-0.5
essai1<-boot.pval(mu,30)
essai2<-boot.pval(mu,50)
essai3<-boot.pval(mu,100)
essai4<-boot.pval(mu,1000)
proba_rejet[3,]<-apply(cbind(essai1,essai2,essai3,essai4),2,function(x) mean(x<.05))
hist(essai1,breaks = seq(0,1,by=.1))
hist(essai2,breaks = seq(0,1,by=.1))
hist(essai3,breaks = seq(0,1,by=.1))
hist(essai4,breaks = seq(0,1,by=.1))
mu<-0.8
essai1<-boot.pval(mu,30)
essai2<-boot.pval(mu,50)
essai3<-boot.pval(mu,100)
essai4<-boot.pval(mu,1000)
proba_rejet[4,]<-apply(cbind(essai1,essai2,essai3,essai4),2,function(x) mean(x<.05))
hist(essai1,breaks = seq(0,1,by=.1))
hist(essai2,breaks = seq(0,1,by=.1))
hist(essai3,breaks = seq(0,1,by=.1))
hist(essai4,breaks = seq(0,1,by=.1))
3.5 Table de puissance
colnames(proba_rejet)<-c("n=30","n=50","n=100","n=1000")
rownames(proba_rejet)<-c("H_0","H_1 avec mu=.2","H_1 avec mu=.5","H_1 avec mu=.8")
kable(proba_rejet,digits=3)| n=30 | n=50 | n=100 | n=1000 | |
|---|---|---|---|---|
| H_0 | 0.051 | 0.050 | 0.053 | 0.052 |
| H_1 avec mu=.2 | 0.185 | 0.283 | 0.516 | 1.000 |
| H_1 avec mu=.5 | 0.759 | 0.934 | 0.999 | 1.000 |
| H_1 avec mu=.8 | 0.988 | 1.000 | 1.000 | 1.000 |
4 Test de Student
4.1 Comparaisons sur deux groupes indépendants
Il existe deux types de tests de Student (var.equal=TRUE or FALSE): - Les deux échantillons ont des variances peu différentes (test de Student usuel) - Les deux échantillons ont des variances différentes (test avec la correction de Welch).
var.test(SGE1~Sexe,data=BD2)##
## F test to compare two variances
##
## data: SGE1 by Sexe
## F = 0.78226, num df = 155, denom df = 118, p-value = 0.152
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5545285 1.0947988
## sample estimates:
## ratio of variances
## 0.7822622t.test(SGE1~Sexe,data=BD2)##
## Welch Two Sample t-test
##
## data: SGE1 by Sexe
## t = 1.1534, df = 236.68, p-value = 0.2499
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0787142 0.3010585
## sample estimates:
## mean in group F mean in group M
## 4.939744 4.828571cohen.d(SGE1~Sexe,data=BD2)##
## Cohen's d
##
## d estimate: 0.1427162 (negligible)
## 95 percent confidence interval:
## lower upper
## -0.09719497 0.382627324.1.1 Construire une table aux normes APA
test.Sudent.Ind<-function(x,group){
group<-factor(group)
descript<-NULL
names_desc<-NULL
for(i in 1:length(levels(group))){
descript<-c(descript,c(mean(x[group==levels(group)[i]],na.rm=TRUE),sd(x[group==levels(group)[i]],na.rm=TRUE)))
names_desc<-c(names_desc,c(paste("Mean_",levels(group)[i]),paste("SD_",levels(group)[i])))
}
if(var.test(x~group)$p.value>.05){
test<-t.test(x~group,var.equal=TRUE)
}else{
test<-t.test(x~group,var.equal=FALSE)
}
effect<-cohen.d(x~group)
res<-c(descript,test$statistic,
test$parameter,
test$p.value,
effect$estimate,
effect$conf.int)
names(res)<-c(names_desc,"t","df","p","d","CI_lower_d","CI_upper_d")
return(res)
}4.1.2 Application
On construit une table aux normes APA qui donne les t.test de SGE1,SGE2,SGE3 entre les filles et les garçons.
tab<-rbind(test.Sudent.Ind(x=BD2$SGE1,group=BD2$Sexe),
test.Sudent.Ind(x=BD2$SGE2,group=BD2$Sexe),
test.Sudent.Ind(x=BD2$SGE3,group=BD2$Sexe))
rownames(tab)<-c("SGE1","SGE2","SGE3")
kable(tab,digits=3)| Mean_ F | SD_ F | Mean_ M | SD_ M | t | df | p | d | CI_lower_d | CI_upper_d | |
|---|---|---|---|---|---|---|---|---|---|---|
| SGE1 | 4.940 | 0.736 | 4.829 | 0.832 | 1.173 | 273 | 0.242 | 0.143 | -0.097 | 0.383 |
| SGE2 | 4.714 | 0.806 | 4.544 | 0.932 | 1.680 | 293 | 0.094 | 0.197 | -0.034 | 0.428 |
| SGE3 | 4.553 | 0.767 | 4.500 | 0.784 | 0.588 | 293 | 0.557 | 0.069 | -0.162 | 0.299 |
4.2 Comparaison de deux temps de mesure
Une même variable quantitative est mesurée deux fois sur les mêmes individus.
Ce test est aussi appelé test de Student sur groupes appariés.
test.Student.App<-function(x1,x2){
x<-na.omit(cbind(x1,x2))
fact<-c(rep("0",dim(x)[1]),rep("1",dim(x)[1]))
X<-c(x[,1],x[,2])
#subj<-factor(rep(1:dim(x)[1],2))
test<-t.test(X~fact,paired = TRUE)
effect<-cohen.d(X~fact,paired = TRUE)
res<-c(mean(x[,1]),sd(x[,1]),mean(x[,2]),sd(x[,2]),test$statistic,test$parameter,test$p.value,effect$estimate,effect$conf.int)
names(res)<-c("M1","SD1","M2","SD2","t","df","p","d","Lower_CI_d","Upper_CI_d")
return(res)
}4.3 Exemple
kable(t(test.Student.App(BD2$SGE1,BD2$SGE2)),digits=3)| M1 | SD1 | M2 | SD2 | t | df | p | d | Lower_CI_d | Upper_CI_d |
|---|---|---|---|---|---|---|---|---|---|
| 4.892 | 0.78 | 4.631 | 0.877 | 5.247 | 274 | 0 | 0.313 | 0.193 | 0.433 |
4.4 Comparaison de plus de 2 moyennes
4.4.1 Anova à 1 facteur
BD2$T2Classe<-factor(BD2$T2Classe,levels=c("CM1","CM2","6ieme","5ieme","4ieme"))
table(BD2$T2Classe)##
## CM1 CM2 6ieme 5ieme 4ieme
## 84 51 71 66 23aov1<-aov(RE2~T2Classe,data=BD2)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## T2Classe 4 31.84 7.960 21.47 1.59e-15 ***
## Residuals 290 107.53 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.4.2 Tests Post Hoc
tk1<-glht(aov1, linfct = mcp(T2Classe = "Tukey"))
summary(tk1)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = RE2 ~ T2Classe, data = BD2)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## CM2 - CM1 == 0 0.19734 0.10810 1.826 0.35382
## 6ieme - CM1 == 0 -0.49943 0.09817 -5.088 < 1e-04 ***
## 5ieme - CM1 == 0 -0.63831 0.10016 -6.373 < 1e-04 ***
## 4ieme - CM1 == 0 -0.55663 0.14330 -3.884 0.00116 **
## 6ieme - CM2 == 0 -0.69677 0.11177 -6.234 < 1e-04 ***
## 5ieme - CM2 == 0 -0.83565 0.11353 -7.361 < 1e-04 ***
## 4ieme - CM2 == 0 -0.75396 0.15295 -4.930 < 1e-04 ***
## 5ieme - 6ieme == 0 -0.13888 0.10412 -1.334 0.66393
## 4ieme - 6ieme == 0 -0.05720 0.14610 -0.391 0.99485
## 4ieme - 5ieme == 0 0.08169 0.14744 0.554 0.98080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)cld(tk1)## CM1 CM2 6ieme 5ieme 4ieme
## "b" "b" "a" "a" "a"4.5 Normes APA
test.postHoc.Ind<-function(x,group){
group<-factor(group)
descript<-NULL
names_desc<-NULL
aov1<-aov(x~group)
tk1<-glht(aov1, linfct = mcp(group = "Tukey"))
res1<-cld(tk1)$mcletters$Letters
for(i in 1:length(levels(group))){
descript<-c(descript,round(c(mean(x[group==levels(group)[i]],na.rm=TRUE),sd(x[group==levels(group)[i]],na.rm=TRUE)),3),res1[i])
names_desc<-c(names_desc,c(paste("Mean_",levels(group)[i]),paste("SD_",levels(group)[i]),"GrPost"))
}
names(descript)<-names_desc
return(descript)
}4.6 Exemple
tab<-rbind(test.postHoc.Ind(BD2$SGE2,BD2$T2Classe),
test.postHoc.Ind(BD2$RE2,BD2$T2Classe),
test.postHoc.Ind(BD2$Sec2,BD2$T2Classe))
row.names(tab)<-c("SGE","RE","Sec")
kable(tab)| Mean_ CM1 | SD_ CM1 | GrPost | Mean_ CM2 | SD_ CM2 | GrPost | Mean_ 6ieme | SD_ 6ieme | GrPost | Mean_ 5ieme | SD_ 5ieme | GrPost | Mean_ 4ieme | SD_ 4ieme | GrPost | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| SGE | 5.024 | 0.779 | b | 4.871 | 0.742 | b | 4.411 | 1.016 | a | 4.448 | 0.693 | a | 3.957 | 0.603 | a |
| RE | 3.026 | 0.648 | b | 3.224 | 0.548 | b | 2.527 | 0.628 | a | 2.388 | 0.61 | a | 2.47 | 0.517 | a |
| Sec | 3.01 | 0.71 | a | 3.071 | 0.706 | ab | 3.118 | 0.67 | ab | 3.448 | 0.462 | c | 3.452 | 0.387 | bc |
5 Tests du CHI2
Ces tests s’appliquent à une ou des variables catégorielles.
5.1 Indépendance
table(BD2$Sexe,BD2$T1StrateF)##
## PRIVE PUBLIC ZEP
## F 36 91 36
## M 24 81 27res<-chisq.test(BD2$Sexe,BD2$T1StrateF)
res##
## Pearson's Chi-squared test
##
## data: BD2$Sexe and BD2$T1StrateF
## X-squared = 1.0208, df = 2, p-value = 0.6003kable(res$expected)| PRIVE | PUBLIC | ZEP | |
|---|---|---|---|
| F | 33.15254 | 95.03729 | 34.81017 |
| M | 26.84746 | 76.96271 | 28.18983 |
kable(res$residuals)| PRIVE | PUBLIC | ZEP | |
|---|---|---|---|
| F | 0.4945371 | -0.4141358 | 0.2016657 |
| M | -0.5495481 | 0.4602031 | -0.2240985 |
Taille d’effet : il s’agit du \(V\) de Cramer :
CramerV = function(x,y) {
V = sqrt(chisq.test(x, y, correct=FALSE)$statistic /
(length(x) * (min(length(unique(x)),length(unique(y))) - 1)))
return(as.numeric(V))
}
CramerV(BD2$Sexe,BD2$T1StrateF)## [1] 0.058823345.2 Adéquation
x<-table(BD2$T1StrateF)
res<-chisq.test(x,p=c(.4,.4,.2))
res##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 53.492, df = 2, p-value = 2.424e-12res$expected## PRIVE PUBLIC ZEP
## 118 118 59