Déterminer la taille des échantillons pour un futur test de Student se fait avec le d de Cohen, une mesure de la taille de l'effet.

La magnitude du d de Cohen s'analyse de la façon suivante : 

• Petit effet : ≈0.2

• Effet moyen : ≈0.5

• Grand effet : ≈0.8

En réalité, cette taille de l'effet peut varier en-dessous de 0 et au-delà de 1 et ne donne pas le même résultat selon que l'on le calcule en faisant l'écart entre la moyen X1 et X2 ou X2 et X1. On va donc prendre une valeur absolue.

Ce faisant, si je compare la croissance d'épis de Maïs ou les dépenses d'établissements, j'aurais du mal à savoir à l'avance quelle différence correspond à quelle taille d'effet. Si mes dépenses ont augmenté de 20%, c'est un grand effet ? un moyen ?

C'est pourquoi, il faut tenter d'estimer la taille de l'effet avant de se lancer dans le calculs de la taille de l'échantillon.

Si on ne veut pas estimer la taille de l'effet, on peut simuler ses données un grand nombre de fois (solution simple), ou sinon, utiliser cette fonction cohens_d() :

cohens_d <- function(mean1, mean2, sd1, sd2, n1=1000, n2=1000) {

  # Vérification des entrées pour s'assurer qu'elles sont valides

  if (n1 <= 1 || n2 <= 1) {

    stop("Les tailles d'échantillon (n1 et n2) doivent être supérieures à 1.")

  }

  if (sd1 <= 0 || sd2 <= 0) {

    stop("Les écarts-types (sd1 et sd2) doivent être positifs.")

  }

  # Calcul de l'écart-type poolé

  pooled_sd <- sqrt(((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / (n1 + n2 - 2))

  # Calcul de Cohen's d

  d <- abs((mean1 - mean2) / pooled_sd)

  return(d)

}

Allez, prenons un exemple. D'après la littérature scientifique, le QI moyen des femmes est de 100, celui des hommes de 2 à 5 points supérieur. En revanche, la variabilité de l'écart-type de 15-15,3 pour les humains impliquant que 95% des QI sont entre 70 et 130, diffère pour les hommes (écart-type 20% plus élevé) de sorte que les QI bas (en-dessous de 70) sont à 66% pour les hommes.

Simulons donc des QI en appliquant une correction pour éviter que les QI partent trop bas.

Bien, maintenant, combien faudrait-il d'hommes et de femmes pour démontrer une différence de QI entre les deux sexes dans 95% des cas ?

1. D'abord, estimons l'effet attendu

cohens_d(100,(100+(2+5)/2), 13.3,17)

Considérons à ce stade que l'on n'a aucune idée du nombre de filles (n1) qu'il faudrait ni de garçons (n2) sachant que tant que n1 et n2 sont équivalents, n n'a pas d'effet sur la détermination sur la taille de l'effet.

On trouve ici un effet très très faible, négatif même, 0.23, ce qui correspond bien à la taille d'effet attendu, petit, autour de 0.2.

2. Ensuite, estimons le nombre d'individus qu'il faudrait pour démontrer, avec une puissance de 95% pour un alpha de 0.05, que les hommes ont un QI moyen supérieur à celui des femmes.

# Installer et activer le package {pwr}

install.packages("pwr") ; library(pwr)

effet <- 0.23 # Taille de l'effet (Cohen's d) - ici très faible. Doit être supérieur à 0. Ne pas mettre d'effet négatif ou trop proche de 0 : il peut être ignoré par la fonction sans message d'erreur.

alpha <- 0.05 # Puissance souhaitée

puissance <- 0.95 # Analyse de puissance pour un test t de Student, ici on cherche 95%

result <- pwr.t.test(d = effet,

     sig.level = alpha,

     power = puissance,

     type = "two.sample") 

# Afficher le résultat

print(result)  

Il nous propose n = 492, ce qui correspond au nombre d'individus recommandé par groupe.

3. Et si on vérifiait ?

pval <- c() ; pop <- 492*2

for (i in 1:10000) {

   IQ <- rnorm(pop,c(100,100+(2+5)/2), c(13.3,17))

   sexe <- rep(c("F","M"),pop/2)

   # Corriger les QI < 50 pour éviter qu'ils descendent trop bas

   IQ[IQ<50] <- IQ[IQ<50] -(IQ[IQ<50]-50)+rnorm(length(IQ[IQ<50]), 0,5)

   # Arrondi des IQ

   IQ <- round(IQ, 0)

   pval <- c(pval,t.test(IQ~sexe,alternative='less')$p.value)

}

quantile(pval, probs=.95)

On notera l'option alternative='less' qui montre que l'on cherche à retenir les p-value qui montrent que les hommes ont une moyenne plus élevée ou plus exactement, que les femmes ont une moyenne plus faible (less). On ne sait pas toujours s'il faut mettre less ou greater pour que la comparaison se femmes entre la catégorie 1 et 2 : un histogramme sur les p-values permet de le confirmer ou le bilan d'un simple test de Student.

Ce résultat montre que le seuil qui sépare 95% des p-values les plus basses des plus hautes est de 0.02.

En réalité, l'on a une puissance de 97,4% :

prop.table(table(ifelse(pval<0.05,"OK","nonOK")))

L'analyse de puissance a être légèrement en décalage par rapport à la la simulation, même si on élimine la correction des QI les plus bas (écrasement de la queue). Une simulation sur bootstrap est-elle plus fiable ? Les résultats restent proches. La simulation montre que l'on pourrait atteindre la puissance souhaité avec un nombre d'individus par groupe de seulement n = 430.

En revanche, dans un situation classique, où l'on fait une estimation du d de Cohen "au doigt mouillé" - j'attends un petit effet, donc d = 0.2 -, on serait arrivé à une suggestion de n de 650. Dans ce cas, une étude de puissance sans estimation précise ou sans simulation, va nous inciter à étudier 50% d'hommes et de femmes en plus, ce qui peut limiter notre capacité à faire des études.

Lorsque l'on veut estimer la taille des groupes pour une future ANOVA, on suit la même démarche qu'en haut : paramétrer la taille de l'effet, le nombre de groupes, l'alpha et la puissance.

Il faut, une fois de plus, avoir idéalement une idée précise de la taille de l'effet, ou revenir sur cette logique de 0.2 (effet petit), 0.5 (moyen), 0.8 (grand). Cette taille d'effet porte ici le nom de f de Cohen.

Le f de Cohen peut s'estimer à partir d'un échantillon grâce à la fonction suivante :

cohen_f <- function(quanti, quali) {

  # Vérifier que quali est un facteur

  if (!is.factor(quali)) {

    stop("La variable qualitative doit être un facteur (utilisez as.factor()).")

  }

  # Ajuster le modèle ANOVA

  anova_model <- aov(quanti ~ quali)

  # Extraire les sommes des carrés

  anova_table <- summary(anova_model)[[1]]

  ss_effect <- anova_table["quali", "Sum Sq"]  # Somme des carrés de l'effet

  ss_error <- anova_table["Residuals", "Sum Sq"]  # Somme des carrés résiduelle

  # Calcul du f de Cohen

  f <- sqrt(ss_effect / ss_error)  

  # Autre façon via eta^2

  #ss_total <- sum(anova_table[, "Sum Sq"])  # Somme totale des carrés

  #eta2 <- ss_effect / ss_total

  # Calcul de f à partir de eta^2

  #f <- sqrt(eta2 / (1 - eta2))

  return(f)

}

Si on l'applique à un jeu de données iris, on peut ainsi développer l'exemple suivant : imaginons que j'ai échantillonné 5 iris de 3 espèces. J'ai leur variabilité, je vais l'utiliser pour estimer le f de Cohen et ainsi le nombre d'iris qu'il me faudrait pour voir dans 95% les différences comme significatives sur la largeur de sépales si elles s'avèrent réellement significatives.

data(iris)

Ech_indice <- c(sample(1:50,5),sample(51:100,5),sample(101:150,5))

Ech <- iris[Ech_indice,]

cohen_f(Ech$Sepal.Width,Ech$Species)

L'effet ici semble plutôt fort, 0.8, allons y pour notre étude de puissance.

effet <- 0.8 # Taille de l'effet (f)

gp <- 3 # Nombre de groupes

alpha <- 0.05 # Limite pour considéré une p-value comme significative

puissance <- 0.95 # Puissance souhaitée

# Analyse de puissance pour une ANOVA 

result <- pwr.anova.test(f = effet,

k = gp, sig.level = alpha,

power = puissance) # Afficher le résultat print(result) 

# Afficher le résultat

print(result)

Le résultat dit que 9 iris de chaque espèce devrait suffir pour atteindre une puissance de 0.95%.

Pour une analyse de contingence, la taille de l'effet semesure  à travers le coefficient de contingence w de Cohen.

calculate_w <- function(tab) {

  if (!is.table(tab)) stop("L'entrée doit être une table de contingence (table(x, y)).")

  # Total des observations

  N <- sum(tab)

  # Valeurs observées (O) et attendues (E)

  observed <- as.vector(tab)

  expected <- as.vector(outer(rowSums(tab), colSums(tab)) / N)

  # Calcul de la statistique

  w <- sqrt(sum((observed - expected)^2 / expected) / N)

  return(w)

}

Imaginons que j'ai une maladie qui touche plus les femmes que les hommes. Je vais estimer ma taille d'effet sur une quarantaine des personnes.

sexe <- sample(c("F","M"),40,replace=T)

maladie <- sample(c("M","NM"),40,replace=T,p=c(0.2,0.8))

maladie[sexe=="F"] <- sample(c("M","NM"),length(maladie[sexe=="F"]),

replace=T,p=c(0.4,0.6))

table(sexe,maladie) -> tab

calculate_w(tab)

J'ai ainsi une taille d'effet w = 0.4.

effet <- 0.4 # Effet attendu (w)

alpha <- 0.05 # Valeur seuil de p-value

puissance <- 0.95 # Puissance souhaitée

ddl <- 1 # degré de liberté - n-1 * m-1 (n et m : nb de lignes et cols)

# Analyse de puissance pour un test du chi-carré 

result <- pwr.chisq.test(w = effet,

df = ddl,

sig.level = alpha,

power = puissance) 

# Afficher le résultat 

print(result) 

Ici, l'on a un N suggéré de 81, soit une population totale de 81 individus pour réussir à atteindre une puisse de 95%.

Prenons l'exemple du lien entre la taille du cerveau et le quotient intellectuel. On trouve une corrélation significative entre les deux de r = 0.24. Combien faudrait-il de personnes pour le démontrer ?

library(pwr)

effet <- 0.24 # Taille de l'effet (r, coeff de corrélation)

alpha <- 0.05 # Puissance souhaitée

puissance <- 0.95 # Analyse de puissance pour un test t de Student

result <- pwr.r.test(r = effet,

sig.level = alpha,

power = puissance) 

# Afficher le résultat

print(result) 

Cette étude suggère une qu'il faut n = 218 valeurs pour atteindre une puissance de 95%.

Pour en savoir plus :

• Déterminer la taille de son échantillon pour une expérience : ce lien dirige vers deux fonctions bootstrap qui permettent rapidement de 1) connaître la taille d'échantillons nécessaire pour atteindre la puissance souhaitée, 2) contrôle la taille d'effet que sera capable de révéler une certaine taille d'échantillon ou 3) contrôler la taille d'échantillon qu'il aurait fallu pour voir une corrélation à partir d'un premier échantillon mesurée. 

• Apprendre à simuler ses données pour déterminer ses tailles d'échantillon sans une analyse de puissance.