################################################################################ # Statistiques et tests d'hypothèses, simulations avec le logiciel R # Gaëlle Lelandais # # Prérequis : (obligatoires) définitions de population, échantillon, individu et variable; # (conseillés) différence entre statistique descriptive et inférentielle. # # Merci aux contributeurs de ce script : # -- Les stagiaires du DU "Création, Analyse et Valorisation de Données Omiques" # -- Thomas Denecker ################################################################################ ## # Sauf mention contraire, ce contenu est mis à disposition selon les # termes de la licence Creative Commons Attribution - Partage dans les mêmes # conditions 4.0 International (CC BY-SA 4.0) ## ################################################################################ # Chargement des librairies et des fonctions R ################################################################################ # Chargement des librairies (installées préalablement) library("tidyverse") library("dplyr") library("gridExtra") library("cowplot") # Chargement des fonctions (nécessaires aux simulations) source("Script_Fonctions.R", encoding="utf-8") ################################################################################ # Partie 1 : Présentation de la première population étudiée : "pop0" ################################################################################ pop0 = readRDS(file = "pop0.rds") # --> Cette population est composée d'un ensemble de mesures d'une variable # d'intérêt, obtenues sur les individus (ou unité statistique) qui compose la # population. # La population est composée de : length(pop0) # 10000 valeurs de la variable d'intérêt. # Nous souhaitons étudier la valeur moyenne de cette variable d'intérêt, # dans la population. Ici il est possible de faire le calcul. mean(pop0) # --> Cette valeur est unique (à noter), elle ne changera pas tant que les individus # qui composent la population n'auront pas changé. # --> Habituellement, cette valeur moyenne est inconnue. C'est la raison pour laquelle # des études statistiques sont réalisées ! ################################################################################ # Illustration de la problématique des fluctuations d'échantillonnage ################################################################################ #------------------------------------------------------------------------------- # Création d'un premier échantillon #------------------------------------------------------------------------------- # Tirage aléatoire de 10 individus dans la population ech0 = echantillonnage(pop0, taille = 10) print("Les valeurs de la variable dans l'échantillon sont :") ech0 # La valeur moyenne des mesures dans l'échantillon est-elle proche de ce que l'on # a dans la population ? mean(ech0) # --> Avons nous été chanceux lors de la création de cet échantillon ? # En effet, plus la valeur moyenne de l'échantillon est proche de la valeur # moyenne de la population, plus nous avons été "chanceux" au moment du tirage. #------------------------------------------------------------------------------- # Création de plusieurs échantillons #------------------------------------------------------------------------------- # Nous allons maintenant créer plusieurs échantillons successivement. # Les échantillons ont tous la même taille et un même individu peut être # sélectionné dans plusieurs échantillons (tirage avec remise). # # Les moyennes, calculées à partir de chaque échantillons sont conservées # dans un vecteur. vecteurMoyennes = echantillonnageMultiple(pop0, taille = 10, hist = T, nbEchantillon = 10, valeurRef = mean(pop0)) # ATTENTION : Sur le graphique, seuls 5 échantillons sont représentés à gauche, # mais l'histogramme (à droite) est bien construit à partir des 10 échantillons. # Ainsi, # --> Selon les échantillons créés, la moyenne est plus ou moins proche de # celle de la population. Parfois, le tirage est plus "heureux" que d'autres # fois. # --> C'est la problématique de "Fluctuations d'échantillonnage" ! # --> Représentation histogramme associée, avec en abscisse les # intervalles de valeurs des moyennes et en ordonnée les effectifs. # Et maintenant, que se passe t-il si un plus grand nombre d'échantillons est créé ? # (l'histogramme est tracé lors de l'appel de la fonction cette fois) echantillonnageMultiple(pop0, taille = 10, hist = T, nbEchantillon = 1000, valeurRef = mean(pop0)) # --> L'histogramme obtenu a une forme de "courbe en cloche" --> Loi normale ? # Si nous répètons un très grand nombre de fois la création d'un échantillon, # la moyenne des moyennes obtenues correspond la moyenne de la population ! # Les résultats précédents montrent un élément intéressant : # --> Il y a moins de chance d'obtenir un échantillon avec une moyenne éloignée # de la valeur moyenne de la population, qu'un échantillon avec une moyenne proche # de la valeur moyenne de la population. # --> Problème : Je ne fais jamais cela en biologie (créer 1000 fois le # même échantillon). Création de UN SEUL échantillon... # Comment faire pour "réduire" mon risque de me tromper avec UN SEUL échantillon ? # Une bonne idée est d'augmenter la taille des échantillons, en effet : vecTaille100 = echantillonnageMultiple(population = pop0, taille = 100, hist = T, nbEchantillon = 1000, valeurRef = mean(pop0)) # --> Les valeurs moyennes des échantillons sont bien plus souvent proches de la moyenne # de la population ! Le risque de se tromper est moindre. # --> Avec de grands échantillons, les fluctuations d'échantillonnage sont plus faibles. #------------------------------------------------------------------------------- # Un peu de théorie : Fonction de densité de probabilité de la moyenne d'un échantillon # Voir cours. #------------------------------------------------------------------------------- # Ainsi la loi normale est très utile, pour réaliser une procédure d'inférence # statistique, CAD extrapoler des observations réalisées à partir d'échantillons sur # une population (normalement inconnue !). # Voir cours. # Dernière question : faut-il pour cela travailler avec une variable qui suit # elle même une loi normale ? #------------------------------------------------------------------------------- # Distribution de la population #------------------------------------------------------------------------------- hist(pop0, main ="Variable distribution (pop0)", col = adjustcolor("red", alpha.f = 0.3), xlab = "Value") # -> Et bien non ! Tout cela grâce au théorème central limite ! # Voir cours. #=============================================================================== # Partie 2 : Tests d'hypothèses : étude des populations A, B et C #=============================================================================== # Nous allons aborder la problématique des tests d'hypothèses, dans le cas de la # comparaison de moyennes. Il s'agira par exemple de comparer les valeurs d'une variable # d'intérêt, entre deux conditions : la condition de référence (notée A) et # une condition testés (notée B ou C). # Les mesures possibles de la variable d'intérêt ont été simulées et sont disponibles # dans les objets R nommés popA (référence), popB et popC. # (ATTENTION :: Ces trois population sont différentes de la population pop0 étudiée dans la # première partie). popA = readRDS(file = "popA.rds") popB = readRDS(file = "popB.rds") popC = readRDS(file = "popC.rds") #------------------------------------------------------------------------------- # Visualisation des mesures de la variable dans les deux premières populations #------------------------------------------------------------------------------- # Cas 1 : Comparer des échantillons issus des populations A et B multiHistogrammes(pop1 = popA, pop2 = popB, names = c("popA", "popB")) # --> Que pensez vous des différences entre les mesures de la variable étudiées, # entre ces deux populations ? # Ces deux populations sont DIFFERENTES. Nous n'avons aucun doute à avoir # ici (intérêt de travailler avec des données simulées). En effet, les moyennes de la # variable dans chacune des populations sont : mean(popA) mean(popB) #------------------------------------------------------------------------------- # Problématique des tests d'hypothèses (attention aux embrouilles !) #------------------------------------------------------------------------------- # La problématique des tests d'hypothèses est la suivante. # Nous disposons de mesures de la variables dans des échantillons. Nous souhaitons # savoir si ces échantillons proviennent : # --> D'une même population (Hypothèse H0) # --> De population différentes (Hypothèse H1) # Exemple, vous disposez de ces deux échantillons "mystères". Le premier # est issu de la population de de référence (popA). Mais qu'en est-il du second ? # Est-il issu de la population de référence aussi ? (hypothèse H0) # Oui bien est-il issu de la population test ? (hypothèse H1) echMystere1 = readRDS(file = "echMystere1.rds") echMystere2 = readRDS(file = "echMystere2.rds") # Calculs des moyennes : mean(echMystere1) mean(echMystere2) # --> Les moyennes calculées à partir des deux échantillons sont # différentes. Oui mais ... # Le problème est l'origine de cette différence : Fluctuation d'échantillonnage (H0) # ou bien populations différentes (en plus des fluctuations d'échantillonnage, H1) ? #------------------------------------------------------------------------------- # Application du test de Student (classique en analyse de données, # les conditions d'application sont correctes) #------------------------------------------------------------------------------- # Posons les hypothèses : # H0 : Les différences observées sont dues UNIQUEMENT aux fluctuations d'échantillonnage. # -> Les deux échantillons proviennent de la même population de référence # (les différences observées # entre les échantillons mystères ne sont pas "significatives"). # H1 : Les différentes observées sont dues aux fluctuations d'échantillonnage ET au fait que les # -> échantillons proviennent de deux populations différentes. t.test(echMystere1, echMystere2) # --> Pour réviser la manière d'écrire le paramètres de Student, voir Cours. # Rappel des risques : # Risque de 1ère espèce : Probabilité de rejeter H0 alors que H0 est vraie. # Risque de 2ème espèce : Probabilité de ne pas rejeter H0 alors que H1 est vraie. # La vraie réponse est .... (à donner en classe !) #------------------------------------------------------------------------------- # Mise en application des calculs des risques. #------------------------------------------------------------------------------- # Grâce aux simulations, nous connaissons la bonne réponse (!). Nous allons évaluer # la capacité du test d'hypothèse à nous donner la bonne réponse... echA1 = echantillonnage(population = popA, taille = 10) echA2 = echantillonnage(population = popA, taille = 10) echB = echantillonnage(population = popB, taille = 10) t.test(echA1, echA2) # -> Obtenons nous la bonne réponse ? (= pas de rejet de H0) t.test(echA1, echB) # -> Obtenons nous la bonne réponse ? (= rejet de H0) #------------------------------------------------------------------------------- # Evaluation empirique des risques #------------------------------------------------------------------------------- # La fonction ci dessous permet de répéter 1) la création des échantillons, # 2) l'application du test de Student et 3) la comparaison de la réponse obtenue # avec la réponse attendue. repetStudent(pop1 = popA, pop2 = popA , tailleEch1 = 10, tailleEch2 = 10, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H0", verbose = T) repetStudent(pop1 = popA, pop2 = popB , tailleEch1 = 10, tailleEch2 = 10, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H1", verbose = T) # Cas plus difficile avec la population C multiHistogrammes(popA, popC, c("popA", "popC")) repetStudent(pop1 = popA, pop2 = popC , tailleEch1 = 10, tailleEch2 = 10, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H1", verbose = T) # Finalement, comment éviter les erreurs ? # N'oubliez pas que dans un cas réel, vous ne connaissez pas (en plus) les écarts # entre les moyennes de la variables entre les populations (on se sait pas si le différentiel # existe entre la condition de référence A, et les populations tests B ou C). # Si les échantillons sont plus grandes : repetStudent(pop1 = popA, pop2 = popA , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H0", verbose = T) repetStudent(pop1 = popA, pop2 = popB , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H1", verbose = T) repetStudent(pop1 = popA, pop2 = popC , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 10, pvalue = 0.05, reponseAttendue = "H1", verbose = T) # Si le nombre de tests est plus grand. Quelle évaluation du risque de première espèce ? # Probabilité de rejeter H0 alors que H0 est vraie :) repetStudent(pop1 = popA, pop2 = popA , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 100, pvalue = 0.05, reponseAttendue = "H0", verbose = T) # Si le nombre de tests est plus grand. Quelle évaluation du risque de deuxième espèce ? # Probabilité de ne pas rejeter H0 alors que H1 est vraie :) repetStudent(pop1 = popA, pop2 = popB , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 100, pvalue = 0.05, reponseAttendue = "H1", verbose = T) # -> Ici le cas est favorable. repetStudent(pop1 = popA, pop2 = popC , tailleEch1 = 100, tailleEch2 = 100, nbrTest = 100, pvalue = 0.05, reponseAttendue = "H1", verbose = T) # -> Ici le cas est défavorable. # BIlan : # -- Effet de la taille des échantillons # -- Effet des différences entre les populations comparées # -- Autres remarques ? #------------------------------------------------------------------------------- # Test de Student répété Taille de l'échantillon #------------------------------------------------------------------------------- # A VS A, B et C vecErreurAA = NULL vecErreurAB = NULL vecErreurAC = NULL for( sizeEch in seq(5,100, 5)){ vecErreurAA = c(vecErreurAA,repetStudent(popA,popA , sizeEch, sizeEch, 100, reponseAttendue = "H0", verbose = F)) vecErreurAB = c(vecErreurAB,repetStudent(popA,popB , sizeEch, sizeEch, 100, reponseAttendue = "H1", verbose = F) ) vecErreurAC = c(vecErreurAC,repetStudent(popA, popC, sizeEch, sizeEch, 100, reponseAttendue = "H1", verbose = F) ) } plot(seq(5,100, 5),vecErreurAA, type = "l", col = "red", las = 2, xlab = "Sample size", ylab = "Percentage of errors", main = "Sample size effect - PopA", ylim = c(0, 100)) lines(seq(5,100, 5),vecErreurAB, col = "royalblue") lines(seq(5,100, 5),vecErreurAC, col = "forestgreen") legend("topright", legend = c("A VS A", "A VS B", "A VS C"), lty = 1, col = c("red", "royalblue", "forestgreen"), box.lty = 0, inset = 0.01) # Idéalement, la taille des échantillons devrait s'adapter à l'intensité de l'effet biologique recherché. # --> Intérêt des manipes préliminaires :)