Come condurre un'ancova in r

Questo tutorial fornisce un esempio di come eseguire un’ANCOVA in R.

Esempio: ANCOVA in R

Effettueremo un’ANCOVA per verificare se la tecnica di studio ha o meno un impatto sui risultati dell’esame utilizzando le seguenti variabili:

• Studio tecnico : la variabile indipendente che desideriamo analizzare
• Voto attuale dello studente : la covariata che desideriamo prendere in considerazione
• Punteggio revisione : le variabili di risposta che vogliamo analizzare

Il seguente set di dati contiene informazioni su 90 studenti divisi casualmente in tre gruppi da 30.

Il set di dati mostra la tecnica di studio utilizzata da ogni studente (A, B o C) , il voto attuale nella classe quando ha iniziato a utilizzare la tecnica di studio e il voto ricevuto all’esame dopo aver utilizzato la tecnica di studio per un mese per prepararsi per l’esame. esame:

 #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)

# technical 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

Passaggio 1: esplorare i dati

Prima di adattare il modello ANCOVA, dobbiamo prima esplorare i dati per comprenderli meglio e verificare che non vi siano valori anomali estremi che potrebbero distorcere i risultati.

Innanzitutto, possiamo visualizzare un riepilogo di ciascuna variabile nel set di dati:

 summary(data)

# technical 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  

Possiamo vedere che ciascun valore della tecnica di studio ( A, B e C) appare 30 volte nei dati.

Possiamo anche vedere come sono stati distribuiti i punteggi attuali degli studenti all’inizio dello studio. Il punteggio minimo della classe era 65,43, il massimo era 93,84 e la media era 78,15.

Allo stesso modo, possiamo vedere che il punteggio minimo ottenuto all’esame è stato 71,17, il punteggio massimo è stato 94,76 e la media è stata 83,38.

Quindi possiamo utilizzare il pacchetto dplyr per trovare facilmente la media e la deviazione standard dei voti attuali e dei risultati degli esami per ciascuna tecnica di studio:

 #load dplyr
library(dplyr)

data %>%
  group_by (technical) %>%
  summarize (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

Possiamo vedere che la media e le deviazioni standard del voto attuale degli studenti che utilizzano ciascuna tecnica di studio sono approssimativamente simili.

Possiamo anche vedere che il punteggio medio dell’esame è significativamente più alto per gli studenti che hanno utilizzato la tecnica di studio A rispetto alle tecniche B e C.

Possiamo anche visualizzare la distribuzione dei risultati degli esami in base alla tecnica di studio utilizzando i boxplot :

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

Allo stesso modo, possiamo anche utilizzare i boxplot per visualizzare la distribuzione dei voti attuali in base alla tecnica di studio:

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

Passaggio 2: verificare le ipotesi del modello

Dopo aver eseguito alcune esplorazioni di base dei dati e acquisito familiarità con i dati, dobbiamo verificare che siano soddisfatte le seguenti ipotesi per ANCOVA:

• La covariata e il trattamento sono indipendenti – è necessario verificare che la covariata ( voto attuale) e il trattamento (tecnica di studio) siano indipendenti l’uno dall’altro, perché aggiungere un termine covariata nel modello ha senso solo se la covariata e il trattamento agisce indipendentemente sulla variabile di risposta ( esame ).
• Omogeneità della varianza : dobbiamo verificare che le varianze tra i gruppi siano uguali

Per verificare che la covariata e il trattamento siano indipendenti, possiamo eseguire un’ANOVA utilizzando il grado corrente come variabile di risposta e la tecnica di studio come variabile predittrice:

 #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)
#technical 2 74 37.21 0.599 0.552
#Residuals 87 5406 62.14    

Il valore p è maggiore di 0,05, quindi la covariata ( grado attuale) e il trattamento ( tecnica di studio ) appaiono indipendenti.

Quindi, per verificare che vi sia omogeneità della varianza tra i gruppi, possiamo eseguire il test di Levene:

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

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

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

Il valore p del test è pari a 0,0001961, il che indica che le varianze tra i gruppi non sono uguali. Anche se potremmo tentare una trasformazione dei dati per correggere questo problema, per il momento non ci preoccuperemo troppo delle differenze di varianza.

Passaggio 3: regolare il modello ANCOVA

Successivamente, adatteremo il modello ANCOVA utilizzando il punteggio dell’esame come variabile di risposta, la tecnica di studio come variabile predittrice (o “trattamento”) e il voto attuale come covariata.

Utilizzeremo la funzione Anova() nel pacchetto perché per fare questo, giusto per poter specificare che vogliamo usare la somma dei quadrati di tipo III per il modello, perché la somma dei quadrati di tipo I dipende dall’ordine in cui i Predittori vengono inseriti nel modello:

 #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") 

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

Possiamo vedere che il valore p per la tecnica è estremamente basso, indicando che la tecnica di studio ha un effetto statisticamente significativo sui punteggi degli esami, anche dopo aver controllato il voto corrente.

Passaggio 4: test post-hoc

Sebbene i risultati dell’ANCOVA ci abbiano detto che la tecnica di studio ha avuto un effetto statisticamente significativo sui punteggi degli esami, dobbiamo condurre test post hoc per determinare quali tecniche di studio differiscono l’una dall’altra.

Per fare ciò, possiamo utilizzare la funzione glht() nel pacchetto multcomp in R per eseguire il test Tukey per confronti multipli:

 #load the multcomp library
library(multicomp)

#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 Assumptions:
#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 Assumptions:
# 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

Dal risultato possiamo vedere che esiste una differenza statisticamente significativa (a α = 0,05) nei risultati dell’esame tra lo studio della tecnica A e lo studio della tecnica B (valore p: 0,000109) nonché tra la tecnica A e la tecnica C (valore p: <1e-04).

Possiamo anche vedere che non vi è alcuna differenza statisticamente significativa (a α = 0,05) tra le tecniche B e C. Anche gli intervalli di confidenza tra le tecniche supportano questi risultati.

Pertanto, possiamo concludere che l’utilizzo della tecnica di studio A porta a un voto d’esame più alto in modo statisticamente significativo per gli studenti rispetto alle tecniche B e C , anche dopo aver controllato il voto attuale dello studente in classe.