Supports de formation en biostatistique avec exemples d'application

Supports de formation en biostatistique avec exemples d'application

…

Après avoir vérifié que le modèle que l’on a construit est bien ajusté aux données (voir fiche 41), l’effet des variables explicatives peut être testé. Il existe trois façons de réaliser ce test, que l’on appelle types I, II et III. S’il n’y a qu’une variable explicative, les trois types de test donnent exactement le même résultat. S’il y a plusieurs facteurs parmi les variables explicatives mais que les effectifs sont équilibrés (i.e. s’il y a autant d’individus dans toutes les classes d’un facteur et dans toutes les combinaisons des facteurs en interaction), il en est de même. Les trois types de test donnent des résultats différents s’il y a plusieurs facteurs et que les effectifs ne sont pas équilibrés. Il est important de comprendre d’où viennent ces différences, car elles reflètent des hypothèse statistiques (et donc biologiques) différentes. On prendra comme exemple un modèle du type reponse~A+B+A:B, où A et B sont des facteurs (voir fiche 40 pour une explication de la formule).

Type I : l’analyse est séquentielle, i.e. les termes sont testés dans l’ordre où ils apparaissent dans la formule (NB : lorsqu’un modèle est créé, R réorganise automatiquement la formule pour placer les facteurs seuls d’abord, puis les interactions d’ordre 2, les interactions d’ordre 3. . .). Dans notre exemple on a donc (i) un test de l’effet de A ; (ii) un test de l’effet de B après avoir retiré l’effet de A ; un test de l’effet de A:B après retiré les effets de A et B. Cela implique que les résultats du test pour les termes A et B sont dépendants de l’ordre dans lequel ils apparaissent dans la formule. Les hypothèses sous-jacentes sont rarement celles que l’on veut tester en biologie, aussi le type I est à oublier (excepté bien sûr si c’est exactement ce que l’on souhaite tester). Les tests de type I sont réalisés par la fonction anova(modele). Attention, selon les modèles ce n’est pas le même test qui est appliqué.

Type II : l’analyse est non séquentielle, i.e. indépendante de l’ordre des termes dans la formule. Elle respecte par contre le principe de marginalité, qui stipule que « lorsqu’une interaction a un effet significatif, l’effet des variables de l’interaction, lorsqu’elles sont impliquées dans des termes d’ordre inférieur, est marginal devant celui de l’interaction ». Dit autrement, ce principe dit que si une interaction entre deux facteurs a un effet significatif, l’information est dans cette interaction et il est inutile de regarder l’effet des facteurs pris seuls (étendu à trois facteurs : si l’interaction entre les trois a un effet significatif, il est inutile de s’intéresser aux facteurs pris seuls et aux interactions d’ordre 2 impliquant ces facteurs). Concrètement, on considère donc que lorsque l’on teste l’effet d’un facteur, tous les termes d’ordre supérieur impliquant ce facteur n’ont pas d’effet. Dans notre exemple, cela donne (i) un test de l’effet de A après avoir retiré l’effet de B ; (ii) un test de l’effet de B après avoir retiré l’effet de A ; (iii) un test de l’effet de A:B après avoir retiré les effets de A et B. Cette approche est celle qui est systématiquement utilisée dans cet aide-mémoire, car dans la très grande majorité des cas c’est celle qui est la plus pertinente vis-à-vis des hypothèses biologiques. Les tests de type II sont réalisés par la fonction Anova(modele,type="II")1 ou plus simplement Anova(modele)1 , le type II étant celui utilisé par défaut. Attention, selon les modèles ce n’est pas le même test qui est appliqué.

Type III : l’analyse est non séquentielle et ne respecte pas le principe de marginalité. Le test de l’effet d’un facteur est donc réalisé en prenant en compte les interactions impliquant ce facteur. Dans notre exemple, cela donne (i) un test de l’effet de A après avoir retiré les effets de B et A:B ; (ii) un test de l’effet de B après avoir retiré les effets de A et A:B ; (iii) un test de l’effet de A:B après avoir retiré les effets de A et B. Cette approche, même si elle peut paraître intéressante au premier abord, est très délicate à utiliser car elle revient à considérer qu’un facteur puisse ne pas avoir d’effet seul tout en modifiant l’effet d’un autre facteur (i.e. tout en ayant un effet ailleurs). Excepté si c’est réellement l’hypothèse que l’on veut tester, le type III est donc à oublier car il est rarement pertinent vis-à-vis des hypothèses biologiques. Les tests de type III sont réalisés par la fonction Anova(modele,type="III")1 . Attention, selon les modèles ce n’est pas le même test qui est appliqué.

Lorsqu’un test sur un modèle a révélé un effet significatif d’un facteur (ou d’une interaction impliquant au moins un facteur), il est nécessaire de réaliser des comparaisons multiples pour voir précisément quelles sont les modalités qui diffèrent. La méthode présentée ici est l’une des plus intéressantes car elle prend en compte l’effet des autres variables explicatives du modèle dans les comparaisons. En d’autres termes les comparaisons se font après avoir retiré la variation due aux autres variables explicatives. Ce ne sont donc plus les moyennes brutes (telles qu’on peut les calculer dans la fiche 35) qui sont comparées, mais des moyennes ajustées en fonction des autres variables explicatives. La méthode en question est celle des moyennes des moindres carrés, ou Least Squares Means (abrégé le plus souvent en LSMeans). Remarque 1 : dans le cas d’une interaction entre un facteur et une covariable, ce ne sont pas des moyennes qui sont comparées mais des pentes (les pentes de la relation entre la variable à expliquer et la covariable, calculées pour chaque modalités du facteur). Remarque 2 : cette fiche concerne tous les modèles présentés dans ce document exceptés ceux analysant une réponse nominale à plus de deux catégories (voir fiche 51) ou un décompte d’individus dans plus de deux catégories (voir fiche 63).

Calcul et enregistrement des moyennes/pentes ajustées

Facteur seul ou interaction entre deux facteurs : moyennes La première étape consiste à calculer les moyennes ajustées, qui sont stockées dans un objet appelé LSM. La syntaxe est du type : LSM<-lsmeans(modele,~facteur)1 où facteur est le nom du facteur ou de l’interaction d’intérêt (attention à ne pas oublier le symbole ~). Pour calculer les moyennes séparément pour chaque niveau d’un autre facteur : LSM<-lsmeans(modele,~facteur1|facteur2)1 où facteur1 est le facteur (ou l’interaction) dont on veut comparer les modalités, et facteur2 le facteur pour lequel on veut réaliser les comparaisons à l’intérieur de chaque modalité. Remarque 1 : pour les GLMMs basés sur une loi binomiale négative (créés avec la fonction glmer.nb() ; voir fiche 55), il est nécessaire de préciser en plus via l’argument data le tableau de données sur lequel est basé le modèle. Dans tous les autres cas ce n’est pas nécessaire. Remarque 2 : pour les MLMs (voir fiche 108), les moyennes ajustées peuvent être calculées séparément pour chaque variable à expliquer. Pour ce faire : LSM<-lsmeans(modele,~facteur|rep.meas)1 . Il est impératif d’utiliser la syntaxe rep.meas.

…

Réalisation des comparaisons

 Pour réaliser toutes les comparaisons deux-à-deux : contrast(à.comparer,"pairwise")1 où à.comparer est soit LSM soit pentes. Dans le cas de comparaisons personnalisées, la procédure se fait en deux étapes : > cont.lsmc <- user.cont(contrastes)2 > contrast(à.comparer,"cont")1 où contrastes est la matrice des contrastes qui spécifie les comparaisons à réaliser. La fonction contrast()1 a l’avantage de donner la p-value associée à chaque comparaison, mais oblige à regrouper les modalités (dans des groupes du type a, ab, b.. .) à la main. Une méthode de regroupement automatique existe cependant : cld(à.comparer) 1 . Cette méthode ne renvoyant pas la p-value de chaque comparaison, elle doit être vue comme complémentaire de la première. Remarque 1 : pour pouvoir utiliser la fonction cld()1 il est nécessaire d’avoir installé le package multcompView. Remarque 2 : dans le cas d’un MLM (voir fiche 108), les comparaisons multiples réalisées prennent en compte l’ensemble des variables à expliquer simultanément ; sauf si l’objet LSM a été créé en utilisant rep.meas, auquel cas les comparaisons sont réalisées séparément pour chaque variable à expliquer.

Tests de significativité des moyennes/pentes

Il est parfois intéressant de tester si les moyennes/pentes sont individuellement différentes de 0. La série de tests est réalisée automatiquement via summary(à.comparer,infer=TRUE).

Récupération des moyennes/pentes ajustées

Facteur seul ou interaction entre deux facteurs : moyennes Représenter sur un diagramme en barres les moyennes ajustées par niveau d’un facteur (ou combinaison de deux facteurs) est souvent plus pertinent que de représenter les moyennes brutes. Cela permet d’illustrer la variation réellement due à ce facteur, après avoir retiré la variation due aux autres variables explicatives du modèle.

Pour tous les modèles sauf ceux analysant un rang dans un classement (voir fiche 52 pour la démarche à utiliser dans ce cas), récupérer les moyennes ajustées se fait de la même façon :

— Si la variable à expliquer n’a pas été transformée avant de créer le modèle : valeurs<-as.data.frame(summary(LSM,type="response")). L’objet valeurs est un tableau dont la (les) première(s) colonne(s) correspond(ent) aux modalités du (des) facteur(s) précisé(s) dans lsmeans()1 . La colonne appelée lsmean, response, prob, rate ou hazard contient les moyennes ajustées. La colonne SE contient l’erreur standard associée à chaque moyenne.

— Si la variable à expliquer a été transformée avant de créer le modèle (ce qui en pratique ne concerne que certains LM(M)s) : valeurs<-back.lsmeans(LSM,transform="transfo",add=valeur) 2 où transfo est le nom de la transformation entre guillemets ("log", "sqrt", "inverse" ou "logit") et valeur une éventuelle constante ajoutée à la variable à expliquer avant transformation. L’objet valeurs est un tableau dont la première colonne correspond aux modalités du (des) facteur(s) précisé(s) dans lsmeans()1 . La colonne Mean donne les moyennes ajustées et les colonnes SE.inf et SE.sup donnent les bornes formées par l’erreur standard, qui est nécessairement asymétrique avec une transformation de la variable à expliquer.

…

Une fois les moyennes et erreurs standards récupérées, elles peuvent être représentées sur un diagramme en barres (voir fiche 34). Remarque : dans le cas d’un MLM (voir fiche 108), les moyennes ajustées n’ont de sens que si elles sont calculées pour chaque variable à expliquer. Il est donc impératif d’utiliser rep.meas pour les récupérer. Interaction entre un facteur et une covariable : pentes Les pentes sont récupérées via : valeurs<-as.data.frame(summary(pentes)). L’objet valeurs est un tableau dont la (les) première(s) colonne(s) correspond(ent) aux modalités du (des) facteur(s) précisé(s) dans lstrends()1 . La colonne dont le nom se termine par .trend contient les pentes ajustées. La colonne SE contient l’erreur standard associée à chaque pente. Remarque : en pratique les pentes ne sont directement interprétables que pour les LM(M)s.

44. Sélection de modèle

On utilise le plus souvent la démarche de sélection de modèle lorsque le modèle servira ensuite à faire de la prédiction. L’objectif est d’aboutir au meilleur compromis entre un bon ajustement aux données (qui croît avec le nombre de variables explicatives) et un petit nombre de variables explicatives (car plus on prend en compte de variables explicatives, plus on perd en capacité à généraliser les résultats ce qui est problématique pour de la prédiction). La sélection est basée sur une valeur (appelée critère) représentant ce compromis. Plus elle est faible, meilleur est le compromis. Il existe plusieurs critères, le plus utilisé est le Critère d’Information d’Akaike (AIC).

En pratique, la sélection est réalisée automatiquement à partir du modèle initial (complet). Pour la réaliser : dredge(modele)1 , où modele est le modèle contenant toutes les variables explicatives. Cette fonction calcule l’AIC de tous les modèles possibles à partir du modèle initial et renvoie un tableau classant tous ces modèles. Parmi les arguments facultatifs de la fonction, deux sont particulièrement intéressants :

— m.max=valeur, où valeur est le nombre maximal de variables explicatives à intégrer dans les modèles à tester

— fixed=variables, où variables est un vecteur contenant le nom des variables explicatives à intégrer dans tous les modèles à tester (entre guillemets).

Attention, si n d f < 40, où n est le nombre d’individus et d f le nombre de paramètres estimés par modèle (renvoyé par le fonction dredge()1 ), il faut utiliser l’AICc (AIC corrigé) et non pas l’AIC. En fait, comme cette situation est courante, la fonction dredge() utilise par défaut l’AICc pour classer les modèles. Dans le cas de modèles avec une loi quasipoisson ou quasibinomial, le critère utilisé est le QAIC, dérivé de l’AIC pour les distributions « quasi ». Il est nécessaire pour ces modèles de renseigner la valeur du paramètre de dispersion (obtenu via summary(modele)) grâce à l’argument chat. Pour utiliser un autre critère que l’AICc pour réaliser la sélection, utiliser l’argument rank (en donnant comme valeur "AIC" ou "QAIC" selon les cas).

…

Lorsque l’on utilise une procédure de sélection automatique, il faut être conscient que le modèle avec l’AIC (ou dérivé) le plus faible n’est pas forcément celui qui biologiquement a le plus de sens. Il ne faut donc pas utiliser cette procédure aveuglément mais toujours réfléchir aux variables et interactions qui sont retenues. Il est ainsi possible d’enlever des termes manuellement si ceux-ci ne sont pas pertinents ou trop complexes à interpréter (comme une interaction d’ordre 3 ou 4).

Séries appariées Test de Wald (paramétrique) Pour réaliser le test : wald.ptheo.test(reponse,blocs,p=proba.theo)1 où reponse est la ré- ponse de chaque individu (sous forme numérique ou d’un facteur, en tout cas binaire) et blocs un facteur (aléatoire) contenant le groupe de chaque individu (dans le même ordre que reponse). Si reponse est codée sous forme 0/1, la probabilité du groupe 1 est testée ; si reponse est un facteur, la probabilité de la 2 nde modalité est testée.

…

Les effectifs théoriques sont obtenus via chisq.bin.exp(reponse~facteur,p=proba.theo)1 où reponse et facteur sont des vecteurs contenant la valeur de chaque individu pour les deux variables (dans le même ordre), et où proba.theo est un vecteur contenant la probabilité théorique dans chaque modalité (dans l’ordre des modalités de facteur). Si reponse est codée sous forme 0/1, la probabilité du groupe 1 est testée ; si reponse est un facteur, la probabilité de la 2nde modalité est testée. Le symbole ~ signifie « expliqué par » ou « en fonction de ». Les effectifs théoriques permettent de vérifier si la règle de Cochran est respectée. Remarque : les probabilités théoriques sont indépendantes entre les modalités de facteur. Elles n’ont pas à donner une somme de 1, puisque les modalités ne sont pas comparées entre elles, mais chacune est testée pour sa propre probabilité théorique. Pour réaliser le test : chisq.theo.bintest(reponse~facteur,p=proba.theo)1 . Une p-value significative indique qu’au moins une probabilité diffère de sa valeur théorique, sans pré- ciser la(les)quelle(s). Il est dans ce cas nécessaire de réaliser des comparaisons deux-à-deux pour identifier la (les) probabilité(s) en question, via prop.bin.multcomp(reponse~facteur,p=proba.theo)1 . Il peut arriver que les comparaisons deux-à-deux n’indiquent aucune différence significative, contrairement au test global. Dans ce cas, la solution la plus prudente est de considérer qu’on ne peut pas savoir quelle modalité est responsable du rejet de l’hypothèse nulle dans le test global.

Séries non appariées

Test exact de Fisher (non paramétrique) Pour réaliser le test : fisher.bintest(reponse~facteur)1 où reponse et facteur sont des vecteurs contenant la valeur de chaque individu pour les deux variables (dans le même ordre). Le symbole ~ signifie « expliqué par » ou « en fonction de ». Si la p-value du test est significative, cela indique qu’au moins deux classes du facteur ont un effet différent sur la variable à expliquer (sans préciser lesquelles). La fonction réalise alors automatiquement toutes les comparaisons deux-à-deux possibles par une série de tests exacts de Fisher. Il peut arriver que les comparaisons deux-à-deux n’indiquent aucune différence significative, contrairement au test global. Dans ce cas, la solution la plus prudente est de considérer qu’on ne peut pas savoir quelles probabilités sont responsables du rejet de l’hypothèse nulle dans le test global.

Test du χ 2 d’homogénéité (non paramétrique) Pour réaliser le test : chisq.bintest(reponse~facteur)1 . Si la p-value est significative, les comparaisons deux-à-deux sont réalisées automatiquement par une série de tests du χ 2 d’homogénéité.

Test G d’homogénéité (non paramétrique) Pour réaliser le test : G.bintest(reponse~facteur)1 . Si la p-value est significative, les comparaisons deux-à-deux sont réalisées automatiquement par une série de tests G d’homogénéité.

Séries appariées

Test Q de Cochran (non paramétrique) Conditions : le plan d’expérience doit être en blocs aléatoires complets sans répétition (voir fiche 13). Pour réaliser le test : cochran.qtest(reponse~fact.fixe|fact.alea)1 où variable, fact.fixe et fact.alea sont des vecteurs contenant la valeur de chaque individu (dans le même ordre) pour la variable à expliquer, le facteur (fixe) dont on veut comparer les modalités et le facteur (aléatoire) servant à définir les séries appariées, respectivement. Si la p-value du test est significative, cela indique qu’au moins deux classes du facteur fixe ont un effet différent sur la variable à expliquer (sans préciser lesquelles). La fonction réalise alors automatiquement toutes les comparaisons deux-à-deux possibles par une série de tests des signes de Wilcoxon.

…

Test(s)

L’effet des variables explicatives est testé via Anova(modele)3 (voir fiche 42 pour une explication détaillée des hypothèses testées). Le test réalisé est un test du rapport des vraisemblances (Likelihood Ratio Test ou LR Test) – en fait un test par terme du modèle (i.e. un par ligne du tableau renvoyé). Si une variable explicative a un effet significatif, cela indique qu’elle influence le rapport entre au moins deux classes de la variable réponse. Autrement dit, la probabilité qu’un individu soit d’une classe donnée par rapport à une autre classe donnée. Pour identifier quel(s) rapport(s) est (sont) significativement influencé(s) par la variable explicative en question : test.multinom(modele,variable)2 où variable est la variable explicative dont on souhaite étudier l’influence en détail. Ce qui est renvoyé par test.multinom()2 dépend de la nature de la variable explicative : — Pour une variable explicative quantitative : la fonction renvoie un tableau où chaque ligne correspond à un rapport entre deux classes de la variable réponse (la syntaxe A|B signifie « probabilité de A par rapport à la probabilité de B »).

Le signe du coefficient indique le sens de la relation entre un rapport et la variable explicative. L’odds ratio est une façon simple de comprendre ce rapport, qui indique de combien il varie pour une augmentation d’une unité de la variable explicative. Remarque : il n’est pas surprenant que les p-values des tests par rapport soient moins nettes que la p-value globale obtenue via Anova(modele)3 . D’abord parce qu’un test global ne se réduit pas à une somme de tests individuels, mais aussi car ce ne sont pas les mêmes tests qui sont réalisés : test du rapport des vraisemblances au niveau global, tests de Wald pour chaque rapport. Le second est moins puissant que le premier (voir fiche 17).