03. Modèles linéaires
2024
Un modèle stochastique (ou aléatoire) combine typiquement
Les modèles sont des “golems” qui servent à obtenir des réponses à nos questions.
Un modèle linéaire est un modèle qui décrit la moyenne d’une variable réponse continue \(Y_i\) d’un échantillon aléatoire de taille \(n\) comme fonction linéaire des variables explicatives (également appelés prédicteurs, régresseurs ou covariables) \(X_1, \ldots, X_p\), \[\begin{align} \underset{\text{moyenne conditionnelle}}{\mathsf{E}(Y_i \mid \boldsymbol{X}_i=\boldsymbol{x}_i)}=\mu_i=\underset{\substack{\text{combinaison linéaire (somme pondérée)}\\ \text{de variables explicatives}}}{\beta_0 + \beta_1x_{i1} + \cdots + \beta_p x_{ip}}\equiv \mathbf{x}_i\boldsymbol{\beta}. \end{align}\] où
Pour l’observation \(i\), on peut écrire \[\begin{align*} \underset{\text{observation}\vphantom{\mu_i}}{Y_i} = \underset{\text{moyenne } \mu_i}{\vphantom{Y_i}\mathbf{x}_i\boldsymbol{\beta}} + \underset{\text{aléa}\vphantom{\mu_i}}{\vphantom{Y_i}\varepsilon_i}, \end{align*}\] où \(\varepsilon_i\) sont des aléas indépendents additifs satisfaisant:
Le modèle linéaire normal \[Y_i \mid \boldsymbol{X}_i=\boldsymbol{x}_i \sim \mathsf{normale}(\mathbf{x}_i\boldsymbol{\beta}, \sigma^2).\]
On peut toujours transformer une variable continue en variable catégorielle pour un modèle pour \(Y\) en fonction de \(X\), mais au prix de coefficients supplémentaires.
Pour simplifier la notation, nous agrégeons les observations en utilisant la notation vectorielle et matricielle suivante: \[\begin{align*} \boldsymbol{Y} = \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{pmatrix}, \; % \boldsymbol{\varepsilon} = %\begin{pmatrix} % \varepsilon_1 \\ % \varepsilon_2 \\ % \vdots \\ % \varepsilon_n % \end{pmatrix}, \; \mathbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix} , \; \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{pmatrix} \end{align*}\] On désigne par \(\mathbf{X}\) la matrice du modèle \(n \times (p+1)\) constituée d’une colonne de uns, concaténée aux \(p\) vecteurs colonnes de variables explicatives.
La \(i\)e ligne de \(\mathbf{X}\) est dénotée \(\mathbf{x}_i\).
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 participants ont été recrutés à l’aide d’un panel en ligne.
Variables:
prodeval
: score moyen d’évaluation du produit sur trois échelles de 9 points (les valeurs les plus élevées sont les meilleures)familiarity
: échelle de Likert de 1 à 7 pour la familiarité avec la marqueconsistency
: groupes d’images et de textes, soit « cohérents », soit « incohérents ».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.
Variables:
group
: facteur pour le groupe expérimental, soit DRTA
, TA
et DR
.pretest1
: score (sur 16) sur le test pré-expérience pour la tâche de détection des erreursposttest1
: score (sur 16) sur le test post-expérience pour la tâche de détection des erreursOn 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. Les données observationnelles college
du paquet hecmodstat
incluent les variables suivantes
salaire
: salaire de professeurs pendant l’année académique 2008–2009 (en 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
ou femme
.service
: nombre d’années de service.annees
: nombre d’années depuis l’obtention du doctorat.L’étude 1 de Moon et VanEpps (2023) (données MV23_S1
, paquet hecedsm
) porte sur la proportion de donateurs à un organisme de charité. 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.
Variables:
before
: est-ce que la personne a déjà donné à l’organisme de charité (0 pour non, 1 pour oui)condition
: facteur pour le groupe expérimental, soit open-ended
pour un montant libre ou quantity
pour le groupe avec montants suggérésamount
: montant du don, NA
si la personne n’a rien donné.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
Voir le Chapitre 11 de Alexander (2023) pour des exemples.
Vérifier
factor
).Considérons un modèle linéaire pour la note moyenne d’évaluation du produit, prodeval
, en fonction de la familiarité de la marque et du facteur expérimental consistency
.
data(LC19_S1, package = "hecedsm")
str(LC19_S1)
## tibble [96 × 5] (S3: tbl_df/tbl/data.frame)
## $ prodeval : num [1:96] 9 8.33 8.67 7.33 9 ...
## $ familiarity: int [1:96] 7 7 7 7 6 5 7 7 4 7 ...
## $ consistency: Factor w/ 2 levels "consistent","inconsistent": 1 1 1 1 1 1 1 1 1 1 ...
## $ gender : Factor w/ 2 levels "male","female": 1 2 1 2 1 1 2 1 1 1 ...
## $ age : int [1:96] 22 26 35 26 39 34 30 33 24 42 ...
length(unique(LC19_S1$prodeval))
## [1] 19
La variable réponse prodeval
est fortement discrétisée, avec seulement 19 valeurs uniques comprises entre 2.33 et 9.
La variable consistency
vaut 0
si la description du texte est cohérente avec l’image, et 1
si elle est incohérente.
Le salaire augmente avec les années de service, mais il y a une plus grande hétérogénéité dans le salaire des professeurs titulaires.
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.
Les variables annees
et service
sont fortement corrélées, avec une corrélation linéaire de 0.91.
Il y a peu de femmes dans l’échantillon. 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 |
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, les 235 observations manquantes seront supprimées.
NA
en zéros.Les variables binaires donate
et before
sont toutes deux des facteurs encodés comme 0
/1
.
Avec des données expérimentales, seules les variables manipulées expérimentalement (affectation aléatoire aux groupes) sont nécessaires.
Pour les données observationnelles, nous avons besoin d’un modèle pour prendre en compte les facteurs de confusion potentiels.
La spécification de la moyenne est \[\begin{align*} \mathsf{E}(Y_i \mid \boldsymbol{X}_i = \boldsymbol{x}_i) = \beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}. \end{align*}\] L’ordonnée à l’origine \(\beta_0\) est la valeur moyenne de \(Y\) lorsque toutes les variables explicatives du modèles sont nulles, soit \(\boldsymbol{x}_i=\boldsymbol{0}_p\). \[\begin{align*} \beta_0 &= \mathsf{E}(Y \mid X_1=0,X_2=0,\ldots,X_p=0) \\ &= \beta_0 + \beta_1 \times 0 + \beta_2 \times 0 + \cdots + \beta_p \times 0 \end{align*}\] Bien sur, il se peut que cette interprétation n’ait aucun sens dans le contexte étudié. Centrer les variables explicatives numériques (pour que leurs moyennes soit zéro) permet de rendre l’ordonnée à l’origine plus interprétable.
En régression linéaire, le paramètre \(\beta_j\) mesure l’effet de la variable \(X_j\) sur la variable \(Y\) une fois que l’on tient compte des effets des autres variables explicatives.
\[\begin{align*} \beta_1 &= \mathsf{E}(Y \mid X_1=x_1+1,X_2=x_2, \ldots, X_p=x_p) \\ & \qquad \qquad - \mathsf{E}(Y \mid X_1=x_1,X_2=x_2,\ldots,X_p=x_p) \\ &= \left\{\beta_1 (x_1+1) + \beta_2 x_2 + \cdots \beta_p x_p \right\} \\ & \qquad \qquad -\left\{\beta_1 x_1 + \beta_2 x_2 + \cdots \beta_p x_p \right\} \end{align*}\]
On définit l’effet marginal comme la dérivée première de la moyenne conditionnelle par rapport à \(X_j\), soit \[\text{effet marginal de }X_j = \frac{\partial \mathsf{E}(Y \mid \boldsymbol{X})}{ \partial X_j}.\]
Le coefficient \(\beta_j\) est aussi l’effet marginal de la variable \(X_j\).
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 \[\begin{align*}
\mathsf{E}(\texttt{amount} \mid \texttt{condition})&= \beta_0 + \beta_1 \mathbf{1}_{\texttt{condition}=\texttt{quantity}}.
\\&= \begin{cases}
\beta_0, & \texttt{condition}=0, \\
\beta_0 + \beta_1 & \texttt{condition}=1.
\end{cases}
\end{align*}\]
open-ended
).quantity
) est \(\beta_0 + \beta_1 = \mu_1\) etopen-ended
et le groupe quantity
.Même si le modèle linéaire définit une droite, cette dernière ne peut être évaluée qu’à \(0\) ou \(1\)
Les données sont fortement discrétisées, avec beaucoup de doublons et de zéros, mais la taille de l’échantillon \((n=869)\) est conséquente.
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
.
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
BSJ92$group[c(1,23,47)]
## [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 \[\mathsf{E}(Y \mid \texttt{group})= \beta_0 + \beta_1\mathbf{1}_{\texttt{group}=\texttt{DRTA}} + \beta_2\mathbf{1}_{\texttt{group}=\texttt{TA}}.\] Puisque la variable group
est catégorielle avec \(K=3\) niveaux, il nous faut mettre \(K-1 = 2\) variables indicatrices.
Avec la paramétrisation en termes de traitements (option par défaut), on obtient
group=DRTA
et zéro sinon,group=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, \(\mu_{\texttt{DR}}=\beta_0\).
(Intercept) | groupDRTA | groupTA | |
---|---|---|---|
DR | 1 | 0 | 0 |
DRTA | 1 | 1 | 0 |
TA | 1 | 0 | 1 |
Si group
=DR
(référence), les deux variables indicatrices binaires groupDRTA
et groupTA
sont nulles. La moyenne de chaque groupe est
Ainsi, \(\beta_1\) est la différence de moyenne entre les groupes DRTA
etDR
, et de la même façon \(\beta_2=\mu_{\texttt{TA}}- \mu_{\texttt{DR}}\).
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
\[\begin{align*} \texttt{salaire} &= \beta_0 + \beta_1 \mathbf{1}_{\texttt{sexe}=\texttt{femme}} +\beta_2 \mathbf{1}_{\texttt{domaine}=\texttt{theorique}} \\&\quad +\beta_3 \mathbf{1}_{\texttt{echelon}=\texttt{aggrege}} +\beta_4 \mathbf{1}_{\texttt{echelon}=\texttt{titulaire}} +\beta_5 \texttt{service} + \varepsilon. \end{align*}\]
\(\widehat{\beta}_0\) | \(\widehat{\beta}_1\) | \(\widehat{\beta}_2\) | \(\widehat{\beta}_3\) | \(\widehat{\beta}_4\) | \(\widehat{\beta}_5\) |
---|---|---|---|---|---|
86596 | -4771 | -13473 | 14560 | 49160 | -89 |
Considérons une matrice de modèle \(\mathbf{X}\) et une équation pour la moyenne du modèle linéaire de la forme \(\mathsf{E}(Y_i) = \mathbf{x}_i\boldsymbol{\beta}\).
Le modèle linéaire comprend
Nous voulons trouver le vecteur de paramètres \(\boldsymbol{\beta} \in \mathbb{R}^{p+1}\) qui minimise l’erreur quadratique moyenne, c’est-à-dire la distance verticale entre les valeurs ajustées \(\widehat{y}_i=\mathbf{x}_i\widehat{\boldsymbol{\beta}}\) et les observations \(y_i\).
Le problème d’optimisation est \[\begin{align*} \widehat{\boldsymbol{\beta}}&=\mathrm{arg min}_{\boldsymbol{\beta} \in \mathbb{R}^{p+1}}\sum_{i=1}^n (y_i-\mathbf{x}_i\boldsymbol{\beta})^2 \\&=(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})^\top(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta}). \end{align*}\]
Si la matrice \(n \times p\) \(\mathbf{X}\) est de plein rang, c’est-à-dire que ses colonnes ne sont pas des combinaisons linéaires les unes des autres, la forme quadratique \(\mathbf{X}^\top \mathbf{X}\) est inversible et nous obtenons \[ \widehat{\boldsymbol{\beta}} = \left(\mathbf{X}^\top \mathbf{X}\right)^{-1}\mathbf{X}^\top \boldsymbol{y}. \qquad(1)\] Cette formule est celle de l’estimateur des moindres carrés ordinaires (MCO). Puisqu’il existe une solution analytique, aucune optimisation numérique n’est requise.
Si \(Y_i \sim \mathsf{normale}(\mathbf{x}_i\boldsymbol{\beta}, \sigma^2)\) et les données sont indépendantes, la log-vraisemblance du modèle linéaire normal s’écrit \[\begin{align*} \ell(\boldsymbol{\beta}, \sigma)&\propto-\frac{n}{2} \ln (\sigma^2) -\frac{(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})^\top(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})}{2\sigma^2}. \end{align*}\]
Maximiser la log-vraisemblance par rapport à \(\boldsymbol{\beta}\) équivaut à minimiser la somme du carré des erreurs \(\sum_{i=1}^n (y_i - \mathbf{x}_i\boldsymbol{\beta})^2\), indépendamment de la valeur de \(\sigma\).
On recouvre l’estimateur des moindres carrés ordinaires \(\widehat{\boldsymbol{\beta}}\).
Pour obtenir les EMV de \(\sigma^2\), on se sert de la log vraisemblance profilée. En excluant les termes constants qui ne dépendent pas de \(\sigma^2\), \[\begin{align*} \ell_{\mathrm{p}}(\sigma^2) &\propto-\frac{1}{2}\left\{n\ln\sigma^2+\frac{1}{\sigma^2}(\boldsymbol{y}-\mathbf{X}\hat{\boldsymbol{\beta}})^\top(\boldsymbol{y}-\mathbf{X}\hat{\boldsymbol{\beta}})\right\}. \end{align*}\]
En différenciant chaque terme par rapport à \(\sigma^2\) et en fixant le gradient à zéro, on obtient \[\begin{align*} \frac{\partial \ell_{\mathrm{p}}(\sigma^2)}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + \frac{(\boldsymbol{y}-\mathbf{X}\hat{\boldsymbol{\beta}})^\top(\boldsymbol{y}-\mathbf{X}\hat{\boldsymbol{\beta}})}{2\sigma^4} = 0 \end{align*}\]
On déduit que l’estimateur du maximum de vraisemblance est la moyenne des carrés des résidus, \[\begin{align*} \widehat{\sigma}^2&=\frac{1}{n}(\boldsymbol{Y}-\mathbf{X}\hat{\boldsymbol{\beta}})^\top(\boldsymbol{Y}-\mathbf{X}\hat{\boldsymbol{\beta}})\\&= \frac{1}{n} \sum_{i=1}^n (y_i - \mathbf{x}_i\widehat{\boldsymbol{\beta}})^2= \frac{\mathsf{SC}_e}{n}; \end{align*}\]
L’estimateur sans biais habituel de \(\sigma^2\) calculé par le logiciel est \[S^2=\mathsf{SC}_e/(n-p-1),\] où le dénominateur est la taille de l’échantillon moins le nombre de paramètres de la moyenne.
Les entrées de la matrice d’information observée sont \[\begin{align*} -\frac{\partial^2 \ell(\boldsymbol{\beta}, \sigma^2)}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^\top} &= \frac{1}{\sigma^2} \frac{\partial \mathbf{X}^\top(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top} = \frac{\mathbf{X}^\top\mathbf{X}}{\sigma^2}\\ -\frac{\partial^2 \ell(\boldsymbol{\beta}, \sigma^2)}{\partial \boldsymbol{\beta}\partial \sigma^2} &=- \frac{\mathbf{X}^\top(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})}{\sigma^4}\\ -\frac{\partial^2 \ell(\boldsymbol{\beta}, \sigma^2)}{\partial (\sigma^2)^2} &= -\frac{n}{2\sigma^4} + \frac{(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})^\top(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta})}{\sigma^6}. \end{align*}\]
Si on évalue l’information observée aux EMV, \[\begin{align*} j(\widehat{\boldsymbol{\beta}}, \widehat{\sigma}^2) = \begin{pmatrix} \frac{\mathbf{X}^\top\mathbf{X}}{\widehat{\sigma^2}} & \boldsymbol{0}_{p+1} \\ \boldsymbol{0}_{p+1}^\top & \frac{n}{2\widehat{\sigma^4}} \end{pmatrix} \end{align*}\] car \(\widehat{\sigma}^2=\mathsf{SC}_e/n\) et les résidus sont orthogonaux à la matrice du modèle.
Puisque \(\mathsf{E}(Y \mid \mathbf{X})=\mathbf{X}\boldsymbol{\beta}\), l’information de Fisher est \[\begin{align*} i(\boldsymbol{\beta}, \sigma^2) = \begin{pmatrix} \frac{\mathbf{X}^\top\mathbf{X}}{\sigma^2} & \boldsymbol{0}_{p+1} \\ \boldsymbol{0}_{p+1}^\top & \frac{n}{2\sigma^4} \end{pmatrix} \end{align*}\]
Puisque la loi asymptotique de l’estimateur est normale, les EMV de \(\sigma^2\) et \(\boldsymbol{\beta}\) sont asymptotiquement indépendants car leur corrélation asymptotique est nulle.
Pourvu que la matrice carrée \((p+1)\), \(\mathbf{X}^\top\mathbf{X}\) soit inversible, la variance asymptotique des estimateurs est