Comment mener une ANCOVA dans R

Ce didacticiel fournit un exemple de la façon d’effectuer une ANCOVA dans R.

Exemple : ANCOVA dans R

Nous réaliserons une ANCOVA pour tester si la technique d’étude a ou non un impact sur les résultats des examens en utilisant les variables suivantes :

• Technique d’étude : La variable indépendante que nous souhaitons analyser
• Note actuelle de l’élève : La covariable que l’on souhaite prendre en compte
• Score d’examen : les variables de réponse que nous souhaitons analyser

L’ensemble de données suivant contient des informations sur 90 étudiants répartis au hasard en trois groupes de 30.

L’ensemble de données montre la technique d’étude utilisée par chaque élève (A, B ou C) , sa note actuelle dans la classe lorsqu’il a commencé à utiliser la technique d’étude et la note qu’il a obtenue à l’examen après avoir utilisé la technique d’étude pendant un mois pour se préparer à l’examen. examen:

#make this example reproducible 
set.seed(10)

#create dataset
data <- data.frame(technique = rep(c("A", "B", "C"), each = 30),
                   current_grade = runif(90, 65, 95),
                   exam = c(runif(30, 80, 95), runif(30, 70, 95), runif(30, 70, 90)))

#view first six lines of dataset
head(data)

#  technique current_grade     exam
#1         A      80.22435 87.32759
#2         A      74.20306 90.67114
#3         A      77.80723 88.87902
#4         A      85.79306 87.75735
#5         A      67.55408 85.72442
#6         A      71.76310 92.52167

Étape 1 : Explorer les données

Avant d’ajuster le modèle ANCOVA, nous devons d’abord explorer les données pour mieux les comprendre et vérifier qu’il n’existe pas de valeurs aberrantes extrêmes susceptibles de fausser les résultats.

Tout d’abord, nous pouvons afficher un résumé de chaque variable de l’ensemble de données :

summary(data)

# technique current_grade        exam      
# A:30      Min.   :65.43   Min.   :71.17  
# B:30      1st Qu.:71.79   1st Qu.:77.27  
# C:30      Median :77.84   Median :84.69  
#           Mean   :78.15   Mean   :83.38  
#           3rd Qu.:83.65   3rd Qu.:89.22  
#           Max.   :93.84   Max.   :94.76  

Nous pouvons voir que chaque valeur de la technique d’étude ( A, B et C) apparaît 30 fois dans les données.

Nous pouvons également voir comment les scores actuels des étudiants étaient répartis au début de l’étude. Le score minimum dans la classe était de 65,43, le maximum était de 93,84 et la moyenne était de 78,15.

De même, nous pouvons voir que la note minimale obtenue à l’examen était de 71,17, la note maximale était de 94,76 et la moyenne était de 83,38.

Ensuite, nous pouvons utiliser le package dplyr pour trouver facilement la moyenne et l’écart type des notes actuelles et des résultats des examens pour chaque technique d’étude :

#load dplyr
library(dplyr)

data %>%
  group_by(technique) %>%
  summarise(mean_grade = mean(current_grade),
            sd_grade = sd(current_grade),
            mean_exam = mean(exam),
            sd_exam = sd(exam))

# A tibble: 3 x 5
#  technique mean_grade sd_grade mean_exam sd_exam                      
#1 A               79.0     7.00      88.5    3.88
#2 B               78.5     8.33      81.8    7.62
#3 C               76.9     8.24      79.9    5.71

Nous pouvons voir que la moyenne et les écarts types de la note actuelle des étudiants utilisant chaque technique d’étude sont à peu près similaires.

Nous pouvons également constater que la note moyenne à l’examen est sensiblement plus élevée pour les étudiants qui ont utilisé la technique d’étude A par rapport aux techniques B et C.

Nous pouvons également visualiser la distribution des résultats des examens en fonction de la technique d’étude en utilisant des boxplots :

boxplot(exam ~ technique,
data = data,
main = "Exam Score by Studying Technique",
xlab = "Studying Technique",
ylab = "Exam Score",
col = "steelblue",
border = "black"
)

De même, nous pouvons également utiliser des boxplots pour visualiser la répartition des notes actuelles en fonction de la technique d’étude :

boxplot(current_grade ~ technique,
data = data,
main = "Current Grade by Studying Technique",
xlab = "Studying Technique",
ylab = "Current Grade",
col = "steelblue",
border = "black"
)

Étape 2 : Vérifiez les hypothèses du modèle

Une fois que nous avons effectué une exploration de base des données et que nous sommes familiers avec les données, nous devons vérifier que les hypothèses suivantes pour l’ANCOVA sont remplies :

• La covariable et le traitement sont indépendants – il faut vérifier que la covariable ( grade actuel) et le traitement (technique d’étude) sont indépendants l’un de l’autre, car l’ajout d’un terme de covariable dans le modèle n’a de sens que si la covariable et le traitement agissent indépendamment sur la variable de réponse ( examen ).
• Homogénéité de la variance – nous devons vérifier que les variances entre les groupes sont égales

Pour vérifier que la covariable et le traitement sont indépendants, nous pouvons effectuer une ANOVA en utilisant le grade actuel comme variable de réponse et la technique d’étude comme variable prédictive :

#fit anova model
anova_model <- aov(current_grade ~ technique, data = data)
#view summary of anova model
summary(anova_model)

#            Df Sum Sq Mean Sq F value Pr(>F)
#technique    2     74   37.21   0.599  0.552
#Residuals   87   5406   62.14    

La valeur p est supérieure à 0,05, donc la covariable ( grade actuel) et le traitement ( technique d’étude ) semblent indépendants.

Ensuite, pour vérifier qu’il existe une homogénéité de variance entre les groupes, nous pouvons effectuer le test de Levene :

#load car library to conduct Levene's Test
libary(car)

#conduct Levene's Test
leveneTest(exam~technique, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
#      Df F value    Pr(>F)    
#group  2  9.4324 0.0001961 ***
#      87   

La valeur p du test est égale à 0,0001961, ce qui indique que les variances entre les groupes ne sont pas égales. Même si nous pourrions tenter une transformation sur les données pour corriger ce problème, nous ne nous soucierons pas trop des différences de variance pour le moment.

Étape 3 : Ajuster le modèle ANCOVA

Ensuite, nous ajusterons le modèle ANCOVA en utilisant le score à l’examen comme variable de réponse, la technique d’étude comme variable prédictive (ou « traitement ») et la note actuelle comme covariable.

Nous utiliserons la fonction Anova() dans le package car pour ce faire, juste pour pouvoir spécifier que nous souhaitons utiliser la somme des carrés de type III pour le modèle, car la somme des carrés de type I dépend de l’ordre dans lequel les prédicteurs sont entrés dans le modèle :

#load car library
library(car)

#fit ANCOVA model
ancova_model <- aov(exam ~ technique + current_grade, data = data)

#view summary of model
Anova(ancova_model, type="III") 

#Response: exam
#              Sum Sq Df  F value    Pr(>F)    
#(Intercept)   7161.2  1 201.4621 < 2.2e-16 ***
#technique     1242.9  2  17.4830 4.255e-07 ***
#current_grade   12.3  1   0.3467    0.5576    
#Residuals     3057.0 86         

Nous pouvons voir que la valeur p pour la technique est extrêmement faible, ce qui indique que la technique d’étude a un effet statistiquement significatif sur les résultats aux examens, même après contrôle de la note actuelle.

Étape 4 : Tests post-hoc

Bien que les résultats de l’ANCOVA nous aient indiqué que la technique d’étude avait un effet statistiquement significatif sur les résultats aux examens, nous devons effectuer des tests post-hoc pour déterminer quelles techniques d’étude diffèrent les unes des autres.

Pour ce faire, nous pouvons utiliser la fonction glht() dans le package multcomp dans R pour effectuer le test de Tukey pour plusieurs comparaisons :

#load the multcomp library
library(multcomp)

#fit the ANCOVA model
ancova_model <- aov(exam ~ technique + current_grade, data = data)

#define the post hoc comparisons to make
postHocs <- glht(ancova_model, linfct = mcp(technique = "Tukey"))

#view a summary of the post hoc comparisons
summary(postHocs)

#Multiple Comparisons of Means: Tukey Contrasts
#
#Fit: aov(formula = exam ~ technique + current_grade, data = data)
#
#Linear Hypotheses:
#           Estimate Std. Error t value Pr(>|t|)    
#B - A == 0   -6.711      1.540  -4.358 0.000109 ***
#C - A == 0   -8.736      1.549  -5.640  < 1e-04 ***
#C - B == 0   -2.025      1.545  -1.311 0.393089    

#view the confidence intervals associated with the multiple comparisons
confint(postHocs)

#	 Simultaneous Confidence Intervals
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#Fit: aov(formula = exam ~ technique + current_grade, data = data)
#
#Quantile = 2.3845
#95% family-wise confidence level
#
#Linear Hypotheses:
#           Estimate lwr      upr     
#B - A == 0  -6.7112 -10.3832  -3.0392
#C - A == 0  -8.7364 -12.4302  -5.0426
#C - B == 0  -2.0252  -5.7091   1.6588

À partir du résultat, nous pouvons voir qu’il existe une différence statistiquement significative (à α = 0,05) dans les résultats des examens entre l’étude de la technique A et l’étude de la technique B (valeur p : 0,000109) ainsi qu’entre la technique A et la technique C ( valeur p : <1e-04).

Nous pouvons également constater qu’il n’y a pas de différence statistiquement significative (à α = 0,05) entre les techniques B et C . Les intervalles de confiance entre les techniques confirment également ces conclusions.

Ainsi, nous pouvons conclure que l’utilisation de la technique d’étude A conduit à une note d’examen statistiquement significativement plus élevée pour les étudiants par rapport aux techniques B et C , même après avoir contrôlé la note actuelle de l’étudiant dans la classe.