Tests non paramétriques

Le diaporama sur les tests non paramétriques :

Slides

Indice de dispersion non paramétrique

On définit le MAD (Median Absolute Deviation) qui peut être vu comme un équivalent non paramétrique à \(\sigma\).

\[ MAD=\text{median}|X_i-\widetilde X|. \] Donc pour \(X\sim \mathcal N(\mu,\sigma)\) on a par définition \(P(|X-\mu|<MAD)=0.5\) ce qui donne \(P(|Z|<\frac{MAD}{\sigma})=0.5\) où \(Z\sim \mathcal N(0,1)\) et donc \(\frac{MAD}{\sigma}=\Phi^{-1}(\frac{3}{4}).\)

Comparaison de deux groupes indépendants

On commence par l’exemple

library(coin)

Loading required package: survival

data("neuropathy")
library(ggplot2)
library(dplyr)


Attaching package: 'dplyr'

The following object is masked from 'package:kableExtra':

    group_rows

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

On représente graphiquement les données

ggplot(neuropathy,aes(x=group,y=pain))+
  geom_boxplot()+
  labs(title="Distribution des résultats selon le groupe",
       y="Mesure de la douleur",
       x="Groupe"
       )

On regarde aussi les paramètres de position et de dispersion

neuropathy |>
  group_by(group)|>
  summarise(M=mean(pain),SD=sd(pain),
            Med=median(pain),mad=mad(pain))

# A tibble: 2 × 5
  group       M    SD    Med   mad
  <fct>   <dbl> <dbl>  <dbl> <dbl>
1 control 0.200 0.515 0.0925 0.205
2 treat   0.430 0.788 0.280  0.523

Définition des deux séries de données :

x=neuropathy |>
  filter(group=="control") |> 
  select(pain) |>
  unlist() # Pour en faire un vecteur

Faire la même chose pour créer le vecteur y

Show the code

y=neuropathy |>
  filter(group=="treat") |> 
  select(pain) |>
  unlist()

Les tests sur les paramètres de position

On regarde tout d’abord les tests paramétriques :

# Student
t.test(x,y,var.equal=T)


    Two Sample t-test

data:  x and y
t = -1.3279, df = 56, p-value = 0.1896
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.5784146  0.1172622
sample estimates:
mean of x mean of y 
0.1995667 0.4301429

# Welch
t.test(x,y,var.equal=F)


    Welch Two Sample t-test

data:  x and y
t = -1.3094, df = 46.051, p-value = 0.1969
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.5850265  0.1238742
sample estimates:
mean of x mean of y 
0.1995667 0.4301429

# test de Yuen
library(robnptests)

Puis les tests non paramétriques :

# Test de Mann-Withney
wilcox.test(x,y)

Warning in wilcox.test.default(x, y): cannot compute exact p-value with ties


    Wilcoxon rank sum test with continuity correction

data:  x and y
W = 357, p-value = 0.3301
alternative hypothesis: true location shift is not equal to 0

trimmed_test(x,y)


    Randomization test based on trimmed means (1000 random permutations)

data:  x and y
trimmed t = -1.3775, df = 34, p-value = 0.1588
alternative hypothesis: true location shift is not equal to 0
sample estimates:
Trimmed mean of x Trimmed mean of y 
        0.1320000         0.3602778

# Test de Behrens-Fisher
library(npsm)

Loading required package: Rfit

fp.test(x,y)

statistic =  0.9124277 , p-value =  0.3653668

# Test de permutation
library(coin)
oneway_test(pain~group,data=neuropathy)


    Asymptotic Two-Sample Fisher-Pitman Permutation Test

data:  pain by group (control, treat)
Z = -1.3191, p-value = 0.1871
alternative hypothesis: true mu is not equal to 0

# Test de Brunner-Manzel
library(brunnermunzel)
brunnermunzel.test(x,y)


    Brunner-Munzel Test

data:  x and y
Brunner-Munzel Test Statistic = 0.94424, df = 41.597, p-value = 0.3505
95 percent confidence interval:
 0.4146606 0.7353394
sample estimates:
P(X<Y)+.5*P(X=Y) 
           0.575

Les tests sur les paramètres d’échelle

# Test de Levene
library(car)

Loading required package: carData


Attaching package: 'car'

The following object is masked from 'package:Rfit':

    subsets

The following object is masked from 'package:dplyr':

    recode

leveneTest(pain~group,data=neuropathy)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value  Pr(>F)  
group  1  4.4733 0.03889 *
      56                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Test de Fligner-Killeen
fk.test(x,y)

statistic =  2.193009 , p-value =  0.02830674 
95  percent confidence interval:
0.8423731 4.273328 
Estimate: 1.897297

Répondre aux questions suivantes :

Que pensez-vous de ces résultats ?
Quel test vous paraît le plus adapté compte tenu des données considérées ?

Un autre exemple :

data("mercuryfish")

On s’interesse aux trois questions suivantes :

Le taux de mercure présent dans le sang est-il significativement supérieur chez le groupe exposé ?
La proportion de cellules avec des anomalies structurales est-elle significativement supérieure chez le groupe exposé ?

Comparaison de deux groupes appariés

Charger les data de ObrienKaiserLong. On va uniquement dans un premier temps s’intéresser à la comparaison pre versus post :

data("OBrienKaiserLong")
data=OBrienKaiserLong |>
  filter(phase %in% c("pre","post")) |>
  select(score,phase)

Peut-on penser que les scores ont significativement évolués entre avant et après ?

Show the code

data |>
  group_by(phase) |>
  summarise(M=mean(score),SD=sd(score),
            Med=median(score),mad=mad(score))

# A tibble: 2 × 5
  phase     M    SD   Med   mad
  <fct> <dbl> <dbl> <dbl> <dbl>
1 pre    4.38  1.87     4  1.48
2 post   5.75  2.33     6  2.97

Show the code

ggplot(data,aes(x=phase,y=score))+
  geom_boxplot()

Show the code

x=data|>
  filter(phase=="pre") |>
  select(score) |>
  unlist()

y=data|>
  filter(phase=="post") |>
  select(score) |>
  unlist()

t.test(x,y,paired=T)


    Paired t-test

data:  x and y
t = -5.7712, df = 79, p-value = 1.476e-07
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -1.8492297 -0.9007703
sample estimates:
mean difference 
         -1.375

Show the code

wilcox.test(x,y,paired=T)


    Wilcoxon signed rank test with continuity correction

data:  x and y
V = 441, p-value = 1.335e-06
alternative hypothesis: true location shift is not equal to 0

Show the code

library(EnvStats)


Attaching package: 'EnvStats'

The following object is masked from 'package:car':

    qqPlot

The following objects are masked from 'package:stats':

    predict, predict.lm

Show the code

oneSamplePermutationTest(x-y)


Results of Hypothesis Test
--------------------------

Null Hypothesis:                 Mean (Median) = 0

Alternative Hypothesis:          True Mean (Median) is not equal to 0

Test Name:                       One-Sample Permutation Test
                                 (Based on Sampling
                                 Permutation Distribution
                                 5000 Times)

Estimated Parameter(s):          Mean = -1.375

Data:                            x - y

Sample Size:                     80

Test Statistic:                  |Sum(x)| = 110

P-value:                         0

ANOVA à un facteur sur Groupes indépendants

Hettmansperger and McKean (2011) étudient l’effet de quatre drogues sur la réduction de cholesterol (LDL) chez les cailles. Attention les groupes sont déséquilibrées !

Analyses descriptives

Faire les analyses descriptives

Show the code

data("quail")
quail |>
  group_by(treat) |>
  summarise(M=mean(ldl),SD=sd(ldl),
            Med=median(ldl),mad=mad(ldl))

# A tibble: 4 × 5
  treat     M    SD   Med   mad
  <fct> <dbl> <dbl> <dbl> <dbl>
1 1      74.5  25.0  68.5 18.5 
2 2      52.3  33.1  35    7.41
3 3      73.8  37.7  62    8.15
4 4      67.6  23.1  65   23.7

Show the code

ggplot(quail,aes(x=treat,y=ldl))+
  geom_boxplot()

Première modélisation : ANOVA paramétrique

mod=aov(ldl~treat,data=quail)
library(car)
Anova(mod,type="III")

Anova Table (Type III tests)

Response: ldl
            Sum Sq Df F value    Pr(>F)    
(Intercept)  55503  1 59.7050 4.519e-09 ***
treat         3189  3  1.1433    0.3451    
Residuals    32536 35                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pour ceux qui préfèrent les contrastes : on peut se dire que le deuxième traitement paraît réduire de façon plus efficace le ldl que les autres. On pourrait tester \(H_0:\mu_2=\frac1{3}(\mu_1+\mu_3+\mu_4)\) versus \(H_1:\mu_2\not =\frac1{3}(\mu_1+\mu_3+\mu_4)\)

# Estimation du modèle linéaire :
mod=lm(ldl~treat,data=quail)
H=as.matrix(c(-1/3,1,-1/3,-1/3))
# Estimation de la différence des moyennes 
mod$coefficients%*%H

          [,1]
[1,] -44.48519

# Estimation de l'erreur standard
sqrt(diag(t(H)%*%vcov(mod)%*%H))

[1] 12.49332

# Calcul de la statistique de stats
z=mod$coefficients%*%H/sqrt(diag(t(H)%*%vcov(mod)%*%H))
# Calcul de la p-value
2*(1-pnorm(abs(z)))

             [,1]
[1,] 0.0003698421

Deuxième modélisation : test de KW

kruskal_test(ldl~treat,data=quail)


    Asymptotic Kruskal-Wallis Test

data:  ldl by treat (1, 2, 3, 4)
chi-squared = 7.1879, df = 3, p-value = 0.06614

Trosième modélisation : permutation

library(lmPerm)
mod=aovp(ldl~treat,data=quail)

[1] "Settings:  unique SS "

Anova(mod,type="III")

Anova Table (Type III tests)

Response: ldl
            Sum Sq Df  F value    Pr(>F)    
(Intercept) 174910  1 188.1536 1.209e-15 ***
treat         3189  3   1.1433    0.3451    
Residuals    32536 35                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Quatrième modélisation ANOVA sur les rangs

library(npsm)
oneway.rfit(y = quail$ldl, g = quail$treat)

Call:
oneway.rfit(y = quail$ldl, g = quail$treat)

Overall Test of All Locations Equal 

Drop in Dispersion Test
F-Statistic     p-value 
   3.916944    0.016394 


    Pairwise comparisons using Rfit 

data:  quail$ldl and quail$treat 

  1      2      3     
2 0.0046 -      -     
3 0.6315 0.0157 -     
4 0.5599 0.0243 0.9069

P value adjustment method: none

Retour sur les contrastes :

H=as.matrix(c(-1/3,1,-1/3,-1/3))
# Estimation de la différence
mod=rfit(ldl~treat-1,data=quail)
mod$coefficients%*%H

     [,1]
[1,]  -22

# Estimation de l'erreur standard
sqrt(diag(t(H)%*%vcov(mod)%*%H))

[1] 6.781191

# Calcul de la statistique de stats
z=mod$coefficients%*%H/sqrt(diag(t(H)%*%vcov(mod)%*%H))
# Calcul de la p-value
2*(1-pnorm(abs(z)))

            [,1]
[1,] 0.001177529

RM Anova à un facteur

data("OBrienKaiserLong")
OBrienKaiserLong |>
  group_by(phase) |>
  summarise(M=mean(score),SD=sd(score),Med=median(score),
            MAD=mad(score))

# A tibble: 3 × 5
  phase     M    SD   Med   MAD
  <fct> <dbl> <dbl> <dbl> <dbl>
1 pre    4.38  1.87     4  1.48
2 post   5.75  2.33     6  2.97
3 fup    6.38  2.20     7  1.48

# Pour ré-ordonner les facteurs de la variable phase

OBrienKaiserLong=OBrienKaiserLong |>
  mutate(phase=factor(phase,levels=c("pre","fup","post")))

ggplot(OBrienKaiserLong,aes(x=phase,y=score))+
  geom_boxplot()

library(lmerTest)

Loading required package: lme4

Loading required package: Matrix


Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

mod=lmer(score~phase+(1|id),data=OBrienKaiserLong)
anova(mod,type=3)

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
phase  167.5   83.75     2   222   38.48 4.487e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

postHoc=emmeans(mod,specs = ~phase)
pwpm(postHoc)

        pre    fup   post
pre  [4.38] <.0001 <.0001
fup  -2.000 [6.38] 0.0215
post -1.375  0.625 [5.75]

Row and column labels: phase
Upper triangle: P values   adjust = "tukey"
Diagonal: [Estimates] (emmean) 
Lower triangle: Comparisons (estimate)   earlier vs. later

On peut aussi faire un contraste par exemple pour tester \(H_0:\mu_{pre}=\frac 1{2}(\mu_{fup}+\mu_{post}).\)

mod1=summary(lmer(score~phase-1+(1|id),data=OBrienKaiserLong))
## Pour l'estimation des moyennes :
moy=mod1$coefficients[,"Estimate"]
## Pour l'estimation de la matrice de covariance
mod1$vcov

3 x 3 Matrix of class "dpoMatrix"
           phasepre  phasefup phasepost
phasepre  0.1857295 0.1585242 0.1585242
phasefup  0.1585242 0.1857295 0.1585242
phasepost 0.1585242 0.1585242 0.1857295

A vous de jouer.

Show the code

H=as.matrix(c(-1,1/2,1/2))

# Estimation de la différence de moyenne :
M=moy%*%H
print(M)
# Estimation de l'erreur standard :
ES=sqrt(diag(t(H)%*%mod1$vcov%*%H))
print(ES)
# Estimation de la statistique de test :
z=M/ES
print(z)
# Estimation de la p-value
2*(1-pnorm(abs(z)))

On peut aussi faire le test de Friedman :

friedman_test(score~phase,data=OBrienKaiserLong)


    Asymptotic Friedman Test

data:  score by
     phase (pre, fup, post) 
     stratified by block
chi-squared = 36.312, df = 2, p-value = 1.303e-08

ou bien de la permutation :

library(permutes)
permutes::perm.lmer(score~phase,data=OBrienKaiserLong,
                    type="anova")

Loading required namespace: buildmer

       Factor df       LRT         F  p
1 (Intercept)  1 396.96809 334.01404 NA
2       phase  2  86.46743  18.26852  0

ANOVA à deux facteurs indépendants

Charger la base serumLH.

On commence par une analyse descriptive des données :

Show the code

data("serumLH")
serumLH=serumLH |>
  mutate(LRF.dose=factor(LRF.dose,levels=c("0","10","50","250","1250")))

serumLH |>
  group_by(light.regime,LRF.dose) |>
  summarise(median(serum))

`summarise()` has grouped output by 'light.regime'. You can override using the
`.groups` argument.

# A tibble: 10 × 3
# Groups:   light.regime [2]
   light.regime LRF.dose `median(serum)`
   <fct>        <fct>              <dbl>
 1 Constant     0                    67 
 2 Constant     10                   80 
 3 Constant     50                  132.
 4 Constant     250                 184.
 5 Constant     1250                202 
 6 Intermittent 0                   101 
 7 Intermittent 10                  167 
 8 Intermittent 50                  284.
 9 Intermittent 250                 390 
10 Intermittent 1250                407

Show the code

ggplot(serumLH,aes(x=LRF.dose,y=serum,color=light.regime))+geom_boxplot()

Anova paramétrique

mod<-lm(serum~LRF.dose*light.regime,data=serumLH)
Anova(mod,type="III")

Anova Table (Type III tests)

Response: serum
                      Sum Sq Df F value   Pr(>F)   
(Intercept)            21600  1  3.5874 0.064012 . 
LRF.dose              130294  4  5.4099 0.001065 **
light.regime            7600  1  1.2623 0.266585   
LRF.dose:light.regime  55099  4  2.2877 0.072866 . 
Residuals             301055 50                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova basée sur les rangs

mod<-raov(serum~LRF.dose*light.regime,data=serumLH)
mod


Robust ANOVA Table
                      DF        RD   Mean RD        F p-value
LRF.dose               4 3027.6735  756.9184 26.86162 0.00000
light.regime           1 1642.3333 1642.3333 58.28334 0.00000
LRF.dose:light.regime  4  451.4553  112.8638  4.00533 0.00678

On veut faire des comparaisons de l’effet dose selon le type de régime (à vous !)

mod=rfit(serum~LRF.dose:light.regime-1,data=serumLH)

Show the code

H=rbind(diag(c(1,1,1,1,1)),diag(c(-1,-1,-1,-1,-1)))
D=mod$coefficients%*%H
ES=sqrt(diag(t(H)%*%vcov(mod)%*%H))
z=D/ES
2*(1-pnorm(abs(z)))

          [,1]       [,2]         [,3]        [,4]         [,5]
[1,] 0.1486039 0.01950028 2.918003e-05 6.16267e-10 7.369958e-10

RM ANOVA à 1 facteur Within et 1 facteur between

On utilise la base de données ObrienKaiserLong. On ajoute le facteur treatment à l’analyse sur ce jeu de données.

Faire les analyses descriptives :

Show the code

OBrienKaiserLong |>
  group_by(phase,treatment) |>
  summarise(M=mean(score),SD=sd(score),Med=median(score),
            MAD=mad(score))

`summarise()` has grouped output by 'phase'. You can override using the
`.groups` argument.

# A tibble: 9 × 6
# Groups:   phase [3]
  phase treatment     M    SD   Med   MAD
  <fct> <fct>     <dbl> <dbl> <dbl> <dbl>
1 pre   control    4.2   1.68     4  1.48
2 pre   A          5     2.15     5  2.97
3 pre   B          4.14  1.80     4  1.48
4 fup   control    4.4   1.80     4  1.48
5 fup   A          7.25  2.20     7  2.97
6 fup   B          7.29  1.43     7  1.48
7 post  control    4     1.80     3  1.48
8 post  A          6.5   2.59     7  2.97
9 post  B          6.57  1.82     6  1.48

Show the code

ggplot(OBrienKaiserLong,aes(x=phase,y=score,color=treatment))+
  geom_boxplot()

mod=lmer(score~phase*treatment+(1|id),data=OBrienKaiserLong)
anova(mod,type=3)

Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
phase           136.789  68.395     2   218 36.7091 1.817e-14 ***
treatment        10.858   5.429     2    13  2.9139   0.09004 .  
phase:treatment  77.000  19.250     4   218 10.3320 1.113e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc=emmeans(mod,specs = ~phase|treatment)
pwpm(posthoc)


treatment = control
       pre    fup   post
pre  [4.2] 0.8626 0.8626
fup   -0.2  [4.4] 0.5549
post   0.2    0.4  [4.0]

treatment = A
        pre    fup   post
pre  [5.00] <.0001 0.0018
fup   -2.25 [7.25] 0.1936
post  -1.50   0.75 [6.50]

treatment = B
        pre    fup   post
pre  [4.14] <.0001 <.0001
fup  -3.143 [7.29] 0.0753
post -2.429  0.714 [6.57]

Row and column labels: phase
Upper triangle: P values   adjust = "tukey"
Diagonal: [Estimates] (emmean) 
Lower triangle: Comparisons (estimate)   earlier vs. later

On peut le faire aussi sous forme de contrastes comme précédemment.

Reprendre l’exemple précédent

On ajoute le facteur group à l’exemple Baumann. Faire les analyses nécessaires.