data(LC19_S1, package = "hecedsm")
<- model.matrix( # Matrice du modèle
modmat ~ familiarity + consistency,
data = LC19_S1)
tail(modmat, n = 5L) # Imprimer les premières 5 lignes
#> (Intercept) familiarity consistencyinconsistent
#> 92 1 6 1
#> 93 1 4 1
#> 94 1 7 1
#> 95 1 7 1
#> 96 1 7 1
dim(modmat) # dimension de la matrice du modèle
#> [1] 96 3
4 Régression linéaire
4.1 Introduction
Le modèle de régression linéaire, ou modèle linéaire, est l’un des outils les plus polyvalents pour l’inférence statistique. La régression linéaire est principalement utilisée pour évaluer les effets des variables explicatives (souvent l’effet d’une manipulation ou d’un traitement dans un cadre expérimental) sur la moyenne d’une variable réponse continue, ou pour la prédiction. Un modèle linéaire est un modèle qui décrit la moyenne d’une variable réponse continue
Dénotons par
Pour simplifier la notation, nous regroupons les observations dans un vecteur
En supposant que la variable réponse provient d’une famille de localisation, nous pouvons réécrire le modèle linéaire en termes de la moyenne plus un aléa,
Le modèle linéaire normal ou gaussien spécifie que les réponses suivent une loi normale, avec
4.1.1 Exemples
Considérons quelques exemples de jeux de données qui serviront à illustrer les méthodes par la suite.
Exemple 4.1 (Cohérence de descriptions de produits) L’étude 1 de Lee et Choi (2019) (base de données LC19_S1
, paquet hecedsm
) considère l’impact sur la perception d’un produit de la divergence entre la description textuelle et l’image. Dans leur première expérience, un paquet de six brosses à dents est vendu, mais l’image montre soit un paquet de six, soit une seule). Les auteurs ont également mesuré la familiarité préalable avec la marque de l’article. Les prodeval
, en fonction de la familiarité de la marque familiarity
, un nombre entier allant de 1 à 7, et une variable binaire pour le facteur expérimental consistency
, codé 0
pour des descriptions d’image/texte cohérentes et 1
si elles sont incohérentes. La matrice du modèle qui en résulte est alors de dimension prodeval
est fortement discrétisée.
Exemple 4.2 (Méthodes d’apprentissage de compréhension de lecture) La base de données BSJ92
du paquet hecedsm
contient les résultats d’une expérience de Baumann, Seifert-Kessell, et Jones (1992) sur l’efficacité de différentes stratégies de lecture sur la compréhension d’enfants.
Soixante-six élèves de quatrième année ont été assignés au hasard à l’un des trois groupes expérimentaux suivants : (a) un groupe « Think-Aloud » (TA), dans lequel les élèves ont appris diverses stratégies de contrôle de la compréhension pour la lecture d’histoires (par exemple : auto-questionnement, prédiction, relecture) par le biais de la réflexion à haute voix; (b) un groupe lecture dirigée-activité de réflexion (DRTA), dans lequel les élèves ont appris une stratégie de prédiction-vérification pour lire et répondre aux histoires; ou (c) un groupe activité de lecture dirigée (DRA), un groupe contrôle dans lequel les élèves se sont engagés dans une lecture guidée non interactive d’histoires.
Les variables d’intérêt sont group
, le facteur pour le groupe expérimental, soit DRTA
, TA
et DR
ainsi que les variables numériques pretest1
et posttest1
, qui donnent le score (sur 16) sur le test pré-expérience pour la tâche de détection des erreurs.
Les données sont balancées puisqu’il y a 22 observations dans chacun des trois sous-groupes. Les chercheurs ont appliqué une série de trois évaluations: le test 1 de détection d’erreurs, le test 2 consistant en un questionnaire de suivi de compréhension, et le test 3 standardisé Degrees of Reading Power). Les tests 1 et 2 ont été administrés à la fois avant et après l’intervention: cela nous permet d’établir l’amélioration moyenne de l’élève en ajoutant le résultat du test pré-intervention comme covariable. Les tests 1 étaient sur 16, mais celui administré après l’expérience a été rendu plus difficile pour éviter les cas d’étudiants obtenant des scores presque complets. La corrélation entre le pré-test et le post-test 1 est
Exemple 4.3 (Discrimination salariale dans un collège américain) On s’intéresse à la discrimination salariale dans un collège américain, au sein duquel une étude a été réalisée pour investiguer s’il existait des inégalités salariales entre hommes et femmes. Le jeu de données college
contient les variables suivantes:
salaire
: salaire de professeurs pendant l’année académique 2008–2009 (en milliers de dollars USD).echelon
: échelon académique, soit adjoint (adjoint
), aggrégé (aggrege
) ou titulaire (titulaire
).domaine
: variable catégorielle indiquant le champ d’expertise du professeur, soit appliqué (applique
) ou théorique (theorique
).sexe
: indicateur binaire pour le sexe,homme
oufemme
.service
: nombre d’années de service.annees
: nombre d’années depuis l’obtention du doctorat.
Exemple 4.4 (Suggestion de montants de dons) L’étude 1 de Moon et VanEpps (2023) (données MV23_S1
, paquet hecedsm
) porte sur la proportion de donateurs à un organisme de charité et le montant de leurs dons. Les participants au panel en ligne avaient la possibilité de gagner 25$ et de faire don d’une partie de cette somme à l’organisme de leur choix. Les données fournies incluent uniquement les personnes qui n’ont pas dépassé ce montant et qui ont indiqué avoir fait un don d’un montant non nul.
Exemple 4.5 (Un emballage en carton supplémentaire est-il considéré comme plus écologique ?) Sokolova, Krishna, et Döring (2023) tient compte des préjugés des consommateurs lorsqu’il s’agit d’évaluer le caractère écologique des emballages. Des produits tels que les céréales sont emballés dans des sacs en plastique, eux-mêmes recouverts d’une boîte. Ils supposent (et constatent) que, paradoxalement, les consommateurs ont tendance à considérer l’emballage comme plus écologique lorsque la quantité de carton ou de carton entourant la boîte est plus importante, ce qui n’est pas le cas. Nous examinons dans la suite les données de l’étude 2A, qui mesure la perception du respect de l’environnement (PEF, variable pef
) en fonction de la proportion
d’emballage en carton (soit aucun, soit la moitié de la surface du plastique, soit la même, soit le double).
4.1.2 Analyse exploratoire des données
L’analyse exploratoire des données est une procédure itérative par laquelle nous interrogeons les données, en utilisant des informations auxiliaires, des statistiques descriptives et des graphiques, afin de mieux informer notre modélisation.
Elle est utile pour mieux comprendre les caractéristiques des données (plan d’échantillonnage, valeurs manquantes, valeurs aberrantes), la nature des observations, qu’il s’agisse de variables réponse ou explicatives et les interrelations entre variables.
Voir le Chapitre 11 de Alexander (2023) pour des exemples. En particulier, il convient de vérifier
- que les variables catégorielles sont adéquatement traitées comme des facteurs (
factor
). - que les valeurs manquantes sont adéquatement déclarées comme telles (code d’erreur, 999, etc.)
- s’il ne vaudrait mieux pas retirer certaines variables explicatives avec beaucoup de valeurs manquantes.
- s’il ne vaudrait mieux pas fusionner des modalités de variables catégorielles si le nombre d’observation par modalité est trop faible.
- qu’il n’y a pas de variable explicative dérivée de la variable réponse
- que le sous-ensemble des observations employé pour l’analyse statistique est adéquat.
- qu’il n’y a pas d’anomalies ou de valeurs aberrantes (par ex., 999 pour valeurs manquantes) qui viendraient fausser les résultats.
Exemple 4.6 (Analyse exploratoire des données college
) Une analyse exploratoire des données est de mise avant d’ébaucher un modèle. Si le salaire augmente au fil des ans, on voit que l’hétérogénéité change en fonction de l’échelon et qu’il y a une relation claire entre ce dernier et le nombre d’années de service (les professeurs n’étant éligibles à des promotions qu’après un certain nombre d’années). Les professeurs adjoints qui ne sont pas promus sont généralement mis à la porte, aussi il y a moins d’occasions pour que les salaires varient sur cette échelle.
Ainsi, le salaire augmente avec les années, mais la variabilité croît également. Les professeurs adjoints qui ne sont pas promus sont généralement mis à la porte, aussi il y a moins d’occasions pour que les salaires varient sur cette échelle. Il y a peu de femmes dans l’échantillon: moins d’information signifie moins de puissance pour détecter de petites différences de salaire. Si on fait un tableau de contingence de l’échelon et du sexe, on peut calculer la proportion relative homme/femme dans chaque échelon: 16% des profs adjoints, 16% pour les aggrégés, mais seulement 7% des titulaires alors que ces derniers sont mieux payés en moyenne.
adjoint | aggrege | titulaire | |
---|---|---|---|
femme | 11 | 10 | 18 |
homme | 56 | 54 | 248 |
Plusieurs des variables explicatives potentielles des données college
sont cat/gorielles (echelon
, sexe
, discipline
), les deux dernières étant binaires. Les variables numériques annees
et service
sont fortement corrélées, avec une corrélation linéaire de 0.91.
Exemple 4.7 (Analyse exploratoire et données manquantes) Il convient de vérifier pour les données de Moon et VanEpps (2023) que la description de la collecte coïncide avec la structure. Puisque les personnes qui n’ont pas donné ne remplissent pas le champ pour le montant, ce dernier indique une valeur manquante. Tous les montants des dons sont entre 0.25$ et 25$.
data(MV23_S1, package = "hecedsm")
str(MV23_S1)
#> tibble [869 × 4] (S3: tbl_df/tbl/data.frame)
#> $ before : int [1:869] 0 1 0 1 1 1 1 0 1 0 ...
#> $ donate : int [1:869] 0 0 0 1 1 0 1 0 0 1 ...
#> $ condition: Factor w/ 2 levels "open-ended","quantity": 1 1 1 1 2 2 2 1 1 1 ...
#> $ amount : num [1:869] NA NA NA 10 5 NA 20 NA NA 25 ...
summary(MV23_S1)
#> before donate condition amount
#> Min. :0.000 Min. :0.00 open-ended:407 Min. : 0.2
#> 1st Qu.:0.000 1st Qu.:0.00 quantity :462 1st Qu.: 5.0
#> Median :1.000 Median :1.00 Median :10.0
#> Mean :0.596 Mean :0.73 Mean :10.7
#> 3rd Qu.:1.000 3rd Qu.:1.00 3rd Qu.:15.0
#> Max. :1.000 Max. :1.00 Max. :25.0
#> NA's :1 NA's :235
Si nous incluons amount
comme variable réponse dans un modèle de régression, les 235 observations manquantes seront supprimées par défaut. Cela ne pose pas de problème si nous voulons comparer le montant moyen des personnes qui ont fait un don, mais dans le cas contraire, nous devons transformer les NA
en zéros. La variable donate
ne doit pas être incluse comme variable explicative dans le modèle, car elle permet de prédire exactement les personnes qui n’ont pas donné.
4.1.3 Spécification du modèle pour la moyenne
La première étape d’une analyse consiste à décider quelles variables explicatives doivent être ajoutées à l’équation de la moyenne, et sous quelle forme. Les modèles ne sont que des approximations de la réalité; la section 2.1 de Venables (2000) affirme que, si nous pensons que la véritable fonction moyenne reliant les variables explicatives
Dans un cadre expérimental, où la condition expérimentale est attribué de manière aléatoire, nous pouvons directement comparer les différents traitements et tirer des conclusions causales (puisque toutes les autres choses sont égales en moyenne constantes, toute différence détectable est due en moyenne à notre manipulation). Bien que nous nous abstenions généralement d’inclure d’autres variables explicatives afin de préserver la simplicité du modèle, il peut néanmoins être utile de prendre en compte certaines variables concomitantes qui expliquent une partie de la variabilité afin de filtrer le bruit de fond et d’augmenter la puissance de l’étude. Par exemple, pour les données de Baumann, Seifert-Kessell, et Jones (1992), l’objectif est de comparer les scores moyens en fonction de la méthode d’enseignement, nous inclurions group
. Dans cet exemple, il serait également logique d’inclure le résultat pretest1
en tant qu’élément explicatif pour posttest1
. De cette façon, nous modéliserons la différence moyenne d’amélioration entre le pré-test et le post-test plutôt que le résultat final.
Dans un contexte observationnel, les participants dans différents groupes ont des caractéristiques différentes et nous devons donc tenir compte de ces différences. Les modèles linéaires utilisés en économie et en finance contiennent souvent des variables de contrôle au modèle pour tenir compte des différences potentielles dues aux variables sociodémographiques (âge, revenu, etc.) qui seraient corrélées à l’appartenance aux groupes. Tout test de coefficients ne prendrait en compte que la corrélation entre le résultat
4.2 Interprétation des coefficients
La spécification de la moyenne est
En régression linéaire, le paramètre
Définition 4.1 (Effet marginal) On définit l’effet marginal comme la dérivée première de la moyenne conditionnelle par rapport à
Les variables indicatrices, qui prennent typiquement des valeurs de
Exemple 4.8 (Modèle linéaire avec une seule variable binaire) Considérons par exemple un modèle linéaire pour les données de Moon et VanEpps (2023) qui inclut le montant (amount
) (en dollars, de 0 pour les personnes qui n’ont pas fait de don, jusqu’à 25 dollars).
L’équation du modèle linéaire simple qui inclut la variable binaire condition
est open-ended
) et quantity
). Un modèle linéaire qui ne contient qu’une variable binaire quantity
) est open-ended
et le groupe quantity
. Cette paramétrisation est commode si on veut tester s’il y a une différence moyenne entre les deux groupes, puisque cette hypothèse nulle correspond à
Même si le modèle linéaire définit une droite, cette dernière ne peut être évaluée qu’à
Même s’il est clair que les données sont fortement discrétisées avec beaucoup de doublons et de zéros, l’échantillon a une taille de 869 observations, donc les conclusions quant aux moyennes de groupe seront fiables.
Considérons des variables catégorielles avec factor
. La paramétrisation par défaut des facteurs se fait en termes de contraste de traitement: le niveau de référence du facteur (par défaut, la première valeur dans l’ordre alphanumérique) sera traité comme la catégorie de référence et assimilé à l’ordonnée à l’origine. Le logiciel créera alors un ensemble de
Exemple 4.9 (Codage binaire pour les variables catégorielles) Considérons l’étude de Baumann, Seifert-Kessell, et Jones (1992) et la seule variable group
. Les données sont classées par groupe : les 22 premières observations concernent le groupe DR
, les 22 suivantes le groupe DRTA
et les 22 dernières le groupe TA
. Si nous ajustons un modèle avec groupe
comme variable catégorielle
data(BSJ92, package = "hecedsm")
class(BSJ92$group) # Vérifier que group est un facteur
#> [1] "factor"
levels(BSJ92$group) # première valeur est la catégorie de référence
#> [1] "DR" "DRTA" "TA"
# Imprimer trois lignes de la matrice du modèle
# (trois enfants de groupes différents)
model.matrix(~ group, data = BSJ92)[c(1,23,47),]
#> (Intercept) groupDRTA groupTA
#> 1 1 0 0
#> 23 1 1 0
#> 47 1 0 1
# Comparer avec les niveaux des facteurs
$group[c(1,23,47)]
BSJ92#> [1] DR DRTA TA
#> Levels: DR DRTA TA
Si nous ajustons un modèle avec groupe
comme variable catégorielle, la spécification de la moyenne du modèle est group
est catégorielle avec
Avec la paramétrisation en termes de traitements (option par défaut), on obtient
sigroup=DRTA
et zéro sinon, sigroup=TA
et zéro sinon.
Étant donné que le modèle comprend une ordonnée à l’origine et que le modèle décrit en fin de compte trois moyennes de groupe, nous n’avons besoin que de deux variables supplémentaires. Avec la paramétrisation en termes de traitements, la moyenne du groupe de référence est l’ordonnée à l’origine. Si group
=DR
(référence), les deux variables indicatrices binaires groupDRTA
et groupTA
sont nulles. La moyenne de chaque groupe est
, et .
Ainsi, DRTA
etDR
, et de la même façon
Remarque 4.1 (Contrainte de somme nulle). La paramétrisation discutée ci-dessus, qui est la valeur par défaut de la fonction lm
, n’est pas la seule disponible. Plutôt que de comparer la moyenne de chaque groupe avec celle d’une catégorie de référence, la paramétrisation par défaut pour les modèles d’analyse de la variance est en termes de contraintes de somme nulle pour les coefficients, où l’ordonnée à l’origine est la moyenne équi-pondérée de chaque groupe, et les paramètres
model.matrix(
~ group,
data = BSJ92,
contrasts.arg = list(group = "contr.sum"))
Dans la contrainte de somme nulle, nous obtenons à nouveau deux variables indicatrices, group1
et group2
, ainsi que l’ordonnée à l’origine. La valeur de group1
est group=DR
, group=DRTA
et group=TA
. ous trouvons
En supprimant l’ordonnée à l’origine, on pourrait inclure trois variables indicatrices pour chaque niveau d’un facteur et chaque paramètre correspondrait alors à la moyenne. Ce n’est pas recommandé dans R car le logiciel traite différemment les modèles sans ordonnée à l’origine et certains résultats seront absurdes (par exemple, le coefficient de détermination sera erroné).
Exemple 4.10 (Interprétation des coefficients) On considère un modèle de régression pour les données college
qui inclut le sexe, l’échelon académique, le nombre d’années de service et le domaine d’expertise (appliquée ou théorique).
Le modèle linéaire postulé s’écrit
L’interprétation des coefficients est la suivante:
- L’ordonnée à l’origine
correspond au salaire moyen d’un professeur adjoint (un homme) qui vient de compléter ses études et qui travaille dans un domaine appliqué: on estime ce salaire à dollars. - toutes choses étant égales par ailleurs (même domaine, échelon et années depuis le dernier diplôme), l’écart de salaire entre un homme et un femme est estimé à
dollars. - ceteris paribus, un(e) professeur(e) qui oeuvre dans un domaine théorique gagne
dollars de plus qu’une personne du même sexe dans un domaine appliqué; on estime cette différence à dollars. - ceteris paribus, la différence moyenne de salaire entre professeurs adjoints et aggrégés est estimée à
dollars. - ceteris paribus, la différence moyenne de salaire entre professeurs adjoints et titulaires est de
dollars. - au sein d’un même échelon, chaque année supplémentaire de service mène à une augmentation de salaire annuelle moyenne de
dollars.
Remarque 4.2 (Polynômes). Il n’est pas toujours possible de fixer la valeur des autres colonnes de
L’utilisation de polynôme, plus flexibles, n’est généralement pas recommendée car ces derniers se généralisent mal hors de l’étendue observée des données. L’utilisation de splines avec une pénalité sur les coefficients, avec des modèles additifs, offre plus de flexibilité.
Exemple 4.11 (Modèle quadratique pour les données automobile) Considérons un modèle de régression linéaire pour l’autonomie d’essence en fonction de la puissance du moteur pour différentes voitures dont les caractéristiques sont données dans le jeu de données automobiles
. Le modèle postulé incluant un terme quadratique est
À vue d’oeil, l’ajustement quadratique est bon: nous verrons plus tard à l’aide de test si une simple droite aurait été suffisante. On voit aussi dans la Figure 4.3 que l’autonomie d’essence décroît rapidement quand la puissance croît entre
La représentation graphique du modèle polynomial de degré 2 présenté dans la Figure 4.3 peut sembler contre-intuitive, mais c’est une projection en 2D d’un plan 3D de coordonnées
4.3 Estimation des paramètres
Considérons un échantillon de
4.3.1 Moindres carrés ordinaires
Soit une matrice de modèle
En d’autres mots, les estimateurs des moindres carrés sont la solution du problème d’optimization convexe
Proposition 4.1 (Moindres carrés ordinaires) L’estimateur des moindres carrés ordinaires résoud le problème d’optimisation non-contraint
4.3.2 Maximum de vraisemblance
Nous pourrions également envisager l’estimation du maximum de vraisemblance. Proposition 4.2 montre que, en supposant la normalité des aléas, les estimateurs des moindres carrés de
Proposition 4.2 (Estimation du maximum de vraisemblance du modèle linéaire normal) Le modèle de régression linéaire spécifie que les observations
On déduit que l’estimateur du maximum de vraisemblance est la moyenne des carrés des résidus,
Proposition 4.3 (Matrices d’information pour modèles linéaires normaux.) Les entrées de la matrice d’information observée du modèle linéaire normal sont les suivantes
Remarque 4.3. Si on suppose que les observations sont normales, alors on peut montrer que
4.3.3 Ajustement des modèles linéaires à l’aide d’un logiciel
Bien que nous puissions construire la matrice du modèle nous-mêmes et utiliser la formule des moindres carrés de l’Équation 4.1, les routines numériques implémentées dans les logiciels sont préférables car plus stables. La fonction lm
dans R ajuste les modèles linéaires, tout comme glm
avec les arguments par défaut. Les objets de la classe lm
ont plusieurs méthodes qui vous permettent d’extraire des objets spécifiques des objets lm
. Par exemple, les fonctions coef
, resid
, fitted
, model.matrix
renvoient les estimations des coefficients
data(BSJ92, package = "hecedsm") # charger les données
str(BSJ92) # vérifier que les variables catégorielles sont "factor"
# Ajustement de la régression linéaire
<- lm(posttest1 ~ pretest1 + group,
linmod data = BSJ92)
<- coef(linmod) # coefficients (betas)
est_beta <- vcov(linmod) # matrice de covariance des betas
vcov_beta summary(linmod) # tableau résumé
<- confint(linmod) # IC de Wald pour betas
beta_ic <- fitted(linmod) # valeurs ajustées
y_adj <- resid(linmod) # résidus ordinaires
e
# Vérifier la formule des moindres carrés ordinaires
<- model.matrix(linmod) # matrice du modèle
X <- college$salary
y isTRUE(all.equal(
c(solve(t(X) %*% X) %*% t(X) %*% y),
as.numeric(coef(linmod))
))
La méthode summary
est sans doute la plus utile: elle affiche les estimations des paramètres de la moyenne ainsi que leurs erreurs type, les valeurs lm
utilise l’estimateur sans biais de la variance
4.4 Prédictions
Une fois les estimations des coefficients obtenues, on peut calculer les valeurs ajustées
Si l’on veut prédire la valeur d’une nouvelle observation, disons
Exemple 4.12 (Prédiction pour une régression linéaire simple) Considérons les données de l’Exemple 4.5. On ajuste un modèle de régression linéaire simple avec
La Figure 4.6 montre les bandes d’incertitude ponctuelles pour une simple régression linéaire des données de Sokolova, Krishna, et Döring (2023) en fonction de la proportion
de carton par rapport au plastique, les valeurs les plus élevées indiquant un emballage avec plus de carton superflu. Le modèle ne tient pas compte du fait que notre réponse provient d’une distribution discrète limitée avec des valeurs entières allant de 1 à 7, et que les ratios testés dans l’expérience sont 0 (pas de carton), 0.5, 1 et 2 uniquement. La droite centrale donne la prédiction des individus lorsque nous faisons varier la proportion carton/plastique. En examinant les formules des intervalles de confiance et de prédiction, il est clair que les bandes ne sont pas linéaires (nous considérons la racine carrée d’une fonction qui implique les prédicteurs), mais il n’est pas évident visuellement que l’incertitude augmente au fur et à mesure que l’on s’éloigne de la moyenne des prédicteurs.
Il est plus facile de s’en rendre compte en reproduisant les courbes potentielles qui auraient pu se produire avec des données différentes: la Figure 4.6 montre les nouvelles pentes potentielles générées à partir de la distribution normale asymptotique des estimateurs pef
/proportion
moyenne, et leur potentiel de déviation est d’autant plus élevé que nous nous éloignons de la moyenne dans chaque direction. Les intervalles de prédiction (gris pâle) sont très larges et couvrent essentiellement l’ensemble des valeurs potentielles de l’échelle de Likert sur la perception du respect de l’environnement, à l’exception de quelques observations. En revanche, les intervalles de confiance pour la moyenne sont assez étroits, en raison de la taille importante de l’échantillon. On constate également que les courbes s’en écartent peu.
Dans R, la fonction générique predict
prend comme arguments un modèle et une nouvelle base de données newdata
contenant un tableau avec la même structure que les données qui ont servi à l’ajustement du modèle (à minima, les colonnes de variables explicatives utilisées dans le modèle).
data(SKD23_S2A, package = "hecedsm") # charger les données
<- lm(pef ~ proportion, data = SKD23_S2A) # régression linéaire simple
lm_simple predict(lm_simple,
newdata = data.frame(proportion = c(0, 0.5, 1, 2)),
interval = "prediction") # intervalles de prédiction
predict(lm_simple,
newdata = data.frame(proportion = c(0, 0.5, 1, 2)),
interval = "confidence") # IC de confiance pour la moyenne
proportion |
prédiction | borne inf. | borne sup. |
---|---|---|---|
0.0 | 2.41 | -0.168 | 4.98 |
0.5 | 2.67 | 0.097 | 5.24 |
1.0 | 2.93 | 0.361 | 5.51 |
2.0 | 3.46 | 0.884 | 6.04 |
moyenne | borne inf. | borne sup. |
---|---|---|
2.41 | 2.27 | 2.55 |
2.67 | 2.57 | 2.77 |
2.93 | 2.84 | 3.02 |
3.46 | 3.30 | 3.62 |
4.5 Tests d’hypothèses
Les tests d’hypothèses dans les modèles linéaires et d’analyse de la variance suivent la procédure usuelle: nous comparons deux modèles emboîtés, dont l’un (le modèle nul) est une simplification d’un modèle plus complexe (modèle alternatif) obtenu en imposant des restrictions sur les coefficients de la moyenne.
Les tests de restrictions pour les composantes de
Dans un contexte inférentiel, il est souvent important de tester si l’effet d’une variable explicative est significatif : si
Proposition 4.4 (Tests de Wald en régression linéaire) Rappelons que la statistique du test de Wald pour l’hypothèse summary
. Outre les estimations des coefficients, il est possible d’obtenir des intervalles de confiance basés sur Wald pour
Exemple 4.13 Considérons les données de Exemple 4.5. Si nous ajustons le modèle de régression linéaire simple, nous pouvons extraire les valeurs -proportion=0
ne peut être nulle. Le coefficient de proportion
suggère une tendance de 0.5 point par unité de ratio, et il est significativement différent de zéro, ce qui indique que le score pef
change avec le ratio carton/plastique.
# tests-t (Wald) pour beta=0 avec valeurs-p
summary(lm_simple)$coefficients
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.407 0.0723 33.31 2.56e-153
#> proportion 0.526 0.0618 8.51 8.40e-17
confint(lm_simple) # intervalles de confiance pour betas
#> 2.5 % 97.5 %
#> (Intercept) 2.266 2.549
#> proportion 0.405 0.648
Pour les variables catégorielles à plus de deux niveaux, tester si
Proposition 4.5 (Tests-F pour comparaison de modèles emboîtés) Considérons le modèle linéaire complet qui contient
La statistique
Quand la
Remarque 4.4 (Tests F versus test du rapport de vraisemblance). Pour la régression linéaire normale, le test du rapport de vraisemblance pour comparer les modèles
4.5.1 Contrastes
Supposons que nous effectuions une analyse de la variance et que le test
La question scientifique qui a justifié l’expérience peut conduire à un ensemble spécifique d’hypothèses, qui peuvent être formulées par les chercheurs comme des comparaisons entre les moyennes de différents sous-groupes. Nous pouvons normalement les exprimer sous la forme de contrastes. Si le test global
Les contrastes encodent la question de recherche : si
Si nous avons
4.5.2 Exemples de tests
Exemple 4.14 (Test du montant des dons) Considérons l’Exemple 4.8, dans lequel nous testons les différences entre les montants libres (open-ended
) et les montants suggérés (quantity
). Le test qui nous intéresse est
data("MV23_S1", package = "hecedsm")
<- MV23_S1 |>
MV23_S1 ::mutate(amount2 = ifelse(is.na(amount), 0, amount))
dplyr<- lm(amount2 ~ condition, data = MV23_S1)
linmod_MV23 # Test Wald avec coefficients
summary(linmod_MV23)
#>
#> Call:
#> lm(formula = amount2 ~ condition, data = MV23_S1)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -8.70 -6.77 -1.77 3.23 18.23
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.771 0.377 17.95 <2e-16 ***
#> conditionquantity 1.929 0.517 3.73 0.0002 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 7.61 on 867 degrees of freedom
#> Multiple R-squared: 0.0158, Adjusted R-squared: 0.0147
#> F-statistic: 13.9 on 1 and 867 DF, p-value: 0.000205
# ANOVA avec tests F
anova(linmod_MV23)
#> Analysis of Variance Table
#>
#> Response: amount2
#> Df Sum Sq Mean Sq F value Pr(>F)
#> condition 1 805 805 13.9 0.0002 ***
#> Residuals 867 50214 58
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Moyennes marginales
<- emmeans::emmeans(linmod_MV23, spec = "condition"))
(emm #> condition emmean SE df lower.CL upper.CL
#> open-ended 6.77 0.377 867 6.03 7.51
#> quantity 8.70 0.354 867 8.01 9.40
#>
#> Confidence level used: 0.95
|> emmeans::contrast(method = "pairwise") # vecteur de contraste (1,-1)
emm #> contrast estimate SE df t.ratio p.value
#> (open-ended) - quantity -1.93 0.517 867 -3.730 0.0002
Exemple 4.15 (Tests et contrastes pour les méthodes de compréhension de la lecture) Nous examinons maintenant les tests pour l’Exemple 4.2 et l’Exemple 4.9, avec une covariable en plus. L’objectif de Baumann, Seifert-Kessell, et Jones (1992) était de faire une comparaison particulière entre des groupes de traitement. Selon le résumé de l’article:
Les analyses quantitatives principales comportaient deux contrastes orthogonaux planifiés: l’effet de l’enseignement (TA + DRTA vs. 2 x DR) et l’intensité de l’enseignement (TA vs. DRTA).
Avec un modèle pré-post, nous allons comparer les moyennes pour une valeur commune de pretest1
, ci-dessous la moyenne globale du score pretest1
.
library(emmeans) # moyennes marginales
data(BSJ92, package = "hecedsm")
<- lm(posttest1 ~ group + pretest1,
mod_post data = BSJ92)
<- lm(posttest1 ~ pretest1,
mod_post0 data = BSJ92)
anova(mod_post0, mod_post) # tests F
#> Analysis of Variance Table
#>
#> Model 1: posttest1 ~ pretest1
#> Model 2: posttest1 ~ group + pretest1
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 64 509
#> 2 62 365 2 143 12.2 0.000035 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<- emmeans(object = mod_post,
emmeans_post specs = "group")
Le résultat du tableau d’analyse de la variance montre qu’il y a bien des différences entre les groupes. On peut donc s’intéresser aux moyennes marginales estimées, qui sont la moyenne de chaque groupe.
Les deux hypothèses et contrastes de Baumann, Seifert-Kessell, et Jones (1992) sont
# Identifier l'ordre de niveau du facteur
with(BSJ92, levels(group))
#> [1] "DR" "DRTA" "TA"
# DR, DRTA, TA (alphabetical)
<- list(
contrastes_list # Contrastes: combo linéaire de moyennes,
# la somme des coefficients doit être nulle
"C1: moy(DRTA+TA) vs DR" = c(-1, 0.5, 0.5),
"C2: DRTA vs TA" = c(0, 1, -1)
)<-
contrastes_post contrast(object = emmeans_post,
method = contrastes_list)
<- summary(contrastes_post) contrastes_summary_post
Nous pouvons examiner ces différences: puisque DRTA
contre TA
est une différence par paire, nous aurions pu obtenir la statistique pairs(emmeans_post)
.
Quelle est la conclusion de notre analyse des contrastes? Il semble que les méthodes impliquant la réflexion à haute voix aient un impact important sur la compréhension de la lecture par rapport à la seule lecture dirigée. Les preuves ne sont pas aussi solides lorsque nous comparons la méthode qui combine la lecture dirigée, l’activité de réflexion et la réflexion à haute voix, mais la différence est néanmoins significative à niveau 5%.
# Extraire les coefficients et les erreurs-type
<- coefficients(mod_post)['pretest1']
beta_pre <- sqrt(c(vcov(mod_post)['pretest1', 'pretest1']))
se_pre <- (beta_pre - 1)/se_pre # test de Wald directionnel
wald # Valeur-p basée sur la référence nulle Student-t avec n-p-1 ddl
<- 2*pt(abs(wald), df = mod_post$df.residual, lower.tail = FALSE)
pval # Comparaison de modèles emboîtés avec appel à 'anova'
<- lm(posttest1 ~ offset(pretest1) + group, data = BSJ92)
mod0 # Le décalage (`offset`) fixe le terme, ce qui équivaut à un coefficient de 1.
<- anova(mod0, mod_post) aov_tab
Une autre hypothèse potentielle intéressante consiste à tester si le coefficient de pretest1
est égal à l’unité. Cela équivaut à l’hypothèse anova
, qui donne une statistique de test de
Exemple 4.16 (Tests et contrastes pour l’effet de l’emballage carton sur la perception) Soit
data(SKD23_S2A, package = "hecedsm")
<- lm(pef ~ proportion, data = SKD23_S2A)
linmod coef(linmod) # extraire coefficients
#> (Intercept) proportion
#> 2.407 0.526
# ANOVA à un facteur
<- lm(pef ~ factor(proportion),
anovamod data = SKD23_S2A)
# Comparer les deux modèles emboîtés
anova(linmod, anovamod) # est-ce que l'effet est linéaire?
#> Analysis of Variance Table
#>
#> Model 1: pef ~ proportion
#> Model 2: pef ~ factor(proportion)
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 800 1373
#> 2 798 1343 2 29.3 8.69 0.00018 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test avec code alternatif (poids pour chaque coefficient)
::linearHypothesis(model = anovamod,
carhypothesis = rbind(c(0, -2, 1, 0),
c(0, 0, -2, 1)))
#> Linear hypothesis test
#>
#> Hypothesis:
#> - 2 factor(proportion)0.5 + factor(proportion)1 = 0
#> - 2 factor(proportion)1 + factor(proportion)2 = 0
#>
#> Model 1: restricted model
#> Model 2: pef ~ factor(proportion)
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 800 1373
#> 2 798 1343 2 29.3 8.69 0.00018 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le résultat montre que les tests anova
.
Les auteurs souhaitaient comparer zéro carton avec d’autres choix: nous nous intéressons aux différences par paire, mais uniquement par rapport à la référence
<- anovamod |>
moymarg ::emmeans(specs = "proportion") # moyennes de groupes
emmeans<- list( # liste de vecteurs de contrastes
contrastlist refvsdemi = c(1, -1, 0, 0),
refvsun = c(1, 0, -1, 0),
refvsdeux = c(1, 0, 0, -1))
# calculer différences relativement à la référence
|> emmeans::contrast(method = contrastlist)
moymarg #> contrast estimate SE df t.ratio p.value
#> refvsdemi -0.749 0.131 798 -5.710 <.0001
#> refvsun -0.901 0.131 798 -6.890 <.0001
#> refvsdeux -1.182 0.129 798 -9.200 <.0001
Les moyennes des groupes rapportées dans le Tableau 4.6 correspondent à celles indiquées par les auteurs dans l’article. Elles suggèrent que la perception du respect de l’environnement augmente avec la quantité de carton utilisée dans l’emballage. Nous avons pu ajuster un modèle de régression simple pour évaluer le changement moyen, en traitant la proportion comme une variable explicative continue. La pente estimée pour le changement du score PEF, qui va de 1 à 7 par incréments de 0.25, est 0.53 point par rapport au carton/plastique. Il y a cependant de fortes indications, compte tenu des données, que le changement n’est pas tout à fait linéaire, puisque l’ajustement du modèle de régression linéaire est significativement plus mauvais que le modèle linéaire correspondant.
Toutes les différences dans le Tableau 4.7 sont significatives et positives, conformément à l’hypothèse des chercheurs.
Exemple 4.17 (Tester la discrimination salariale dans l’enseignement supérieur) Considérons l’exemple des données college
et le modèle linéaire associé avec echelon
, sexe
, années de service
et domaine
comme variables explicatives.
data(college, package = "hecmodstat")
<- lm(salaire ~ sexe + domaine + echelon + service, data = college)
mod1_college <- lm(salaire ~ domaine + echelon + service, data = college)
mod0_college # F-test avec "anova" comparant les modèles emboîtés
<- anova(mod0_college, mod1_college)
aov_tab_college # Test t de Wald
<- summary(mod1_college)$coefficients[2,]
wald_college # Test du rapport de vraisemblance avec approx khi-deux
<- pchisq(q = as.numeric(2*(logLik(mod1_college) - logLik(mod0_college))),
pval_lrt df = 1, lower.tail = FALSE)
Le seul test qui nous intéresse ici est
term | estimation | erreur-type | stat de Wald | valeur-p |
---|---|---|---|---|
(Intercept) | 86.596 | 2.96 | 29.25 | < 0.001 |
sexe [femme] | -4.771 | 3.878 | -1.23 | 0.22 |
domaine [théorique] | -13.473 | 2.315 | -5.82 | < 0.001 |
échelon [agrégé] | 14.56 | 4.098 | 3.55 | < 0.001 |
échelon [titulaire] | 49.16 | 3.834 | 12.82 | < 0.001 |
service | -0.089 | 0.112 | -0.8 | 0.43 |
4.6 Plans factoriels et interactions
Le modèle additif pour la moyenne spécifie que l’effet marginal d’une variable ( y compris pour les variables catégorielles) est indépendant des autres. Nous pouvons souhaiter assouplir cette hypothèse en incluant des termes d’interaction.
Définition 4.2 (Interaction) On parle d’interaction lorsque des combinaisons de variables explicatives affectent la variable réponse différemment que lorsqu’elles sont considérées individuellement. Si
On s’attarde au cas où au moins une des variables est catégorielle (facteur).
Exemple 4.18 (Primes d’assurances et interactions) On considère la relation entre fumeur
et l’indice de masse corporel pour la détermination de primes d’assurance. Les fumeurs dont l’indice de masse corporelle (IMC) est égal ou supérieur à 30 paient une prime élevée, mais il semble que le montant de la prime augmente de façon linéaire en fonction de l’IMC. Cette tarification ne semble pas s’appliquer aux non-fumeurs.
Exemple 4.19 (Intention d’achat) On considère un exemple avec des données bidons interaction
. Le modèle additif (sans interaction) a pour moyenne
L’effet de la variable continue fixation
est identique pour les deux sexes. De même, l’effet de la variable binaire est supposé être le même pour toutes les valeurs possibles de la variable continue. Nous pouvons le voir sur le graphique, car la différence entre les lignes représente l’effet de
Pour ajuster une pente différente par sexe, on crée une nouvelle variable égale au produit
Selon la valeur de
est l’intention d’achat moyenne lorsque le temps de fixation est nul pour les hommes, est la différence d’ordonnée à l’origine entre les femmes et les hommes (différence d’intention d’achat moyenne entre femmes et hommes quand le temps de fixation est nul), est l’augmentation unitaire de l’intention d’achat par seconde de fixation pour les hommes, est la différence de pente entre les femmes et les hommes (différence d’intention d’achat moyenne femmes vs hommes pour une augmentation d’une seconde de fixation).
Tester la significativité de l’interaction revient à vérifier si
data(interaction, package = "hecmodstat")
# Pour spécifier une interaction, utiliser :
<- lm(intention ~ sexe + fixation + sexe:fixation,
mod data = interaction)
# Un raccourci est sexe*fixation, qui donne la même chose
summary(mod)$coefficients
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.741 0.282 9.73 1.02e-16
#> sexe 1.312 0.380 3.45 7.74e-04
#> fixation 0.504 0.153 3.29 1.33e-03
#> sexe:fixation 2.135 0.200 10.69 5.61e-19
Le modèle avec interaction est significativement meilleur, ce qui signifie que l’effet du temps de fixation sur l’intention d’achat varie en fonction du sexe.
Remarque 4.5 (Principe de marginalité). Tous les termes d’ordre inférieurs devraient être inclus si l’interaction est présente.
Par exemple, on ne retirera pas
Comme le choix de catégorie de référence est arbitraire, changer la variable indicatrice pour
Le concept d’interaction se généralise à des variables catégorielles avec plus de deux niveaux. Dans ce cas, on doit considérer la statistique
Définition 4.3 (Analyse de variance) Une analyse de variance est un modèle linéaire dans lequel la moyenne est une fonction de variables explicatives catégorielles. Si nous disposons de données pour toutes les combinaisons différentes de facteurs, les facteurs sont croisés et nous pouvons envisager d’inclure leurs interactions.
Considérons un modèle d’analyse de variance à deux facteurs, disons
est le e mesure du e niveau du facteur et du e niveau du facteur , disons est la moyenne théorique du sous-groupe sont des termes d’aléas indépendants de moyenne zéro et d’écart-type .
Dans un devis factoriel complet avec interaction, on peut écrire l’espérance de la variable réponse
Nous pouvons l’exprimer de manière équivalente en termes d’ordonnée à l’origine, d’effets principaux de l’une ou l’autre variable et de termes d’interaction. Le modèle additif, sans interaction, a une moyenne pour la cellule
Nous pouvons envisager des simplifications de modèle de bas en haut. La suppression de l’interaction conduit à un modèle avec
Bien que des tests formels soient nécessaires pour vérifier les interactions, le concept peut être mieux compris en examinant des graphiques.
Définition 4.4 (Diagramme d’interaction) Nous pouvons essayer de détecter visuellement les interactions en traçant la réponse (moyenne) en fonction de l’une des covariables, en utilisant ce que l’on appelle un diagramme d’interaction. Lorsqu’il y a plus de deux variables catégorielles, nous pouvons utiliser des couleurs, des symboles ou des panneaux pour représenter les catégories. L’absence d’interaction dans ces diagrammes implique des lignes parallèles, mais il faut tenir compte de l’incertitude.
Définition 4.5 (Effets simples et effets principaux) Lorsqu’il n’y a pas d’interactions, il est logique de faire abstraction d’une ou plusieurs variables et de considérer les effets marginaux, obtenus en regroupant les données des facteurs omis, en calculant la moyenne équipondérée des sous-groupes. Supposons, par exemple, que nous soyons intéressés par la comparaison des niveaux de
- effets simples: différences entre les niveaux d’un élément dans une combinaison fixe d’autres éléments. Les effets simples consistent à comparer les moyennes des cellules dans une ligne ou une colonne donnée.
- effets principaux : différences par rapport à la moyenne pour chaque condition d’un facteur. Les effets principaux sont des moyennes de lignes ou de colonnes.
Exemple 4.20 (Effet psychologique d’emprunt) L’étude complémentaire 5 de Sharma, Tully, et Cryder (2021) vérifie la perception psychologique de l’emprunt en fonction du terme employé pour désigner le montant monétaire. Les auteurs ont effectué une comparaison deux par deux inter-sujets (ANOVA à deux facteurs) en faisant varier le type de dette (si l’argent était proposé sous forme de crédit (credit
) ou prêt (loan
), de même que le type d’achat pour lequel l’argent serait utilisé, soit dépenses discrétionnaires (discretionary
) ou pour combler à des besoins essentiels (need
). La réponse est la moyenne de (a) la probabilité et de (b) l’intérêt pour le produit, toutes deux mesurés à l’aide d’une échelle de Likert allant de 1 à 9.
Le modèle pour la moyenne avec interaction peut être écrit en utilisant la paramétrisation usuelle comme suit
On peut calculer la moyenne de chaque groupe et déduire l’interprétation des coefficients:
pourpurchase=discretionnary
etdebttype=credit
pourpurchase=need
etdebttype=credit
pourpurchase=discretionnary
etdebttype=loan
pourpurchase=need
etdebttype=loan
Ainsi,
Sharma, Tully, et Cryder (2021) a ajusté un modèle avec deux facteurs, chacun avec deux niveaux, et leur interaction. Comme il y a une moyenne globale et deux effets principaux (différence supplémentaire de moyenne pour les deux facteurs debttype
et purchase
), l’interaction a un degré de liberté puisque nous passons d’un modèle avec trois paramètres à un modèle qui a une moyenne différente pour chacun des quatre sous-groupes.
La raison pour laquelle on teste d’abord la présence d’interaction est que, si l’effet d’un facteur dépend du niveau de l’autre, comme le montre Figure 4.9, il faut alors comparer l’étiquette du type de dette séparément pour chaque type d’achat et vice-versa à l’aide d’effets simples. Si l’interaction n’est pas significative, nous pouvons regrouper les observations pour obtenir une ANOVA à un facteur, ce qui donne les comparaisons marginales avec les effets principaux.
L’ajustement du modèle incluant l’interaction entre les facteurs garantit que nous conservons l’hypothèse d’additivité et que nos conclusions ne sont pas trompeuses: le prix à payer est l’estimation de paramètres moyens supplémentaires, ce qui n’est pas un problème si vous recueillez suffisamment de données, mais peut être critique lorsque la collecte de données est extrêmement coûteuse et que seules quelques observations sont disponibles pour chaque sous-groupe.
Dans R, on inclut les deux facteurs et leurs interactions avec reponse ~ facteurA * facteurB
, le symbole *
indiquant que les deux interagissent, un raccourci pour facteuA + facteurB + facteurA:facteurB
; dans un modèle avec seulement les effets principaux, on sépare les facteurs par un +
pour obtenir le modèle additif.
# Données de l'étude supp. 5
# de Sharma, Tully, et Cryder (2021)
data(STC21_SS5, package = "hecedsm")
# La fonction 'aov' sert à ajuste des ANOVA
# Équivalent à "lm" avec variables catégorielles, contrasts somme nulle
<- aov(
modlin_STC21 ~ purchase*debttype,
likelihood data = STC21_SS5)
# Calculer le décompte par sous-catégorie (données débalancées)
xtabs(~ purchase + debttype, data = STC21_SS5)
#> debttype
#> purchase credit loan
#> discretionary 392 359
#> need 361 389
# Calcul de la moyenne globale/lignes/colonnes/cellules
<- model.tables(modlin_STC21, type = "means")
moy_groupes
# Tableau d'ANOVA avec effets de type 2
::Anova(modlin_STC21, type = 2)
car#> Anova Table (Type II tests)
#>
#> Response: likelihood
#> Sum Sq Df F value Pr(>F)
#> purchase 752 1 98.21 < 2e-16 ***
#> debttype 92 1 12.04 0.00054 ***
#> purchase:debttype 14 1 1.79 0.18171
#> Residuals 11467 1497
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L’interaction n’étant pas significative, nous pouvons n’interpréter que l’effet principal de la fixation.
Cette différence de moyenne conditionnelle est appelée effet marginal, car elle est obtenue en calculant la moyenne de toutes les sous-catégories pour un même niveau du facteur. Le modèle estime cependant la variance sur la base des résidus du modèle d’interaction complet avec quatre moyennes de cellules, et diffère donc de celui obtenu en exécutant (à tort) un modèle avec seulement purchase
comme variable explicative.
Dans le tableau d’analyse de la variance, nous nous concentrons exclusivement sur la dernière ligne avec la somme des carrés pour l’interaction purchase:debttype
. La statistique
On peut donc regrouper les données et étudier uniquement l’effet du libellé (prêt loan
ou credit
) en combinant les données pour les types d’achats, une des comparaisons planifiées des annexes en ligne. Pour ce faire, on utilise la fonction emmeans
du paquet éponyme en spécifiant le nom du ou des facteurs d’intérêts (ceux que l’on veut conserver) avec l’argument specs
. Par défaut, on calcule les moyennes marginales estimées, l’argument contr = "pairwise"
indique que l’on veut en plus les différences deux à deux, ici le seul contraste possible pour des différences entre deux facteurs.
# Comparisons deux à deux au sein des niveaux de "purchase"
# Effets simples
::emmeans(modlin_STC21,
emmeansspecs = "purchase",
contr = "pairwise")
#> $emmeans
#> purchase emmean SE df lower.CL upper.CL
#> discretionary 4.17 0.101 1497 3.97 4.36
#> need 5.58 0.101 1497 5.39 5.78
#>
#> Results are averaged over the levels of: debttype
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> discretionary - need -1.42 0.143 1497 -9.910 <.0001
#>
#> Results are averaged over the levels of: debttype
# Diagramme d'interaction
::emmip(modlin_STC21,
emmeans~ debttype,
purchase CIs = TRUE) +
theme_minimal()
Remarque 4.6 (Décomposition de sommes des carrés). Il existe différentes décompositions de la somme des carrés (type I, II et III) pour la comparaison des modèles emboîtés dans les tableaux d’analyse de la variance.
Ces décompositions testent différents modèles à l’aide des statistiques
Tableau 4.8 montre les différentes sommes des erreurs quadratiques des modèles, avec les termes entre parenthèses indiquant quels termes sont inclus (
La décomposition de type I, la valeur par défaut du générique anova
, utilise l’ordre dans lequel les termes sont spécifiés, disons
La décomposition de type II considère les termes de même niveau dans la hiérarchie, de sorte que les tests pour les effets principaux sont
La décomposition de type III, popularisée par SAS et souvent choisie par défaut dans les logiciels, prend en compte tous les autres termes, et testerait donc les effets principaux comme
Les trois méthodes donnent le même résultat et la même comparaison pour le dernier niveau avec l’interaction.
Toutes les discussions relatives à une ANOVA à deux voies s’appliquent à des plans d’expérience à
Exemple 4.21 (Perception d’appropriation culturelle par idéologie politique) On considère un modèle d’ANOVA à trois facteurs de Lin et al. (2024). Leur étude 4 s’intéresse à l’appropriation culturelle pour une recette de soul food, un courant afro-américain, parue dans le livre du “Chef Dax”. Les auteurs manipulent l’ethnicité du chef, afro-Américain ou pas, et la façon dont la recette a été obtenue (furtivement, en demandant la permission ou sans mention dans le cas du groupe contrôle). Les auteurs ont postulé que la perception de l’appropriation varie selon l’idéologie politique. L’étude utilise un devis expérimental
data(LKUK24_S4, package = "hecedsm")
# Vérifier la répartition d'observations en sous-groupes
xtabs(~politideo + chefdax + brandaction,
data = LKUK24_S4)
#> , , brandaction = peeking
#>
#> chefdax
#> politideo not black black
#> conservative 33 36
#> liberal 87 84
#>
#> , , brandaction = permission
#>
#> chefdax
#> politideo not black black
#> conservative 42 34
#> liberal 77 84
#>
#> , , brandaction = control
#>
#> chefdax
#> politideo not black black
#> conservative 38 32
#> liberal 79 85
# Les facteurs sont croisés et il y a des réplications
# On ajuste un modèle d'ANOVA à trois facteurs (avec toutes les interactions)
<- lm(appropriation ~ politideo * chefdax * brandaction,
mod data = LKUK24_S4)
# Calculer les estimations des moyennes marginales
# pour chacun des 12 sous-groupes
<- emmeans(mod,
emm specs = c("chefdax", "brandaction", "politideo"))
# Créer un diagramme d'interaction
emmip(object = emm,
formula = brandaction ~ chefdax | politideo,
CIs = TRUE) +
::scale_color_met_d(name = "Hiroshige") +
MetBrewerlabs(y = "prédicteur linéaire",
x = "niveaux de chefdax")
# Tableau d'ANOVA (type 2)
<- car::Anova(mod, type = 2) anova_tab
Pour l’ANOVA à
terme | somme des carrés | ddl | stat | valeur-p |
---|---|---|---|---|
politideo | 48.49 | 1 | 21.35 | <0.001 |
chefdax | 473.72 | 1 | 208.61 | <0.001 |
brandaction | 34.24 | 2 | 7.54 | <0.001 |
politideo:chefdax | 65.00 | 1 | 28.63 | <0.001 |
politideo:brandaction | 1.56 | 2 | 0.34 | 0.71 |
chefdax:brandaction | 0.62 | 2 | 0.14 | 0.87 |
politideo:chefdax:brandaction | 0.66 | 2 | 0.15 | 0.86 |
Residuals | 1587.33 | 699 |
Si on considère le Tableau 4.9, nous constatons qu’il n’y a pas d’interaction à trois voies et, si l’on omet cette dernière et que l’on se concentre sur les niveaux inférieurs, une seule interaction à deux voies entre l’idéologie politique et la race du chef Dax.Nous ne pouvons pas interpréter la valeur brandaction
, mais nous pouvons examiner les moyennes marginales.
Sur la base des données, nous réduirons les données à une ANOVA à une voie comparant les trois niveaux de brandaction
et à une ANOVA à deux facteurs, chefdax
et politideo
. Les résultats sont obtenus en calculant la moyenne pour le facteur manquant, mais en estimant l’écart-type du modèle complet.
Nous souhaitons comparer la perception de la race du chef Dax (noir ou non), car la cuisine soul est plus susceptible d’être associée à l’appropriation culturelle si le chef Dax n’est pas noir. Nous procédons avec emmeans
en calculant les moyennes marginales séparément pour chacune des quatre sous-catégories, mais nous comparons la race du chef Dax séparément pour les libéraux et les conservateurs en raison de la présence de l’interaction.
data(LKUK24_S4, package = "hecedsm")
library(emmeans)
<- lm(appropriation ~ politideo * chefdax * brandaction,
mod data = LKUK24_S4)
# Moyennes marginales pour
# idéologie politique/ethnicité du Chef Dax
# Calculer les effets simples séparément
# pour chaque idéologie politique
emmeans(mod,
specs = "chefdax",
by = "politideo",
contrast = "pairwise")
#> politideo = conservative:
#> chefdax emmean SE df lower.CL upper.CL
#> not black 2.38 0.1425 699 2.11 2.66
#> black 1.68 0.1494 699 1.38 1.97
#>
#> politideo = liberal:
#> chefdax emmean SE df lower.CL upper.CL
#> not black 3.60 0.0968 699 3.41 3.79
#> black 1.57 0.0947 699 1.38 1.75
#>
#> Results are averaged over the levels of: brandaction
#> Confidence level used: 0.95
Nous constatons que les libéraux sont beaucoup plus susceptibles de considérer le livre de cuisine du chef Dax comme un exemple d’appropriation culturelle s’il n’est pas noir; il y a peu de preuves d’une différence entre les conservateurs et les libéraux lorsque le chef Dax est noir. On peut calculer les effets marginaux pour idéologie (afro-Américain ou pas). Les deux différences sont statistiquement significatives, mais la différence est beaucoup plus marquée pour les répondants de gauche.
Pour brandaction
, nous supposons que les participants verront le fait de copier furtivement moins favorablement que si le chef Dax demandait l’autorisation de publier la recette. Il est difficile de connaître l’effet du groupe contrôle, car on ne mentionne pas comment la recette a été acquise.
# Effet marginal (effet principal) pour brandaction
<- emmeans(mod, specs = c("brandaction"))
emm_brand
emm_brand#> brandaction emmean SE df lower.CL upper.CL
#> peeking 2.56 0.107 699 2.35 2.77
#> permission 2.29 0.105 699 2.09 2.50
#> control 2.07 0.108 699 1.86 2.28
#>
#> Results are averaged over the levels of: politideo, chefdax
#> Confidence level used: 0.95
# Test F conjoint pour l'effet principale de brandaction
|> pairs() |> joint_tests()
emm_brand #> model term df1 df2 F.ratio p.value
#> contrast 2 699 5.090 0.0064
Un test
4.7 Géométrie des moindres carrés
Remarque 4.9.
Remarque 4.9 (Géométrie). Le vecteur de valeurs ajustées
Corollaire 4.1 (Orthogonalité des résidus et des valeurs ajustées) Une conséquence directe de ces résultats est le fait que la corrélation linéaire entre
Corollaire 4.2 (Moyenne des résidus nulle) Puisque le produit scalaire est zéro, la moyenne de
<- lm(salaire ~ sexe + echelon + domaine + service, data = college)
mod # Corrélation nulle
with(college, cor(resid(mod), service))
#> [1] 1.26e-17
cor(resid(mod), fitted(mod))
#> [1] 2.32e-17
# Moyenne des résidus ordinaires est nulle
mean(resid(mod))
#> [1] -5.6e-16
Puisque les aléas avaient moyenne théorique de zéro, on veut forcer les résidus ordinaires à avoir une moyenne empirique de zéro en incluant l’ordonnée à l’origine.
Remarque 4.8 (Invariance). Une conséquence directe de l’expression de l’estimateur des MCO en terme de matrice de projection est que les valeurs ajustées
<- lm(salaire ~ sexe + echelon + service, data = college)
modA <- lm(salaire ~ 0 + sexe + echelon + service, # Enlever l'ordonnée à l'origine
modB data = college |>
::mutate(service = scale(service)), # Centrer-réduire une variable
dplyrcontrasts = list(echelon = contr.sum)) # changer la paramétrisation
head(model.matrix(modA), n = 3L)
#> (Intercept) sexefemme echelonaggrege echelontitulaire service
#> 1 1 0 0 1 18
#> 2 1 0 0 1 16
#> 3 1 0 0 0 3
head(model.matrix(modB), n = 3L)
#> sexehomme sexefemme echelon1 echelon2 service
#> 1 1 0 -1 -1 0.0296
#> 2 1 0 -1 -1 -0.1241
#> 3 1 0 1 0 -1.1237
# Invariance du modèle
isTRUE(all.equal(fitted(modA), fitted(modB)))
#> [1] TRUE
La valeur de
4.7.1 Résidus
Les résidus sont les prédictions des aléas
Toutes les observations ne contribuent pas de la même manière à l’ajustement de l’hyperplan ajusté. La géométrie des moindres carrés montre que les résidus sont orthogonaux aux valeurs ajustées, et que
Si les aléas sont indépendants et homoscédastiques, les résidus ordinaires
Nous concluons donc que les résidus ordinaires n’ont pas tous le même écart-type et qu’ils ne sont pas indépendants. Ceci est problématique, car nous ne pouvons pas faire de comparaisons de leurs lois: les points ayant un faible effet de levier s’écartent davantage du modèle ajusté que les autres. Pour pallier ce problème, nous pouvons normaliser les résidus de façon à ce que chacun ait la même variance sous l’hypothèse nulle d’erreurs homoscédastiques indépendantes — les termes d’effet de levier
La seule question qui subsiste est celle de l’estimation de la variance. Si nous utilisons la rstudent
.
Quand utiliser quels résidus? Par construction, le vecteur des résidus ordinaires
Bien que les résidus studentisés externes
Définition 4.6 (Coefficient de détermination) Lorsque nous spécifions un modèle, les aléas
La somme totale des carrés, définie comme la somme des carrés des résidus du modèle à ordonnée à l’origine uniquement, sert de comparaison — le modèle le plus simple que nous puissions trouver impliquerait chaque observation par la moyenne de l’échantillon de la réponse, ce qui donne la variance expliquée
Il est important de noter que le
En outre, il est possible de gonfler la valeur de
4.7.2 Colinéarité
Le postulat de linéarité peut être interprétée au sens large comme signifiant que toutes les covariables pertinentes ont été incluses et que leur effet est correctement spécifié dans l’équation de la moyenne. L’ajout de covariables superflues à un modèle a un impact limité: si la corrélation (partielle) entre un vecteur colonne
Dans les études observationnelles, il est néanmoins préférable d’inclure davantage de variables que d’oublier des variables explicatives clés: si nous omettons une variable prédictive importante, son effet peut être capturé par d’autres variables confondantes, qui sont corréléws avec à la fois les variables omises et la réponse et qui causent les deux. Par exemple, le modèle linéaire simple (ou le test
Un modèle linéaire n’est pas un [modèle causal] (https://xkcd.com/552/): il ne fait que capturer la corrélation linéaire entre une variable explicative et la réponse. Lorsqu’il y a plus d’une variable explicative, l’effet de
L’un des inconvénients de la colinéarité est la diminution de la précision des estimateurs de paramètres. En présence de variables explicatives colinéaires, de nombreuses combinaisons linéaires des covariables représentent presque aussi bien la réponse. En raison du manque (ou presque) d’identifiabilité, les coefficients estimés deviennent numériquement instables, ce qui entraîne une augmentation des erreurs-type des paramètres. Les valeurs prédites ou ajustées ne sont pas affectées. En général, les coefficients de régression peuvent changer radicalement lorsque de nouvelles observations sont incluses dans le modèle, ou lorsque nous incluons ou supprimons des variables explicatives. Les coefficients
Le diagramme de régression partielle montre à l’aide d’un nuage de points la relation entre la réponse
Une idée similaire peut être utilisée pour voir quelle part de
Que doit-on faire s’il y a de la colinéarité ? Si l’objectif de l’étude est de développer un modèle prédictif et que nous ne sommes pas intéressés par les paramètres eux-mêmes, alors nous n’avons rien à faire. La colinéarité n’est pas un problème pour le modèle global: c’est seulement un problème pour les effets individuels des variables. Leur effet conjoint est toujours présent dans le modèle, quelle que soit la manière dont les effets individuels sont combinés.
Si nous nous intéressons aux estimations des paramètres individuels, par exemple pour voir comment (et dans quelle mesure) les variables prédictives expliquent le comportement de
- essayer d’obtenir plus de données, afin de réduire les effets de colinéarité apparaissant dans des échantillons spécifiques ou dus à la petite taille de l’échantillon.
- créer un score composite en combinant d’une manière ou d’une autre les variables présentant une colinéarité.
- supprimer une ou plusieurs des variables colinéaires. Vous devez en revanche faire attention à ne pas vous retrouver avec un modèle mal spécifié.
- utiliser la régression pénalisée si
n’est (presque) pas inversible, cela peut restaurer l’unicité de la solution. Les pénalités introduisent un biais, mais peuvent réduire la variance des estimateurs . Les choix populaires incluent la régression en crête (avec une pénalité de ), lasso (pénalité de ), mais ceux-ci requièrent un ajustement pour l’inférence post-sélection
Quelle que soit la méthode utilisée, il est important de comprendre qu’il peut être très difficile (et parfois impossible) d’isoler l’effet individuel d’une variable explicative fortement corrélée avec d’autres.
Exemple 4.22 (Colinéarité des données de college
) On considère l’analyse des données sur l’inéquité salariale dans un college, en incluant cette fois annees
, le nombre d’années depuis l’obtiention du doctorat. On peut penser que, à moins qu’un(e) professeur(e) ait entamé sa carrière dans une autre institution d’enseignement, le nombre d’années de service sera fortement lié à ces derniers. De fait, la corrélation linéaire entre service
et annees
est 0.91. Cette corrélation n’est pas problématique puisque le
4.7.3 Levier et aberrances
L’effet de levier
Les points à fort effet de levier sont ceux qui présentent des combinaisons inhabituelles de variables explicatives. Une observation influente (
Il est important de faire la distinction entre les observations influentes (qui ont une valeur
Si les observations influentes peuvent être détectées en inspectant le levier de chaque observation, les valeurs aberrantes sont plus difficiles à diagnostiquer. Une valeur aberrante se distingue du reste des observations, soit parce qu’elle a une valeur de réponse habituelle, soit parce qu’elle se situe loin de la surface de régression. En gros, une valeur aberrante est une valeur inhabituelle de
Les observations aberrantes et influentes ne doivent pas être négligées parce qu’elles ne sont pas conformes au modèle, mais doivent faire l’objet d’un examen plus approfondi. Elles peuvent motiver une modélisation plus poussée des caractéristiques non prises en compte. Il est également utile de vérifier les erreurs d’enregistrement dans les données (qui peuvent être écartées sans risque). Dans les très grands échantillons, l’impact d’une seule valeur aberrante est, espérons-le, limité. Les transformations de la réponse peuvent contribuer à réduire le caractère aberrant. Sinon, il est possible d’utiliser d’autres fonctions objectives que le critère des moindres carrés ordinaires (telles que celles employées dans la régression robuste); celles-ci pondèrent les observations extrêmes, au détriment de l’efficacité.
4.8 Postulats du modèle et diagnostics
Cette section passe en revue les postulats du modèle énoncés pour permettre l’inférence statistique à l’aide du modèle linéaire et des différents résidus qui servent d’éléments de base pour les diagnostics graphiques. Nous étudions les conséquences de la violation de ces postulats et décrivons des stratégies d’atténuation potentielles, dont beaucoup sont abordées dans d’autres chapitres.
Jusqu’à présent, nous avons ajusté des modèles et testé la significativité des coefficients sans valider notre modèle. La fiabilité des valeurs
- linéarité et additivité: la moyenne de
est , - homoscédasticité: la variance des observations
est constante, - indépendence des observations (conditionnellement aux covariables),
- normalité.
Lorsque nous effectuons un test d’hypothèse et que nous ne rejetons pas l’hypothèse nulle, c’est soit parce qu’elle est vraie, soit par manque de preuves. Il en va de même pour la vérification de la validité des postulats du modèle: le raisonnement scientifique veut que nous ne puissions pas savoir avec certitude si ces derniers sont vrais. Notre stratégie consiste donc à utiliser les implications des hypothèses du modèle linéaire pour créer des outils de diagnostics graphiques, afin de s’assurer qu’il n’y a pas de violation flagrante des postulats. Toutefois, il est important de se garder de surinterpréter les diagnostics graphiques: l’oeil humain est très doué pour trouver des schémas inexistants.
4.8.1 Postulat de linéarité et d’additivité
Le postulat de linéarité signifie que le modèle moyen est correctement spécifié, que toutes les covariables pertinentes ont été incluses et que leur effet est correctement spécifié (y compris les effets non linéaires et les interactions). L’additivité sous-tend que le modèle peut être exprimé comme la somme de moyenne plus aléa. Pour vérifier que la surface de réponse du modèle linéaire est adéquate, nous dessinons un nuage de points de
S’il existe une structure résiduelle dans les graphiques des résidus ordinaires en fonction (a) des valeurs ajustées ou (b) des variables explicatives, un modèle plus complexe peut être ajusté, y compris un contenant des interactions, des fonctions non linéaires, etc. Si l’effet d’une variable explicative est clairement non linéaire et compliqué, des termes de lissage peuvent être ajoutés (nous ne couvrirons pas les modèles additifs généralisés dans ce cours).
La représentation graphique des résidus en fonction des variables explicatives omises peut également servir à vérifier que tout le pouvoir explicatif de la covariable omise est déjà expliqué par les colonnes de
Si une variable importante a été omise et n’est pas disponible dans l’ensemble de données, l’effet de cette variable est capturé à la fois par les erreurs (la partie orthogonale à la matrice du modèle
4.8.2 Postulat d’homoscédasticité
Si la variance des aléas est la même pour toutes les observations (homoscédasticité), celle des observations
Si les aléas (ou variables réponses) sont hétéroscédastiques (variance non constante), les effets estimés des variables (les paramètres
L’examen de nuages de points des résidus studentisés externes en fonction des régresseurs (ou des valeurs ajustées), appelé diagrammes de niveau et de dispersion, est instructif — par exemple, nous voyons souvent un modèle en entonnoir lorsqu’il y a une augmentation de la variance dans le tracé des résidus studentisés externes en fonction de la valeur ajustée, ou encore dans les boîtes à moustache pour une variable catégorielle comme dans la Figure 4.15. Cependant, si nous voulons ajuster un lissage local pour observer les tendances, il est préférable de tracer la valeur absolue des résidus
Une extension évidente du modèle linéaire consiste à permettre à la variance de varier en fonction des variables explicatives, généralement des covariables catégorielles. Cela est facile à faire en modifiant la vraisemblance et nous couvrirons cette approche plus en détail.
Nous pouvons effectuer des tests d’hypothèse pour l’hypothèse d’homogénéité (égalité) de la variance. Les tests les plus couramment utilisés sont le test de Bartlett, un test du rapport de vraisemblance sous l’hypothèse que les données sont tirées d’une loi normale, avec une correction de Bartlett pour améliorer l’approximation
Exemple 4.23 (Violation du postulat d’homoscédasticité)
Qu’arrive-t-il si la variance est inégale? Considérons un problème simple de comparaison de moyennes entre deux groupes, qui revient à effectuer un
Souvent, une variance inégale se produit parce que le modèle n’est pas additif. Vous pouvez utiliser des transformations stabilisant la variance (par exemple, pour des effets multiplicatifs) afin de garantir une variance à peu près égale dans chaque groupe. Dans ce cas de figure, une transformation logarithmique (ou une transformation de Box–Cox) peut aider à stabiliser la variance, mais il faut que la réponse soit positive. Une autre option consiste à utiliser un modèle adapté au type de réponse que vous avez (y compris les données de décompte et les données binaires). Enfin, il peut être nécessaire de modéliser explicitement la variance dans des modèles plus complexes (y compris les mesures répétées) lorsqu’il y a un effet d’apprentissage au fil du temps et que la variabilité diminue en conséquence. Consultez un expert si nécessaire.
Les économistes utilisent fréquemment des estimateurs sandwich (White 1980), en remplaçant l’estimateur usuel de la matrice de covariance des
4.8.3 Postulat d’indépendance
Habituellement, l’indépendance des observations découle directement du type d’échantillonnage utilisé — cette hypothèse est implicitement vraie si les observations représentent un échantillon aléatoire de la population. Ce n’est généralement pas le cas pour les données longitudinales, qui contiennent des mesures répétées des mêmes individus au fil du temps. De même, les séries temporelles ne sont pas constituées d’observations indépendantes. Si nous voulons inclure tous les points temporels dans l’analyse, nous devons tenir compte de l’éventuelle dépendance (corrélation) entre les observations.
Quel est l’impact de la dépendance entre les mesures? D’un point de vue heuristique, les mesures corrélées contiennent moins d’informations que les mesures indépendantes. Dans le cas le plus extrême, il n’y a pas d’information supplémentaire et les mesures sont identiques, mais le fait de les ajouter plusieurs fois gonfle indûment la la taille de l’échantillon. Si nous ignorons la corrélation, les erreurs-type estimées sont trop petites, car la taille effective de l’échantillon est inférieure au nombre d’observations. Cela enfle la statistique et conduit à des rejets plus fréquents des hypothèses nulles, par erreur.
Le manque d’indépendance peut également avoir des conséquences dramatiques sur l’inférence et conduire à de fausses conclusions: la Figure 4.17 montre un exemple avec des échantillons corrélés au sein d’un groupe (ou, de manière équivalente, des mesures répétées d’individus) avec cinq observations par groupe. L’axe des ordonnées montre la proportion de fois où l’hypothèse nulle est rejetée alors qu’elle ne devrait l’être en principe que 5% du temps. Ici, comme les données sont générées à partir du modèle nul (moyenne égale) avec une variance égale, l’inflation du nombre d’erreurs de type I est alarmante et l’inflation du niveau du test est substantielle même avec une corrélation très limitée entre les mesures.
La première source de dépendance est constituée par les données groupées, c’est-à-dire les mesures prises sur des sujets qui ne sont pas indépendants les uns des autres (famille, groupes, etc.) On distingue entre données longitudinales, qui sont des mesures répétées sont effectuées sur les mêmes sujets (quelques points temporels) et séries chronologiques, dont la longueur et la fréquence d’échantillonage est plus élevée. Les séries temporelles nécessitent des modèles spécifiques qui ne sont pas abordés dans ce cours. En raison de l’autocorrélation, les erreurs positives ont tendance à être suivies d’erreurs positives, etc. Nous pouvons tracer les résidus en fonction du temps et un nuage de points des résidus retardés
Toutefois, les diagrammes de résidus décalés ne montrent que la dépendance au premier décalage entre les observations. Pour les séries temporelles, nous pouvons utiliser un corrélogramme, c’est-à-dire un diagramme à bandes de la corrélation entre deux observations distantes de
Pour
Si la série est corrélée, l’autocorrélation de l’échantillon se situera probablement en dehors des intervalles de confiance ponctuels, comme le montre la Figure 4.19. La présence d’autocorrélation nécessite de modéliser la corrélation entre les observations de manière explicite à l’aide d’outils spécifiques issus de la littérature sur les séries temporelles. Nous examinerons toutefois les modèles
Lorsque les observations sont positivement corrélées, les erreurs-types estimées indiquées par le logiciel sont trop petites. Cela signifie que nous sommes trop confiants et que nous rejetterons l’hypothèse nulle plus souvent que nous ne le devrions si l’hypothèse nulle était vraie (erreur de type I gonflée, plus de faux positifs).
Comment prendre en compte la dépendance entre observations? L’idée principale est de modéliser la corrélation et la variance, en partant du vecteur complet d’observations et en supposant que
4.8.4 Postulat de normalité
Le postulat de normalité des aléas est commode, mais pas strictement nécessaire dans la majorité des cas pour la validité des tests sur les coefficients ou les énoncés reliés à la moyenne prédite. Si les aléas suivent une loi normale, les estimateurs des moindres carrés et du maximum de vraisemblance de
Parfois, des transformations peuvent améliorer la normalité : si les données sont asymétriques à droite et que la variable réponse est strictement positive, un modèle log-linéaire peut être plus adéquat Section 4.8.5. Ceci peut être évalué en regardant le diagramme quantile-quantile des résidus studentisés externes. Si la réponse
Si l’aléa
Les graphiques quantile-quantile sont abordés dans la Définition 1.14, mais leur interprétation nécessite de la pratique. Par exemple, la Figure 4.21 montre de nombreux scénarios courants qui peuvent être détectés à l’aide de diagrammes quantile-quantile. Les données discrètes sont responsables des motifs en escalier, les données asymétriques à droite ont des quantiles bas trop élevés et des quantiles hauts trop bas par rapport aux positions de tracé, les données à ailes lourdes ont des observations élevées de part et d’autre et les données bimodales conduisent à des sauts dans le tracé.
Exemple 4.24 (Diagnostics graphiques pour l’inéquité salariale) On considère le modèle pour les données college
avec toutes les observations (sans interaction) pour voir si les postulats tiennent la route.
Sur la base des graphiques de Figure 4.22, nous constatons qu’il existe une hétéroscédasticité au sein des échelons. Comme le nombre d’années pour les adjoint(e)s est limité et que tous les professeurs assistants ont été embauchés au cours des six dernières années, il y a moins de disparité dans leurs revenus. Il est important de ne pas confondre le modèle sur l’axe
On effectue quelques tests avec les résidus studentisés externes pour les données college
pour valider ce que les diagnostics graphiques indiquent.
<- rstudent(modlin1_college)
r # Test F de Levene
::leveneTest(r ~ echelon, center = "mean", data = college)
car#> Levene's Test for Homogeneity of Variance (center = "mean")
#> Df F value Pr(>F)
#> group 2 50 <2e-16 ***
#> 394
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test du score avec Breusch-Pagan
::ncvTest(modlin1_college, var.formula = ~ echelon)
car#> Non-constant Variance Score Test
#> Variance formula: ~ echelon
#> Chisquare = 70.2, Df = 2, p = 6e-16
Pour les données de collège, on spécifie donc plutôt
library(nlme) # modèles mixtes et structures de corrélation
<- nlme::gls(
modlin.college2 model = salaire ~ echelon + domaine + sexe + service, # spécification de la moyenne
weights = nlme::varIdent(form = ~1 | echelon), # variance constante par échelon
data = college)
plot(modlin.college2)
::Anova(modlin.college2)
car#> Analysis of Deviance Table (Type II tests)
#>
#> Response: salaire
#> Df Chisq Pr(>Chisq)
#> echelon 2 363.50 <2e-16 ***
#> domaine 1 91.04 <2e-16 ***
#> sexe 1 1.80 0.180
#> service 1 2.97 0.085 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le modèle est ajusté par maximum de vraisemblance restreint avec la fonction gls
du paquet nlme
. La modification semble suffisante pour capturer l’hétéroscédasticité dans le diagramme des résidus standardisés vs valeurs ajustées.
On pourrait aussi essayer d’utiliser la matrice sandwich, en remplaçant l’estimateur de
# Matrice de covariance sandwich
<- car::hccm(modlin1_college)
vcov_HCE # Statistiques de Wald
::coeftest(modlin1_college, vcov. = vcov_HCE) lmtest
4.8.5 Transformation de la variable réponse
Si la réponse est strictement positive, certaines options peuvent atténuer le manque d’additivité, plus particulièrement les relations multiplicatives entre la moyenne et la variance. Si les données sont asymétriques et que la réponse est strictement positive, un modèle log-linéaire peut être plus approprié et les paramètres peuvent être interprétés.
Écrivons le modèle log-linéaire
Comparez le rapport
Parfois, on veut considérer une transformation à la fois de la réponse et d’une variable explicative positive continue, un modèle log-log. Considérons le cas où on prend le logarithme de
Exemple 4.25 (Modèle de Cobb-Douglas) Considérons par exemple la fonction de production Cobb–Douglas (Douglas 1976), qui spécifie que la production économique
Proposition 4.6 (Transformation de Box–Cox) Parfois, le postulat de normalité de l’erreur dans une régression linéaire ne tient pas. Avec des données strictement positives, on peut envisager une transformation de Box–Cox (Box et Cox 1964),
Si on suppose que
La log-vraisemblance profilée est donc
Nous ne pouvons pas comparer les modèles ajustés à
La transformation de Box–Cox n’est pas une solution miracle et doit être réservée aux cas où la transformation réduit l’hétéroscédasticité (variance inégale) ou crée une relation linéaire entre les explications et la réponse. La théorie fournit une explication convaincante des données avec, par exemple, la fonction de production Cobb-Douglas utilisée en économie qui peut être linéarisée par une transformation logarithmique. Plutôt que de choisir une transformation ad hoc, on pourrait choisir une transformation logarithmique si la valeur 0$ est incluse dans l’intervalle de confiance à 95%, car cela améliore l’interprétabilité.
Exemple 4.26 (Transformation de Box–Cox pour données sur les poisons) Box et Cox (1964) modélisent le temps de survie de 48 animaux sur la base d’un essai aléatoire. Les données sur les poisons
sont équilibrées, 3 poisons ayant été administrés avec 4 traitements à 4 animaux chacun. Nous pourrions envisager une ANOVA à deux facteurs sans interaction, étant donné le peu d’observations pour chaque combinaison. Le modèle s’écrit alors
Le tracé des valeurs ajustées par rapport aux résidus montre que le modèle n’est pas additif (panneau du milieu de la Figure 4.23); il y a également des indications que la variance augmente avec la réponse moyenne. Le modèle est inadéquat: les temps de survie les plus faibles sont sous-estimés, ce qui signifie que les résidus sont positifs, de même que les réponses moyennes. Un test formel de non-additivité indique également la non-additivité (Davison 2003, Exemple 8.24). Dans l’ensemble, l’ajustement du modèle est médiocre et toute conclusion tirée de celui-ci est douteuse.
On pourrait envisager d’utiliser un Box–Cox pour trouver une transformation appropriée des résidus afin d’améliorer la normalité. Une analyse des résidus dans les quatre premiers graphiques de Figure 4.23 montre des signes d’hétéroscédasticité en fonction du poison et du traitement. Ceci est évident en regardant le graphique des résidus ordinaires, qui montre une augmentation de la variance avec le temps de survie. Le tracé quantile-quantile dans le tracé du milieu à droite montre quelques signes d’écart par rapport à la normalité, mais la non-linéarité et l’hétéroscédasticité le masquent.
L’intervalle de confiance à 95% basé sur la log-vraisemblance profilée pour le paramètre du modèle Box–Cox contient
4.9 Remarque finale
La régression linéaire est le modèle statistique le plus connu et le plus utilisé. Le nom peut sembler réducteur, mais de nombreux tests statistiques (t-tests, ANOVA, Wilcoxon, Kruskal–Wallis) [peuvent être formulés à l’aide d’une régression linéaire] (https://lindeloev.github.io/tests-as-linear/linear_tests_cheat_sheet.pdf), tandis que [des modèles aussi divers que les arbres, les composantes principales et les réseaux neuronaux profonds ne sont que des modèles de régression linéaire déguisés] (https://threadreaderapp.com/thread/1286420597505892352.html). Ce qui change sous le capot d’un modèle à l’autre, c’est la méthode d’optimisation (par exemple, les moindres carrés ordinaires, l’optimisation sous contrainte ou la descente de gradient stochastique) et le choix des variables explicatives entrant dans le modèle (bases de splines pour la régression non paramétrique, variable indicatrice sélectionnée via une recherche glouton pour les arbres, fonctions d’activation pour les réseaux neuronaux).