Test des contrastes de Scheffé

plot of chunk
Lise_Vaudor_headband

Le test de Scheffé est un test qu’on applique souvent après une ANOVA: on parle de test post-hoc (au même titre, par exemple, qu’un test de Tukey). En effet, l’ANOVA à 1 facteur permet de mettre en évidence (le cas échéant) le fait qu’au moins un groupe a une moyenne différente des autres. Si on a affaire à 3 groupes ou plus, une question se pose alors: quels sont les groupes différents deux à deux?

Le test de Scheffé, en pratique:

plot of chunk
exemple_donnees_ANOVA_3_niveaux

Testons la significativité de l’effet de X sur Y à travers une ANOVA:

mymod=lm(Y~X)
anova(mymod)

## Analysis of Variance Table
## 
## Response: Y
##           Df Sum Sq Mean Sq F value Pr(>F)  
## X          2     61    30.7    2.55  0.087 .
## Residuals 59    710    12.0                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On rejette l’hypothèse d’égalité de toutes les moyennes au seuil de α = 10%.

Il n’existe pas à l’heure actuelle de test de Scheffé sous R.

_{^{(en fait, si, il y a une fonction scheffe.test dans un package qui s’appelle agricolae, mais il me semble après l’avoir testée qu’il y ait une petite erreur dans cette fonction, alors en attendant que son auteur la corrige -ça ne devrait pas être bien long- on va faire comme s’il n’y avait rien…).}}

Je vous propose d’utiliser une petite fonction que j’ai écrite (vous trouverez des explications sur le code dans le paragraphe suivant) et disponible ici.

Voici comment utiliser cette fonction:

source("test_Scheffe_maison.R")
Scheffe.test(X,Y,alpha=0.1)

##     contrast sigma    LCI    UCI    pval
## A-B   -2.555 1.133 -5.034 -0.076 0.08718
## A-C   -0.868 1.014 -3.087  1.352 0.69507
## B-C    1.688 1.161 -0.854  4.229 0.35448

Les résultats du test de Scheffé montrent que la différence de moyenne entre A et B est significative au seuil de α = 10%. En revanche on n’est pas en mesure d’affirmer que B est différent de C ou que A est différent de C.

Le test de Scheffé, en théorie

Maintenant que vous avez vu comment (et pourquoi) effectuer un test de Scheffé, revenons « en arrière » pour approfondir le principe de ce test.

On considère r groupes, et l’on définit μ_i comme étant la moyenne du groupe numéro i. Un contraste se définit de la manière suivante:

    $$C = \sum_{i=1}^r c_i\mu_i \\
      avec \quad \sum_{i=1}^r c_i = 0
    $$

Cela signifie que si l’on a trois groupes, on peut considérer (entre autres) les contrastes suivants:

    $$C_{A-B}=\mu_A-\mu_B \\
      C_{A-C}=\mu_B-\mu_C \\
      C_{B-C}=\mu_A-\mu_C
    $$

La valeur estimée d’un contraste est simplement:

    $$\hat{C} = \sum_{i=1}^r c_i\bar{Y}_i$$

Donc les estimations des trois contrastes qui nous intéressent sont:

estimated_mu=as.vector(tapply(Y,X,"mean"))
contrasts=c(estimated_mu[1]-estimated_mu[2],
            estimated_mu[1]-estimated_mu[3],
            estimated_mu[2]-estimated_mu[3])

print(contrasts)

## [1] -2.5551 -0.8675  1.6876

La variance associée à cet estimateur est:

    $$s_{\hat{C}}^2 = \hat{\sigma}^2\sum_{i=1}^r \frac{c_i^2}{n_i}$$

Donc l’écart-type associée à chacun des trois estimateurs considérés est:

n=as.vector(tapply(Y,X,"length"))
N=sum(n)
r=3
estimated_sigma=sqrt(sum(mymod$residuals^2)/(N-r))
sigma_contrasts=c(sqrt(estimated_sigma^2*(1/n[1]+1/n[2])),
                  sqrt(estimated_sigma^2*(1/n[1]+1/n[3])),
                  sqrt(estimated_sigma^2*(1/n[2]+1/n[3])))
print(sigma_contrasts)

## [1] 1.133 1.014 1.161

Ainsi il y a une probabilité de 1 − α que C se trouve dans l’intervalle suivant:

    $$\hat{C}\pm \quad s_\hat{C}\sqrt{\left(r-1\right)F_{\alpha;r-1;N-r}}$$

On peut donc calculer les intervalles de confiance comme suit:

alpha=0.1
delta=sigma_contrasts*sqrt(2*qf(1-alpha,3-1,length(X)-3))


results=data.frame(round(contrasts,3),
                   round(sigma_contrasts,3),
                   round(contrasts-delta,3),
                   round(contrasts+delta,3))
row.names(results)=c("A-B","A-C","B-C")
colnames(results)=c("contrast","sigma","LCI","UCI")

Par ailleurs la statistique suivante:

    $$T_{Scheffe}=\frac{C^2}{(r-1)s_{\hat{C}}^2}
    $$

suit une loi de Fisher sous hypothèse que tous les contrastes sont nuls.

On peut donc calculer les p-values comme suit:

T_Scheffe=((contrasts^2)/(sigma_contrasts^2))/(3-1)
pvalues=1-pf(T_Scheffe, 3-1,length(X)-3)
print(data.frame(results,pvalues))

##     contrast sigma    LCI    UCI pvalues
## A-B   -2.555 1.133 -5.034 -0.076 0.08718
## A-C   -0.868 1.014 -3.087  1.352 0.69507
## B-C    1.688 1.161 -0.854  4.229 0.35448

R-atique