La différence entre glm et lm dans R

L’analyse statistique moderne et la modélisation prédictive au sein de l’environnement de programmation statistique R reposent sur un ensemble de paradigmes algorithmiques et mathématiques rigoureusement structurés. Au cœur de cet édifice se trouvent deux fonctions fondamentales, omniprésentes tant dans les travaux de recherche académique que dans les applications industrielles de pointe : lm() et glm(). Si la première de ces fonctions est universellement reconnue comme le point d’ancrage de la modélisation linéaire classique par la méthode des moindres carrés, la seconde constitue une généralisation théorique et computationnelle d’une portée considérablement plus vaste, formalisée par l’architecture des modèles linéaires généralisés (GLM). Comprendre la frontière conceptuelle, computationnelle et méthodologique qui sépare ces deux outils ne constitue pas un simple exercice d’érudition informatique, mais une nécessité absolue pour tout chercheur ou praticien soucieux de la validité de ses inférences statistiques.

Dans la pratique contemporaine de la modélisation, une confusion persiste fréquemment quant aux contextes d’application respectifs de ces deux fonctions. Il arrive régulièrement que lm() soit employée de manière inappropriée sur des données discrètes ou bornées, conduisant à des estimations biaisées, à des résidus hétéroscédastiques et à des prédictions dépourvues de sens logique ou probabiliste. Réciproquement, la fonction glm() est parfois perçue à tort comme un algorithme fondamentalement distinct et indépendant, alors même que le modèle linéaire canonique s’inscrit mathématiquement comme un cas particulier de la famille exponentielle de distributions sous une fonction de lien identité. L’objectif du présent traité est de déconstruire avec une rigueur exhaustive les fondements, les mécanismes sous-jacents, les divergences d’implémentation et les répercussions empiriques de l’arbitrage entre ces deux piliers du calcul statistique.

À travers une exploration systématique des dimensions épistémologiques, syntaxiques, algorithmiques et diagnostiques, cet article offre un guide définitif pour naviguer entre lm() et glm(). En s’appuyant sur des illustrations appliquées à la psychométrie et aux sciences comportementales — disciplines particulièrement vulnérables aux erreurs de spécification distributionnelle —, nous mettrons en lumière les mécanismes internes du moteur d’estimation de R, depuis la décomposition QR directe jusqu’aux itérations du score de Fisher. Cette analyse fournira aux méthodologistes et aux analystes de données les clés conceptuelles requises pour maximiser la robustesse, la précision et la reproductibilité de leurs investigations quantitatives.

1. Fondements épistémologiques et théoriques de la modélisation statistique dans R

1.1 La régression linéaire classique et l’approche des moindres carrés ordinaires

La régression linéaire classique s’enracine dans la théorie de l’estimation linéaire sans biais optimale formalisée par le théorème de Gauss-Markov. Dans le cadre défini par la fonction lm() de R, le modèle postule que la variable dépendante continue observée est générée par une combinaison linéaire de prédicteurs fixes ou aléatoires, à laquelle s’additionne une perturbation aléatoire non observable. Les postulats gaussiens fondamentaux requièrent impérativement la linéarité des paramètres, l’indépendance statistique des observations, l’homoscédasticité — c’est-à-dire une variance d’erreur conditionnelle constante à travers l’ensemble des niveaux des régresseurs — et, pour l’inférence exacte en échantillon fini, la normalité multivariée des erreurs stochastiques.

D’un point de vue géométrique, l’approche des moindres carrés ordinaires (MCO ou OLS pour Ordinary Least Squares) consiste à projeter orthogonalement le vecteur des observations multidimensionnelles sur le sous-espace vectoriel linéaire engendré par les colonnes de la matrice de plan d’expérience (design matrix). L’objectif explicite de cette minimisation est d’annuler la dérivée première de la somme des carrés des erreurs résiduelles par rapport au vecteur des coefficients, fournissant une solution algébrique directe et fermée. Cette méthode ne nécessite pas initialement la spécification intégrale d’une distribution de probabilité pour identifier le point de minimisation géométrique, bien que l’hypothèse d’une loi gaussienne sous-jacente devienne incontournable dès lors que l’on procède à des tests d’hypothèses basés sur les lois de Student et de Fisher.

L’évolution historique de cette implémentation algorithmique s’étend sur plusieurs décennies. Développée originellement au sein des laboratoires Bell par John Chambers et ses collaborateurs dans le langage S, l’implémentation de la régression linéaire a été pensée pour découpler la formulation symbolique des modèles scientifiques de leur résolution numérique. Lors de la transition vers le langage R au début des années 1990 par Ross Ihaka et Robert Gentleman, la routine lm() a été refondue pour intégrer les standards les plus exigeants de l’analyse numérique matricielle, s’appuyant sur des sous-programmes de la bibliothèque Fortran LINPACK puis LAPACK, assurant une fiabilité absolue du calcul des projections orthogonales.

1.2 L’émergence des modèles linéaires généralisés selon le cadre de Nelder et Wedderburn

En 1972, une rupture paradigmatique majeure s’est opérée dans la littérature statistique avec la publication de l’article séminal de John Nelder et Robert Wedderburn au sein du Journal of the Royal Statistical Society. Les auteurs ont démontré que la régression linéaire classique, la régression logistique, le modèle probit et les modèles log-linéaires pour tableaux de contingence ne constituaient pas des familles théoriques disjointes, mais relevaient tous d’une classe analytique unifiée : les modèles linéaires généralisés (GLM). Cette unification repose sur l’affranchissement de l’hypothèse restrictive de normalité des résidus par l’adoption de la famille exponentielle de distributions, englobant notamment les lois gaussienne, binomiale, poissonienne, gamma et gaussienne inverse.

L’architecture structurelle d’un modèle linéaire généralisé repose sur trois composantes interdépendantes rigoureusement articulées. Premièrement, la composante aléatoire spécifie la distribution de probabilité de la variable dépendante, appartenant nécessairement à la famille exponentielle univariée ou de dispersion. Deuxièmement, le prédicteur linéaire, conventionnellement noté eta, combine de manière strictement additive les variables explicatives pondérées par leurs coefficients respectifs. Troisièmement, la fonction de lien mathématique établit une transformation monotone différentiable reliant l’espérance conditionnelle de la réponse au prédicteur linéaire, permettant d’assurer que les prédictions générées respectent scrupuleusement le domaine de définition théorique de la variable réponse.

Cette formalisation a révolutionné la théorie statistique en permettant de traiter des phénomènes dont la variance est intrinsèquement dépendante de l’espérance conditionnelle — une caractéristique observée de façon systématique dans les données discrètes ou fortement asymétriques. Le cadre des GLM transcende ainsi les limitations de la linéarité additive brute en autorisant des relations non linéaires entre les moyennes et les prédicteurs, tout en préservant les propriétés d’optimalité asymptotique de l’estimation par le maximum de vraisemblance.

1.3 Pertinence méthodologique en psychométrie et sciences comportementales

Dans les disciplines vouées à l’évaluation psychologique, aux neurosciences cognitives et aux sciences comportementales, les variables dépendantes se plient rarement aux prérequis idéalisés de la modélisation gaussienne continue. Les chercheurs sont confrontés quotidiennement à des protocoles de mesure employant des échelles d’évaluation de Likert polytomiques, des temps de réaction fortement asymétriques tronqués à zéro, des réussites ou échecs binaires à des items d’aptitude cognitive, ou des fréquences de manifestations symptomatiques. Ces métriques se caractérisent par des distributions intrinsèquement discrètes, bornées et non homoscédastiques, remettant directement en cause l’usage instinctif de la fonction lm().

L’application aveugle de modèles linéaires classiques sur des scores de nature ordinale ou dichotomique génère des risques analytiques substantiels. L’erreur la plus pernicieuse réside dans la violation structurelle de l’homoscédasticité : pour une variable binaire, la variance résiduelle dépend mathématiquement du produit de la probabilité de succès et de son complément, atteignant son maximum pour une probabilité de 0,5 et tendant vers zéro aux extrémités de l’échelle. L’emploi des moindres carrés ordinaires dans ce contexte conduit à une sous-estimation systématique des erreurs-types, générant des taux de faux positifs inacceptables et compromettant irrémédiablement l’inférence statistique sur les effets d’intérêt.

Par conséquent, l’arbitrage méthodologique entre le paradigme linéaire restreint de lm() et le cadre élargi de glm() conditionne directement la validité de construit et la précision inférentielle des recherches comportementales. Recourir à une fonction de lien adéquate permet d’ajuster des modèles psychométriques qui respectent la nature probabiliste des processus latents mesurés, préservant la signification clinique des paramètres estimés face aux distorsions induites par des transformations ad hoc de variables.

2. Analyse syntaxique et structurelle : comparaison formelle des signatures lm() et glm()

2.1 Structure d’appel et décomposition de la formule dans lm()

L’invocation canonique de la fonction lm() dans R est régie par une grammaire formelle héritée de la notation de Wilkinson et Rogers. La signature minimale s’articule autour de la formule symbolique reliant la variable réponse aux variables explicatives par le biais de l’opérateur tilde. Au sein de cette formule, R déploie un moteur d’interprétation syntaxique sophistiqué capable de développer automatiquement les contrastes factoriels, d’inclure implicitement un terme d’interception à l’origine — sauf exclusion explicite via l’opérateur de soustraction — et de générer l’ensemble des interactions d’ordre supérieur à travers les opérateurs de multiplication ou de croisement factoriel.

La gestion native des contrastes par lm() illustre la fluidité du traitement des variables qualitatives. Lorsqu’un facteur est introduit, R examine l’attribut contrasts associé à la variable ou recourt aux options globales de l’environnement, le comportement par défaut étant l’application de contrastes de traitement qui comparent chaque modalité non basale à une modalité de référence. Parallèlement, l’argument optionnel weights permet de pondérer chaque observation par un scalaire strictement positif, implémentant ainsi les moindres carrés pondérés (WLS) pour corriger une hétéroscédasticité connue a priori, tandis que l’argument subset extrait conditionnellement un sous-échantillon sans requérir d’opération de filtrage préliminaire de la matrice de données.

L’objet formule passé à lm() conserve également des méta-informations critiques au sein de son environnement parent, ce qui garantit que l’évaluation des symboles s’effectue dans le contexte lexical approprié. Cette rigueur structurelle permet à lm() d’extraire la matrice d’expérience $X$ en mobilisant en coulisse les fonctions de bas niveau model.frame() et model.matrix(), étapes indispensables pour préparer les matrices numériques destinées aux sous-programmes compilés en langage C et Fortran.

2.2 L’argument discriminateur family dans glm() et sa valeur par défaut

La distinction syntaxique la plus critique entre les deux fonctions réside dans la présence de l’argument family au sein de glm(). Cet argument attend un objet de classe formelle family, encapsulant non seulement le nom de la distribution choisie, mais également un ensemble complet de fonctions vectorisées implémentant la fonction de lien, sa réciproque, la dérivée du lien, la fonction de variance exprimée en fonction de la moyenne, ainsi que la fonction de déviance résiduelle. L’absence d’un tel argument dans lm() traduit son ancrage exclusif dans le paradigme gaussien à lien identité.

Lorsque la fonction glm() est appelée sans que l’utilisateur ne renseigne explicitement l’argument family, R applique par défaut la valeur gaussian(link = "identity"). Cette configuration par défaut engendre une convergence mathématique remarquable : elle aligne les spécifications distributionnelles du GLM sur celles de la régression linéaire classique. Il est toutefois possible d’instancier explicitement des générateurs de familles statistiques diversifiés, tels que binomial(link = "logit"), poisson(link = "log"), ou Gamma(link = "inverse"), en passant soit la chaîne de caractères correspondante, soit l’appel de fonction paramétré.

L’architecture modulaire de l’objet family confère à glm() une extensibilité exceptionnelle. Un utilisateur avancé peut construire une famille entièrement personnalisée en définissant manuellement les méthodes mathématiques requises par l’algorithme d’optimisation, ce qui permet d’ajuster des modèles hors du champ conventionnel, une flexibilité totalement absente de la fonction rigide et monolithique lm().

2.3 Gestion des arguments optionnels et contrôle des itérations d’optimisation

Puisque l’estimation des modèles linéaires généralisés requiert généralement la résolution numérique d’équations non linéaires par des méthodes d’approximation successives, glm() introduit des arguments spécifiques d’une importance capitale pour la stabilité de l’algorithme. L’argument control accepte une liste de paramètres régulateurs générés par la fonction auxiliaire glm.control(). Celle-ci permet de moduler la tolérance de convergence numérique — mesurée par le changement relatif de déviance entre deux itérations —, de fixer le seuil maximal d’itérations autorisées avant interruption, et de solliciter l’affichage détaillé de l’historique d’optimisation.

Par ailleurs, pour pallier les écueils de divergence numérique fréquents dans les espaces de paramètres complexes ou mal conditionnés, glm() autorise l’initialisation déterministe de la recherche numérique grâce aux arguments start, etastart et mustart. Ces vecteurs permettent à l’analyste d’injecter des valeurs a priori plausibles pour les coefficients du modèle, pour le prédicteur linéaire ou directement pour l’espérance mathématique, évitant ainsi que les algorithmes ne dérivent vers des régions indéfinies de la surface de log-vraisemblance où les gradients s’annulent ou deviennent infinis.

Enfin, la modélisation des processus de taux ou de comptage temporellement indexés requiert la prise en compte d’un décalage d’exposition connu, généralement désigné sous le terme d’offset. L’argument formel offset de glm() permet d’intégrer dans le prédicteur linéaire un vecteur numérique dont le coefficient associé est rigidement contraint à la valeur unité ($1.0$). Bien que lm() supporte également la spécification d’un offset, son utilisation dans glm() est conceptuellement distincte, puisqu’elle permet de modéliser des taux d’occurrence par unité de temps, de volume ou de population au sein d’une structure logarithmique, formalisant la décomposition multiplicative essentielle aux processus d’intensité poissonnienne.

3. Mécanismes mathématiques et algorithmes d’estimation sous-jacents

3.1 Moindres carrés ordinaires (OLS) via décomposition QR dans lm()

Le cœur computationnel de la fonction lm() repose sur un algorithme matriciel non itératif garantissant une stabilité numérique exceptionnelle : la décomposition QR de la matrice d’expérience $X$. Au lieu d’inverser directement la matrice de produit scalaire $X^T X$ — une opération extrêmement sensible aux erreurs d’arrondi numérique dès lors que des prédicteurs présentent des degrés de corrélation modérés —, R factorise la matrice $X$, de dimension $n \times p$, sous la forme d’un produit $X = Q R$. Dans cette décomposition, la matrice $Q$ est une matrice orthogonale de dimension $n \times p$ vérifiant $Q^T Q = I_p$, et $R$ est une matrice triangulaire supérieure de dimension $p \times p$, inversible si le rang de $X$ est complet.

En substituant cette factorisation dans l’équation normale des moindres carrés, l’expression formelle du vecteur des estimateurs $\hat{\beta} = (X^T X)^{-1} X^T Y$ se simplifie de manière spectaculaire :
$$R \hat{\beta} = Q^T Y$$
Cette structure triangulaire permet de déterminer les coefficients $\hat{\beta}$ par un algorithme élémentaire de substitution arrière (back-substitution). Ce procédé évite tout recours à une inversion matricielle explicite, réduisant considérablement la propagation des erreurs en arithmétique flottante conforme à la norme IEEE 754.

L’implémentation de bas niveau s’opère via la routine interne dqrls, qui intègre une stratégie de pivotement de colonnes. Ce mécanisme détecte automatiquement la colinéarité exacte : si une colonne de la matrice de plan s’avère linéairement dépendante des précédentes, son élément diagonal dans la matrice $R$ s’annule, conduisant le programme à affecter la valeur NA au coefficient associé et à éliminer la variable redondante du calcul des degrés de liberté, préservant ainsi l’unicité mathématique de la solution retenue pour le sous-espace identifiable.

3.2 Moindres carrés récursivement repondérés (IRLS) et score de Fisher dans glm()

Contrairement au cas linéaire classique, la maximisation de la log-vraisemblance d’un modèle linéaire généralisé ne débouche généralement pas sur une formulation analytique explicite pour le vecteur $\hat{\beta}$. La fonction glm() résout ce problème d’optimisation non linéaire par le biais de l’algorithme des moindres carrés itérativement repondérés (IRLS pour Iteratively Reweighted Least Squares), qui se révèle mathématiquement équivalent à l’algorithme du score de Fisher — une variante de la méthode de Newton-Raphson où la matrice hessienne empirique est remplacée par son espérance mathématique, la matrice d’information de Fisher.

Le déroulement séquentiel de l’IRLS repose sur l’actualisation conjointe, à chaque itération $t$, d’une variable dépendante de travail ajustée, notée $z_t$, et d’une matrice diagonale de poids d’ajustement, notée $W_t$. La variable ajustée est obtenue par linéarisation de premier ordre de la fonction de lien autour de l’estimation courante de l’espérance conditionnelle :
$$z_t = \eta_t + (Y – \mu_t) \left(\frac{d\eta}{d\mu}\right)_t$$
Les poids de travail reflètent quant à eux l’inverse de la variance asymptotique de cette variable linéarisée :
$$w_{ii, t} = \frac{1}{V(\mu_{i, t}) \left(\frac{d\eta_i}{d\mu_i}\right)_t^2}$$
où $V(\mu)$ représente la fonction de variance caractéristique de la famille distributionnelle choisie.

À chaque pas itératif, R calcule une nouvelle estimation des coefficients par la résolution d’un problème de moindres carrés pondérés standard :
$$\hat{\beta}_{t+1} = (X^T W_t X)^{-1} X^T W_t z_t$$
Ce cycle computationnel se répète jusqu’à ce que la déviance résiduelle du modèle stabilise sa variation relative en deçà du seuil de tolérance numérique spécifié, ou jusqu’à ce que le nombre maximal d’itérations soit épuisé. En présence d’une séparation quasi-complète des données dans les régressions logistiques, l’algorithme peut échouer à converger, les coefficients tendant asymptotiquement vers l’infini, ce qui se manifeste dans R par l’émission d’avertissements explicites relatifs à l’occurrence de probabilités ajustées égales à 0 ou 1.

3.3 Complexité algorithmique et consommation des ressources informatiques

Les divergences fondamentales dans la conception algorithmique de lm() et de glm() dictent des profils d’allocation de mémoire vive et d’utilisation du processeur profondément dissemblables. Pour un jeu de données comportant $n$ observations et $p$ prédicteurs, l’exécution de lm() présente une complexité computationnelle théorique dominée par la décomposition QR, requérant un ordre de grandeur de $O(2np^2 – \frac{2}{3}p^3)$ opérations en virgule flottante. Ce calcul est exécuté en une passe unique non itérative, ce qui confère à la fonction une rapidité d’exécution remarquable et une consommation mémoire déterministe, les structures intermédiaires étant allouées de manière contiguë.

À l’inverse, la fonction glm() multiplie cette charge computationnelle par le nombre d’itérations nécessaires à la convergence de l’IRLS, le nombre d’itérations typique variant de 4 à 25 selon la qualité de la distribution et l’éloignement des valeurs initiales. À chaque cycle de recalcul, l’algorithme doit reconstruire la variable dépendante ajustée $z$, recalculer l’ensemble des poids diagonaux $W$, et exécuter à nouveau une décomposition QR pondérée sur la matrice $\sqrt{W}X$. Par conséquent, sur des volumétries de données massives — telles que celles issues des protocoles d’acquisition continue en imagerie cérébrale ou de traces numériques d’apprentissage —, la surcharge computationnelle induite par glm() peut engendrer un facteur multiplicatif de temps de traitement allant de 5 à plus de 50 par rapport à lm().

Face au problème de la colinéarité sévère, les deux fonctions réagissent selon des dynamiques distinctes. Tandis que le pivotement QR de lm() isole instantanément les colonnes collinearisées au rang près, la convergence de l’IRLS dans glm() devient extrêmement lente et instable lorsque la matrice de variance $X^T W X$ approche de la singularité. La surface de log-vraisemblance aplatie conduit à des oscillations de l’estimateur, risquant de déclencher des alertes d’arrêt prématuré avant l’atteinte d’un véritable extremum local ou global.

4. Le cas de convergence parfaite : lm() versus glm(family = gaussian)

4.1 Identité mathématique stricte des coefficients estimés

L’un des théorèmes unificateurs les plus élégants de la théorie statistique stipule que, sous l’hypothèse de perturbations stochastiques identiquement et indépendamment distribuées selon une loi normale d’espérance nulle et de variance constante $\sigma^2$, l’estimateur des moindres carrés ordinaires s’identifie rigoureusement à l’estimateur du maximum de vraisemblance. En d’autres termes, minimiser la somme des carrés des écarts euclidiens revient exactement à maximiser le noyau de la fonction de densité gaussienne multivariée.

Lorsque l’on paramètre la fonction glm() avec l’argument explicite family = gaussian(link = "identity"), la fonction de lien est la transformation identité $g(\mu) = \mu$, de sorte que la dérivée $deta/d\mu$ est universellement égale à l’unité ($1$). Par ailleurs, la fonction de variance pour la distribution gaussienne est constante : $V(\mu) = 1$. En injectant ces simplifications dans les équations de travail de l’IRLS présentées précédemment, on constate immédiatement que la variable ajustée $z_t$ devient identique à la variable réponse originale $Y$, et que la matrice des poids $W_t$ se réduit à une matrice identité scalaire pondérée.

Dès lors, l’algorithme IRLS converge dès la première itération, réalisant exactement le même calcul que l’équation normale résolue par lm(). Il s’ensuit une égalité numérique absolue, jusqu’à la limite de précision d’arrondi de la machine (généralement de l’ordre de $10^{-16}$ en double précision IEEE), entre le vecteur des coefficients $\hat{\beta}$ produit par lm() et celui retourné par glm(..., family = gaussian). Aucun écart substantiel d’estimation ponctuelle ne peut être observé, indépendamment du nombre de régresseurs ou de la taille de l’échantillon.

4.2 Divergences subtiles dans l’estimation de la variance et de la dispersion

En dépit de l’identité absolue constatée au niveau des estimations ponctuelles des coefficients de régression, une divergence méthodologique et computationnelle subtile sépare les sorties de lm() et de glm() quant à l’estimation de la variance résiduelle du modèle, avec des répercussions directes sur les erreurs-types associées aux régresseurs. Cette divergence prend sa source dans le traitement des degrés de liberté de l’échantillon pour corriger le biais de l’estimateur de la variance.

La théorie classique des moindres carrés implémentée dans lm() recourt à un estimateur sans biais de la variance d’erreur résiduelle, noté $s^2$, qui divise la somme des résidus au carré par les degrés de liberté résiduels exacts, soit $n – k$ (où $k$ représente le nombre total de paramètres libres estimés, y compris la constante d’interception) :
$$s^2 = \frac{\sum_{i=1}^n (y_i – \hat{y}_i)^2}{n – k}$$
À l’inverse, l’estimateur pur du maximum de vraisemblance, ancré dans l’asymptotique des GLM, divise conventionnellement la déviance résiduelle par la taille totale de l’échantillon $n$, produisant un estimateur biaisé mais convergent de la variance résiduelle :
$$\hat{\sigma}_{ML}^2 = \frac{\sum_{i=1}^n (y_i – \hat{y}_i)^2}{n}$$

Dans l’environnement R, les concepteurs ont toutefois atténué cette dichotomie au sein de la fonction d’affichage summary.glm() : pour la famille gaussienne, le paramètre de dispersion est ajusté en utilisant la correction par les degrés de liberté résiduels $n – k$. Néanmoins, des divergences d’interfaçage demeurent visibles dans le calcul des métriques de déviance et dans l’extraction de certains critères d’information, l’objet glm encapsulant prioritairement la logique asymptotique du maximum de vraisemblance.

4.3 Implications pratiques pour la réplication empirique en psychologie

Ces nuances formelles ont des répercussions tangibles sur les protocoles de réplication et la transparence méthodologique au sein des revues scientifiques de psychologie quantitative. Lorsqu’un analyste utilise indifféremment lm() ou glm(..., family = gaussian), le lecteur non averti ou l’évaluateur de la revue peut percevoir l’usage de glm() comme l’indice d’une modélisation non standard, supposant implicitement une distribution d’erreur asymétrique ou une fonction de lien non linéaire.

Sur le plan de la reproductibilité logicielle, documenter un modèle gaussien standard via glm() alourdit inutilement le script et induit des dépendances de classes d’objets distinctes lors de l’exportation des résultats vers des outils de mise en forme automatisés tels que Stargazer ou modelsummary. Bien que les coefficients et les tests de Student rapportés par summary() s’avèrent numériquement superposables dans la plupart des configurations contemporaines, le recours à lm() demeure le standard absolu et universellement attendu par la communauté académique pour la régression linéaire ordinaire.

Le choix programmatique doit donc respecter un principe de parcimonie expressive : réserver rigoureusement lm() à l’ensemble des modélisations fondées sur l’hypothèse de résidus gaussiens homoscédastiques, et ne basculer vers glm() que lorsque la violation théorique ou distributionnelle de ces postulats justifie explicitement l’adjonction d’une famille statistique spécifique ou d’une fonction de lien non canonique.

5. La fonction de lien : modélisation non linéaire du prédicteur linéaire

5.1 Rôle théorique et formulation mathématique du lien

La composante architecturale qui confère aux modèles linéaires généralisés toute leur puissance expressive est la fonction de lien, notée $g(\cdot)$. Sur le plan formel, la fonction de lien établit une correspondance fonctionnelle continue, strictement monotone et différentiable entre l’espérance mathématique conditionnelle de la variable réponse, $\mu = E(Y|X)$, et le prédicteur linéaire scalaire $eta$ :
$$\eta = g(\mu) = \sum_{j=0}^p \beta_j X_j$$
La propriété de monotonicité garantit la bijectivité locale de la transformation, ce qui permet de projeter univoquement les prédictions depuis l’espace non contraint du prédicteur linéaire vers l’espace naturel des observations à l’aide de la fonction de lien inverse :
$$\mu = g^{-1}(\eta)$$

Il est fondamental de distinguer conceptuellement la modélisation de la fonction de lien, $g(E(Y|X))$, de la transformation directe de la variable dépendante, $E(g(Y)|X)$, cette dernière étant fréquemment pratiquée de façon heuristique avec la fonction lm(). Dans un modèle linéaire où l’on applique une transformation logarithmique préalable à la variable $Y$, l’équation estimée correspond à :
$$E(\ln(Y)|X) = X\beta$$
Dans un modèle linéaire généralisé avec lien log, la spécification structurelle s’écrit quant à elle :
$$\ln(E(Y|X)) = X\beta iff E(Y|X) = \exp(X\beta)$$
Ces deux formulations ne sont en aucun cas interchangeables : l’espérance d’une fonction non linéaire d’une variable aléatoire ne correspond pas à la fonction de l’espérance de cette variable, une asymétrie directement dictée par l’inégalité de Jensen.

La fonction de lien permet ainsi de préserver la linéarité des paramètres d’intérêt tout en contraignant l’espace de prédiction dans des bornes physiquement ou probabilistiquement cohérentes. Pour une proportion, le lien restreint $\mu$ à l’intervalle ouvert $(0, 1)$ ; pour un comptage, il garantit des valeurs strictement positives dans $(0, +\infty)$.

5.2 Fonctions de lien canoniques versus fonctions de lien non canoniques

Chaque distribution appartenant à la famille exponentielle naturelle dispose d’une fonction de lien privilégiée qualifiée de canonique. Le lien canonique survient de manière organique lorsque le prédicteur linéaire est égalé directement au paramètre naturel de la formulation canonique de la densité de probabilité. Sur le plan statistique, l’adoption d’une fonction de lien canonique confère des propriétés d’optimalité analytique remarquables : elle garantit l’existence de statistiques exhaustives de dimension minimale pour le vecteur des coefficients $\beta$, simplifie la matrice d’information de Fisher en l’identifiant à la matrice hessienne observée, et assure la convergence globale monotone de l’algorithme d’optimisation en éliminant les risques de minima locaux sur la surface de déviance.

Le tableau suivant illustre la cartographie rigoureuse reliant les distributions majeures de la famille exponentielle, leurs paramètres naturels et leurs fonctions de lien canoniques respectives :

• Distribution Gaussienne : Lien canonique identité, $g(\mu) = \mu$, préservant l’additivité brute des effets sans déformation de l’espace métrique.
• Distribution Binomiale / Bernoulli : Lien canonique logit, $g(\mu) = \ln\left(\frac{\mu}{1-\mu}\right)$, modélisant le logarithme du rapport de cotes.
• Distribution de Poisson : Lien canonique logarithmique, $g(\mu) = \ln(\mu)$, contraignant l’espérance à l’ensemble des réels strictement positifs.
• Distribution Gamma : Lien canonique inverse, $g(\mu) = \mu^{-1}$, liant l’échelle des taux à la dégradation temporelle des processus continus.

Le chercheur n’est cependant nullement contraint d’employer systématiquement le lien canonique. Des impératifs théoriques ou substantiels justifient régulièrement le recours à des liens non canoniques. À titre d’exemple, pour la distribution binomiale, la fonction de lien probit — correspondant à l’inverse de la fonction de répartition de la loi normale centrée réduite — s’impose naturellement en psychométrie lorsque l’on postule l’existence d’une propension latente continue régie par un processus gaussien non observé sous-jacent à la réponse catégorielle manifeste. De même, la fonction log-log complémentaire (cloglog), $g(\mu) = \ln(-\ln(1-\mu))$, est mobilisée dès lors que les mécanismes d’apparition d’un comportement présentent une asymétrie prononcée, les taux de transition vers l’événement approchant l’unité beaucoup plus lentement qu’ils ne quittent la borne inférieure.

5.3 Écueils méthodologiques des transformations manuelles de variables dans lm()

Face à des variables fortement asymétriques ou hétéroscédastiques, une pratique historique largement répandue en sciences appliquées consiste à transformer manuellement la variable dépendante à l’aide de fonctions mathématiques usuelles (racine carrée, logarithme népérien, transformation inverse) avant d’appliquer une régression linéaire standard via lm(). Si cette démarche vise à rétablir la normalité et la stabilité de la variance, elle induit des biais méthodologiques majeurs qui sont totalement évités par le recours aux fonctions de lien au sein de glm().

Le premier écueil fondamental réside dans l’inégalité de Jensen, qui établit formellement que pour toute fonction strictement concave comme le logarithme, l’espérance de la transformée est strictement inférieure à la transformée de l’espérance : $E(ln(Y)) < ln(E(Y))$. Par conséquent, lorsque l'analyste applique l'exponentielle aux prédictions issues d'un modèle lm(log(Y) ~ X) pour tenter de retrouver l’échelle métrique originale de sa mesure psychologique, il sous-estime systématiquement la moyenne arithmétique conditionnelle réelle de la population d’intérêt. La correction de ce biais nécessite l’adjonction de facteurs de dispersion correctifs complexes (tels que le facteur d’ajustement de Duan), souvent omis dans les publications empiriques.

Le second écueil méthodologique concerne l’indéfinition mathématique des transformations élémentaires face aux observations nulles ou négatives. Pour contourner l’impossibilité de calculer le logarithme de zéro, une pratique courante consiste à ajouter une constante arbitraire (usuellement notée $ln(Y + c)$ avec $c=1$ ou $c=0.1$). Cette manipulation ad hoc altère fondamentalement la topologie des données, modifie substantiellement l’ampleur des coefficients estimés et rend les inférences hautement dépendantes de la constante arbitrairement choisie. glm() résout nativement cette impasse : la fonction de lien transforme l’espérance mathématique modélisée $\mu$, qui demeure strictement positive pour des distributions adaptées comme la loi de Poisson ou Gamma, autorisant la présence naturelle d’observations individuelles nulles ($y_i = 0$) sans distorsion de l’espace métrique.

6. Modélisation des variables binaires et catégorielles : limites de lm() et primauté de glm()

6.1 Le modèle de probabilité linéaire via lm() et ses défaillances

L’estimation d’une variable réponse dépendante binaire codée de manière dichotomique ($Y in {0, 1}$) par la méthode des moindres carrés ordinaires à l’aide de la fonction lm() est traditionnellement qualifiée de Modèle de Probabilité Linéaire (MPL). Dans ce cadre, la valeur ajustée par le modèle linéaire est directement interprétée comme la probabilité conditionnelle que l’événement survienne, soit $\hat{Y} = P(Y=1|X) = X\beta$. Bien que séduisante par son apparente simplicité d’interprétation des coefficients sous forme de variations directes de probabilités absolues, cette approche accumule des défaillances statistiques et logiques rédhibitoires.

La faille la plus immédiatement visible réside dans la violation intégrale des contraintes axiomatiques élémentaires de Kolmogorov régissant les probabilités. L’opérateur linéaire $X\beta$ n’étant soumis à aucune restriction de bornage supérieur ou inférieur, le modèle de probabilité linéaire est intrinsèquement capable de générer des prédictions absurdes strictement négatives ($P(Y=1) < 0$) ou excédant l’unité ($P(Y=1) > 1$). Ce phénomène survient fréquemment dès lors que les prédicteurs prennent des valeurs extrêmes ou que la prévalence globale du phénomène clinique ou comportemental étudié est fortement déséquilibrée au sein de la population.

Sur le plan inférentiel, le MPL viole de manière irréductible l’hypothèse fondamentale d’homoscédasticité des résidus stipulée par le théorème de Gauss-Markov. Par construction, la variable réponse $Y$ ne pouvant prendre que les valeurs $0$ ou $1$, la distribution du terme d’erreur résiduel $e_i = Y_i – X_i\beta$ est elle-même strictement discrète et conditionnée aux régresseurs, prenant la valeur $1 – X_i\beta$ avec une probabilité $X_i\beta$, et la valeur $-X_i\beta$ avec une probabilité $1 – X_i\beta$. La variance de cette erreur résiduelle conditionnelle s’exprime rigoureusement par :
$$Var(e_i|X_i) = X_i\beta (1 – X_i\beta)$$
Cette variance dépend donc structurellement des valeurs prises par les prédicteurs, induisant une hétéroscédasticité généralisée. De surcroît, la distribution d’échantillonnage de ces résidus ne peut en aucun cas converger vers une loi normale, rendant les tests statistiques de Student et de Fisher délivrés par summary.lm() non fiables en échantillon fini.

6.2 Régression logistique binaire avec glm(family = binomial(link = ‘logit’))

Pour s’affranchir définitivement des impasses logiques et distributionnelles du modèle de probabilité linéaire, la modélisation standard repose sur l’instanciation de la régression logistique binaire par l’appel glm(family = binomial(link = "logit")). Cette formalisation s’appuie sur la loi de probabilité de Bernoulli, où le paramètre fondamental modélisé est la probabilité conditionnelle de succès $\pi = P(Y=1|X)$. L’adoption de la fonction de lien canonique logit transforme cette probabilité en son équivalent logarithmique de rapport de cotes (log-odds) :
$$\ln\left(\frac{\pi}{1 – \pi}\right) = \beta_0 + \beta_1 X_1 + dots + \beta_p X_p$$

L’estimation des paramètres s’opère par la maximisation numérique de la fonction de log-vraisemblance globale de Bernoulli, définie sur l’ensemble de l’échantillon par l’expression :
$$ell(\beta) = \sum_{i=1}^n \left[ y_i \ln(\pi_i) + (1 – y_i) \ln(1 – \pi_i) \right]$$
Cette méthode garantit que les probabilités prédites obtenues via la fonction de lien réciproque, connue sous le nom de fonction sigmoïde ou fonction logistique standard :
$$\pi_i = \frac{1}{1 + \exp(-X_i\beta)}$$
demeurent asymptotiquement confinées dans l’intervalle ouvert strict $]0, 1[$, éradiquant toute possibilité de prédictions probabilistes aberrantes ou impossibles.

Dans le domaine de la psychopathologie clinique — par exemple lors de l’investigation des déterminants associés à la rechute dépressive binaire face à une rémission stable —, la régression logistique offre une métrique d’interprétation substantielle puissante à travers le rapport de cotes (odds ratio ou OR). En appliquant une simple exponentiation au coefficient estimé, $\exp(\hat{\beta}_j)$, le clinicien quantifie précisément la multiplication relative de la cote de survenue de l’événement symptomatique induite par l’accroissement d’une unité de la variable explicative, préservant une validité diagnostique rigoureuse exempte de biais de bornage.

6.3 Extensions polytomiques et alternatives probit

Lorsque la variable de réponse comportementale excède le cadre purement binaire pour embrasser des échelles polytomiques comptant trois modalités ou davantage, les limitations de l’architecture native de glm() commencent à se manifester. Pour les variables ordinales ordonnées — telles que les stades de sévérité d’un trouble anxieux (léger, modéré, sévère) —, l’application de lm() en assignant des scores numériques équidistants ($1, 2, 3$) postule arbitrairement une distance psychologique identique entre chaque échelon, faussant les estimations de pente.

Bien que le moteur standard de glm() permette d’ajuster des modèles binaires avec un lien probit (family = binomial(link = "probit")) — particulièrement adapté à la psychométrie en vertu de son postulat d’une propension latente distribuée selon une courbe normale unitaire de Gauss —, il ne gère pas nativement les réponses polytomiques nominales ou ordinales complexes. Pour ces architectures de données, l’analyste doit transiter vers des packages spécialisés de l’écosystème R étendant le cadre des GLM, tels que le package MASS avec la fonction polr() pour les régressions logistiques à cotes proportionnelles, ou le package nnet avec multinom() pour les choix nominaux non ordonnés.

Ces extensions avancées permettent d’estimer des paramètres de seuils latents (cutpoints) délimitant les zones de transition probabiliste entre les différents états comportementaux observables. Cette étape est essentielle pour modéliser adéquatement la sensibilité et la spécificité des instruments d’évaluation clinique standardisés, soulignant une nouvelle fois la distance méthodologique infranchissable qui sépare l’approche générale des GLM des simplifications excessives de la régression linéaire ordinaire.

7. Traitement des données de comptage et de fréquences comportementales

7.1 Inadéquation des moindres carrés pour les distributions discrètes positives

L’enregistrement quantitatif des fréquences d’apparition de comportements spécifiques au cours d’une période d’observation standardisée engendre des structures de données dénommées données de comptage. Qu’il s’agisse du nombre d’erreurs commises lors d’un test d’attention soutenue, du nombre d’interactions agressives manifestées par un sujet en milieu scolaire, ou de la fréquence d’épisodes de panique hebdomadaires, ces métriques se caractérisent par un support de valeurs strictement restreint à l’ensemble des entiers naturels $\mathbb{N} = {0, 1, 2, dots}$. L’application d’un modèle linéaire standard lm() sur de telles métriques constitue une entorse méthodologique majeure aux postulats statistiques.

Le premier défaut intrinsèque réside dans la troncature naturelle et structurelle des données à zéro. Une fréquence de comportements ne pouvant être négative, la distribution empirique est asymétrique à droite, avec une concentration importante de masses de probabilité sur les faibles valeurs, voire sur le zéro absolu. Dès lors que l’espérance mathématique globale de la variable est modérée, l’hypothèse de symétrie et de normalité des résidus de lm() s’effondre irrémédiablement, induisant des queues de distribution non conformes et générant de nouveau des valeurs d’ajustement prédites totalement aberrantes situées dans les réels négatifs.

Le second écueil, plus structurel encore, est la corrélation systématique entre l’espérance conditionnelle et la variabilité intrinsèque du phénomène. Dans la quasi-totalité des mécanismes comportementaux générateurs de comptages, les sujets ou conditions expérimentales caractérisés par une moyenne d’occurrence élevée présentent également des fluctuations d’une amplitude considérablement plus large que ceux caractérisés par des fréquences rares. Cette dépendance fonctionnelle entre moyenne et variance détruit le postulat d’homoscédasticité, privant les moindres carrés ordinaires de leur efficience et biaisant les inférences relatives aux coefficients de prédiction.

7.2 La régression de Poisson via glm(family = poisson(link = ‘log’))

La formulation théorique rigoureuse pour l’analyse de ces variables d’occurrence repose sur la modélisation du taux sous-jacent d’intensité par une distribution de Poisson via l’instruction glm(family = poisson(link = "log")). La distribution de Poisson assigne à l’observation discrète $Y = y$ la probabilité formelle :
$$P(Y = y|\mu) = \frac{\mu^y \exp(-\mu)}{y!}$$
où $\mu$ représente l’espérance conditionnelle d’occurrence du comportement d’intérêt. La mobilisation de la fonction de lien canonique logarithmique garantit que le prédicteur linéaire préserve la positivité stricte de l’espérance :
$$\ln(\mu) = X\beta iff \mu = \exp(X\beta)$$

L’estimation est exécutée par maximisation de la log-vraisemblance spécifique à la loi de Poisson :
$$ell(\beta) = \sum_{i=1}^n \left[ y_i X_i\beta – \exp(X_i\beta) – \ln(y_i!) \right]$$
Sur le plan de l’interprétation substantielle, cette fonction de lien logarithmique induit une relation purement multiplicative entre les prédicteurs et l’espérance du comptage. L’exponentiation des coefficients de régression, $\exp(\hat{\beta}_j)$, livre une métrique reconnue internationalement en épidémiologie et en psychologie sous l’appellation de Rapport des Taux d’Incidence (IRR pour Incidence Rate Ratio).

Un IRR de $1,25$ associé à un traitement thérapeutique spécifique signifiera rigoureusement que chaque accroissement d’une unité de ce prédicteur élève de 25 % l’incidence moyenne du comportement cible. Cette lecture multiplicative respecte la mécanique comportementale réelle, où les facteurs environnementaux et génétiques amplifient ou réduisent proportionnellement les seuils d’activation du comportement, là où lm() tenterait maladroitement d’ajuster une pente additive uniforme constante sur l’ensemble de l’étendue des mesures.

7.3 Gestion empirique de la surdispersion : quasi-poisson et approches alternatives

L’un des postulats fondamentaux les plus rigides et restrictifs de la distribution de Poisson standard est l’équidispersion stricte : la moyenne conditionnelle de la variable aléatoire doit être mathématiquement égale à sa variance conditionnelle ($Var(Y|X) = E(Y|X) = \mu$). Cependant, dans la pratique empirique des sciences comportementales, cette égalité est très fréquemment invalidée par le phénomène de surdispersion (overdispersion), dans lequel la variance observée excède substantiellement l’espérance conditionnelle prédite. Cette situation est souvent causée par une hétérogénéité individuelle non observée, une agrégation spatio-temporelle ou une corrélation intra-sujet non modélisée.

Le diagnostic primaire de la surdispersion s’obtient au sein de R en examinant le ratio entre la déviance résiduelle affichée dans le sommaire du modèle et les degrés de liberté résiduels :
$$\hat{\phi} = \frac{\text{Déviance résiduelle}}{n – k}$$
Si ce ratio s’avère nettement supérieur à $1,0$ (dépassant couramment des seuils de $1,5$ ou $2,0$), l’application aveugle de glm(..., family = poisson) engendre une sous-estimation massive des erreurs-types des coefficients, entraînant une inflation majeure du risque d’erreur de première espèce (rejet erroné de l’hypothèse nulle d’absence d’effet).

Pour remédier à cette défaillance sans délaisser le cadre des GLM, R offre une solution méthodologique robuste via la quasi-vraisemblance, accessible par family = quasipoisson(link = "log"). Ce modèle conserve strictement inchangées les estimations ponctuelles des coefficients $\hat{\beta}$, mais introduit un paramètre de dispersion libre $phi$ qui gonfle proportionnellement la matrice de variance-covariance des paramètres :
$$Var(Y|X) = \phi \mu$$
Lorsque la surdispersion est particulièrement sévère ou découle d’un mécanisme biologique sous-jacent spécifique, l’analyste peut recourir à la modélisation binomiale négative (via la fonction glm.nb() du package MASS) qui postule une variance quadratique de la forme $Var(Y) = \mu + \alpha \mu^2$, ou déployer des modèles à inflation de zéros (Zero-Inflated Poisson) pour dissocier le processus d’absence absolue de comportement du processus de comptage contingent.

8. Diagnostics de modèles, évaluation de l’ajustement et sélection

8.1 Typologie comparée des résidus : ordinaires, standardisés, de Pearson et de déviance

L’évaluation rigoureuse de l’adéquation d’un modèle statistique ajusté passe impérativement par l’analyse minutieuse de ses résidus. Dans le contexte de la régression linéaire classique sous lm(), la notion de résidu est univoque : le résidu ordinaire, ou brut, représente la différence euclidienne élémentaire entre l’observation et la valeur ajustée ($e_i = y_i – \hat{y}_i$). Pour corriger l’effet des points de fort levier, ces résidus peuvent être divisés par l’écart-type de leur distribution d’échantillonnage pour constituer les résidus standardisés ou studentisés.

Dans l’architecture des GLM régie par glm(), la variabilité conditionnelle de la variable réponse dépendant intrinsèquement de la moyenne, la simple différence arithmétique brute $y_i – \hat{\mu}_i$ est dépourvue de propriétés comparables entre différentes zones de l’espace des prédicteurs. R implémente donc une typologie diversifiée de résidus répondant à des fonctions diagnostiques spécialisées :

• Résidus de Pearson : Définis par la déviation brute pondérée par la racine carrée de la fonction de variance :
  $$r_{P, i} = \frac{y_i – \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)}}$$
  Ces résidus mesurent l’écart relatif en unités de variance attendue et servent de base directe au calcul de la statistique de Pearson globale.
• Résidus de déviance : Définis comme la contribution signée de chaque observation individuelle à la fonction de déviance globale du modèle :
  $$d_i = \text{sign}(y_i – \hat{\mu}_i) \sqrt{2 \left[ ell(y_i; y_i) – ell(\hat{\mu}_i; y_i) \right]}$$
  où $ell(y_i; y_i)$ représente la log-vraisemblance maximale atteignable dans un modèle saturé parfait. Les résidus de déviance présentent un comportement distributionnel beaucoup plus proche de la normalité asymptotique, en faisant l’outil privilégié pour la vérification graphique de l’ajustement.

Pour le calcul des mesures d’influence et des leviers, les fonctions diagnostiques de R adaptent la formulation classique du chapeau matriciel (hat matrix) en incorporant la matrice des poids de travail finaux $W$, ce qui permet d’isoler avec exactitude les valeurs aberrantes ayant un poids disproportionné sur les algorithmes d’ajustement non linéaires.

8.2 Métriques d’ajustement global : R² ajusté contre Déviance et AIC

L’une des distinctions les plus évidentes entre les affichages terminaux de summary.lm() et summary.glm() réside dans l’abandon systématique du coefficient de détermination $R^2$ au profit du concept de déviance. Dans le modèle linéaire classique, le $R^2$ et son pendant ajusté quantifient avec précision la proportion de la variance totale de la réponse expliquée par l’hyperplan de régression, en vertu de l’orthogonalité exacte de la décomposition de la somme des carrés conférée par le théorème de Pythagore dans l’espace euclidien.

Cette décomposition additive de la variance devient mathématiquement caduque dès lors que l’on manipule des distributions non gaussiennes et des fonctions de lien non linéaires. Dès lors, summary.glm() ne fournit délibérément aucun coefficient $R^2$. Tenter de calculer un $R^2$ ordinaire sur les valeurs brutes d’une régression logistique ou poissonienne constitue une erreur conceptuelle manifeste. Les analystes peuvent mobiliser des pseudo-$R^2$ (tels que ceux de McFadden, Cox & Snell, ou Nagelkerke), mais ces indices n’offrent aucune propriété d’orthogonalité stricte et doivent être interprétés avec une grande réserve.

Dans glm(), la mesure cardinale de l’inadéquation globale est la déviance résiduelle, qui mesure la distance statistique séparant le modèle ajusté du modèle saturé (modèle parfait allouant autant de paramètres qu’il existe d’observations). La déviance résiduelle s’analyse comparativement à la déviance nulle, représentative d’un modèle ne comportant que l’interception constante. Pour la sélection de modèles concurrents emboîtés ou non emboîtés, la pratique méthodologique repose prioritairement sur les critères de pénalisation de la complexité : le critère d’information d’Akaike (AIC) et le critère bayésien de Schwarz (BIC) :
$$\text{AIC} = -2ell(\hat{\beta}) + 2k$$
Ces métriques permettent d’arbitrer rigoureusement entre le gain d’ajustement à la log-vraisemblance et la surparamétrisation du modèle prédictif.

8.3 Tests d’hypothèses globales : test F de Fisher versus test du rapport de vraisemblance

L’évaluation statistique globale de l’hypothèse nulle stipulant que l’ensemble des coefficients de régression (hormis l’interception) sont simultanément nuls s’opère selon des logiques différentes dans les deux architectures. Dans le résultat d’un appel à lm(), R délivre de façon canonique la statistique du test $F$ de Fisher-Snedecor. Ce test bénéficie d’une distribution d’échantillonnage exacte sous l’hypothèse de normalité des résidus :
$$F = \frac{(\text{SST} – \text{SSR}) / p}{\text{SSR} / (n – p – 1)}$$
où $\text{SST}$ représente la somme totale des carrés et $\text{SSR}$ la somme des carrés résiduels. Ce ratio garantit un contrôle parfait du risque d’erreur de première espèce en échantillon fini, pourvu que les postulats gaussiens soient scrupuleusement honorés.

Dans l’univers des GLM, le test $F$ universel n’est théoriquement disponible que dans le cas particulier gaussien ou quasi-vraisemblable avec estimation de dispersion. La stratégie canonique d’évaluation globale s’appuie sur le Test du Rapport de Vraisemblance (LRT pour Likelihood Ratio Test), formalisé sous R via l’analyse de déviance par la méthode anova(modele, test = "Chisq"). Ce test repose sur le théorème de Wilks, démontrant que la différence de déviance entre deux modèles emboîtés converge asymptotiquement vers une loi du Chi-deux dont les degrés de liberté correspondent à la différence du nombre de paramètres estimés :
$$\Delta D = D_{\text{réduit}} – D_{\text{complet}} \sim \chi^2(\Delta k)$$

À côté du LRT, qui impose l’ajustement consécutif de deux modèles distincts, figurent le test de Wald et le test du Score (Lagrange Multiplier). Si le test de Wald est directement inspectable via les erreurs-types de la sortie univariée de summary.glm(), il est largement documenté dans la littérature statistique que le test du rapport de vraisemblance présente une puissance statistique systématiquement supérieure et une robustesse accrue, particulièrement dans les échantillons de taille modérée ou en présence de configurations extrêmes où le test de Wald souffre de l’aberration de Hauck-Donner (écrasement pathologique de l’erreur-type pour de très grands effets).

9. Inférence statistique, incertitude et tests d’hypothèses spécifiques

9.1 Distributions d’échantillonnage : test t de Student contre statistique z de Wald

L’inférence univariée sur la significativité statistique d’un coefficient individuel $\hat{\beta}_j$ illustre de manière concrète la dichotomie entre inférence exacte en échantillon fini et inférence asymptotique en grand échantillon. Lorsque le chercheur inspecte la table des coefficients générée par summary.lm(), la quatrième colonne affiche sans exception une statistique $t$ associée à une p-valeur issue d’une loi de Student à $n – k$ degrés de liberté. Cette distribution d’échantillonnage découle du rapport entre l’estimateur gaussien du coefficient et l’estimateur sans biais de son erreur-type basé sur la racine carrée de la variance résiduelle exacte :
$$t = \frac{\hat{\beta}_j – 0}{\text{SE}(\hat{\beta}_j)} \sim \mathcal{T}(n – k)$$

À l’inverse, l’affichage généré par summary.glm() pour les familles non gaussiennes (telles que les familles binomiale et de Poisson) délaisse le test $t$ pour présenter une statistique $z$ adossée à une loi normale centrée réduite asymptotique $\mathcal{N}(0, 1)$ :
$$z = \frac{\hat{\beta}_j – 0}{\text{SE}_{asy\mp}(\hat{\beta}_j)} \sim \mathcal{N}(0, 1)$$
Ce recours à la loi normale unitaire découle directement des propriétés de normalité asymptotique des estimateurs du maximum de vraisemblance, postulées par le théorème central limite sous les conditions de régularité de Cramér-Rao.

Cette différence induit des conséquences directes pour l’expérimentateur travaillant sur de petits échantillons en psychologie cognitive ou en neurosciences. Alors que la loi de Student prend en compte l’incertitude supplémentaire liée à l’estimation conjointe du paramètre d’échelle inconnu $\sigma$ en épaississant les queues de distribution, la statistique $z$ d’un GLM ignore cet ajustement en postulant implicitement que la taille de l’échantillon approche de l’infini. En conséquence, sur de très faibles effectifs ($n < 30$), les tests de Wald univariés rapportés par summary.glm() ont tendance à fournir des p-valeurs optimistes, nécessitant une vigilance accrue et l’éventuel recours à des méthodes d’inférence par bootstrap ou permutations.

9.2 Matrices de variance-covariance et robustesse aux violations

L’évaluation de l’incertitude entourant le vecteur des estimateurs s’incarne dans la matrice de variance-covariance des coefficients, extractible sous R à l’aide de la fonction générique vcov(). Dans un modèle lm(), cette matrice est le produit scalaire direct de la variance de résidu estimée et de la matrice d’expérience inversée :
$$\text{Cov}(\hat{\beta}) = s^2 (X^T X)^{-1}$$
Dans l’environnement glm(), la covariance est calculée par inversion de la matrice d’information de Fisher attendue, pondérée par les poids de l’algorithme IRLS :
$$\text{Cov}(\hat{\beta}) = \phi (X^T W X)^{-1}$$
où $phi$ est le paramètre de dispersion (fixé à l’unité dans les distributions de Poisson et binomiale pures).

Cette formulation révèle l’extrême sensibilité des erreurs-types des GLM à la spécification correcte de la fonction de variance $V(\mu)$. Si la structure de moyenne est parfaitement modélisée mais que la dispersion réelle des données diverge des hypothèses théoriques de la famille déclarée, les erreurs-types standards issues de vcov() se trouvent structurellement biaisées, invalidant les tests de significativité.

Pour immuniser l’inférence statistique contre ces violations d’hétéroscédasticité ou de mauvaise spécification de la dispersion, les méthodologistes mobilisent des estimateurs sandwich de variance-covariance hétéroscédastiques-consistants (HC), largement popularisés par le package sandwich de R. En couplant ces estimateurs robustes (tels que HC3 ou HC4) avec la fonction coeftest() du package lmtest, il devient possible de mener des inférences valides sur les paramètres d’un modèle lm() ou glm() même lorsque les postulats distributionnels d’homogénéité sont substantiellement enfreints.

9.3 Calcul des intervalles de confiance et prédictions probabilistes

L’estimation de l’incertitude prédictive constitue un point de divergence opérationnel majeur entre les méthodes predict.lm() et predict.glm(). Dans le modèle linéaire classique, le prédicteur et l’espérance de la variable dépendante évoluent dans le même espace euclidien continu non borné. Le calcul d’un intervalle de confiance pour l’espérance conditionnelle via predict(modele_lm, interval = "confidence") s’effectue par une simple combinaison symétrique de la marge d’erreur autour de la valeur ajustée : $\hat{y} \pm t_{alpha/2} \cdot \text{SE}$.

Dans un modèle linéaire généralisé, le calcul d’intervalles de confiance impose une rigueur méthodologique bien plus sophistiquée en raison de la non-linéarité de la fonction de lien. Une erreur couramment constatée chez les praticiens consiste à calculer les erreurs-types directement sur l’échelle de la réponse brute, générant des intervalles de confiance symétriques susceptibles de dépasser les bornes logiques (par exemple des bornes inférieures de probabilités négatives pour un GLM logistique). La démarche rigoureuse consiste impérativement à extraire les prédictions et leurs erreurs-types sur l’échelle du prédicteur linéaire non contraint via :
$$pred = \text{predict}(modele_glm, type = « link », se.fit = TRUE)$$
L’intervalle de confiance est d’abord calculé de manière symétrique sur cette échelle latente :
$$\eta_{\text{\inf, \sup}} = \hat{\eta} \pm z_{alpha/2} \cdot \text{SE}_{\eta}$$
Ce n’est qu’ensuite que les bornes résultantes sont projetées sur l’espace d’observation naturel par l’application de la fonction de lien réciproque :
$$\mu_{\text{\inf, \sup}} = g^{-1}(\eta_{\text{\inf, \sup}})$$

Cette séquence mathématique produit des intervalles de confiance intrinsèquement asymétriques sur l’échelle des probabilités ou des taux, reflétant la compression progressive de l’incertitude statistique au fur et à mesure que la moyenne s’approche des bornes asymptotiques structurelles du domaine de validité empirique.

10. Architecture logicielle sous R : anatomie des classes d’objets lm et glm

10.1 Hiérarchie et héritage dans le système de classes S3 de R

L’écosystème de programmation statistique R s’appuie de manière prépondérante sur le système de programmation orientée objet informel S3. L’inspection des métadonnées de classe révèle une structure d’héritage fondamentale pour la compréhension du fonctionnement logiciel : un objet généré par l’exécution de la fonction glm() possède un vecteur d’appartenance de classe composé de deux éléments : c("glm", "lm"). Cette double déclaration signifie explicitement que la classe glm hérite directement de la classe mère lm.

En vertu de ce mécanisme d’héritage polymorphique, lorsqu’une fonction générique (telle que residuals(), coefficients() ou model.matrix()) est invoquée sur un objet produit par glm(), le dispatcheur de méthodes S3 recherche initialement l’existence d’une variante spécialisée portant le suffixe .glm. Si cette méthode spécifique n’est pas implémentée dans l’espace de noms des packages chargés, R délègue automatiquement l’exécution à la méthode générique de niveau supérieur portant le suffixe .lm, préservant une cohérence fonctionnelle exemplaire.

Néanmoins, l’anatomie interne de la liste chaînée sous-jacente d’un objet glm intègre des composants additionnels inexistants dans un simple objet lm. Parmi ces attributs cruciaux figurent la composante family (encapsulant l’ensemble des mécanismes différentiels de la distribution), la valeur scalaire de deviance, le critère aic, la mesure de null.deviance, ainsi que le vecteur logique converged indiquant si l’algorithme des moindres carrés itérativement repondérés a validé son critère d’arrêt numérique avant d’atteindre le plafond d’itérations alloué.

10.2 Comportement différentiel des méthodes génériques

Malgré cette communauté d’héritage, le comportement polymorphique de certaines méthodes génériques universelles présente des bifurcations fondamentales selon que l’instance traitée est un objet purement lm ou un objet dérivé glm. L’exemple le plus prégnant s’incarne dans la méthode summary() : summary.lm() mobilise un sous-programme compilé destiné à évaluer l’estimateur sans biais de $\sigma$ et construit une table de contrastes adossée aux distributions de Student et de Fisher, tandis que summary.glm() exécute l’analyse asymptotique du maximum de vraisemblance, calcule les déviances relatives et sollicite les lois du Chi-deux et de la normale centrée réduite.

De façon tout aussi marquante, la fonction graphique générique de diagnostic plot() active des comportements visuels contrastés :

• Dans plot.lm(), les quatre graphiques diagnostiques standards représentent les résidus ordinaires studentisés en regard des valeurs prédites, le diagramme quantile-quantile de normalité des erreurs brutes, l’échelle de dispersion-localisation homoscédastique, et l’espace de levier pondéré par les distances de Cook.
• Dans plot.glm(), R substitue automatiquement aux résidus ordinaires les résidus de déviance standardisés, adaptant la lecture des écarts extrêmes à la courbure de la fonction de lien et à la variance attendue sous la famille sélectionnée, évitant ainsi d’induire le chercheur en erreur sur la prétendue non-normalité d’un phénomène purement catégoriel.

Ces divergences de comportement polymorphique s’étendent aux packages avancés d’inférence et d’effets marginaux, tels que emmeans ou effects, qui doivent impérativement identifier si le modèle d’origine opère sur une échelle transformée par un lien non linéaire afin de moduler l’échelle d’affichage des moyennes marginales ajustées présentées à l’utilisateur.

10.3 Extraction programmatique et intégration dans des pipelines de calcul

Dans l’ingénierie contemporaine de données sous R, l’interfaçage programmatique des sorties de modélisation avec les structures de données rectangulaires du tidyverse repose quasi-exclusivement sur les verbes normalisés du package broom : tidy(), glance() et augment(). L’utilisation de ces verbes révèle les spécificités de la structure sous-jacente de chaque modèle.

Appliquée à un modèle issu de lm(), la fonction glance() retourne un tableau à ligne unique renseignant des métriques purement linéaires : $R^2$, $R^2$ ajusté, statistique $F$, et variance de l’erreur résiduelle standard. Lorsqu’elle traite un modèle issu de glm(), glance() purge son affichage de tout indicateur de variance expliquée euclidienne pour renseigner la déviance résiduelle, la déviance du modèle nul et l’AIC. La tentative d’automatiser l’extraction d’une colonne r.squared au sein d’une boucle itérative mixant indifféremment des modèles linéaires et des modèles logistiques provoquera inévitablement une interruption d’exécution pour cause de valeur manquante ou nulle.

Par ailleurs, dans le cadre de déploiements d’architectures de calcul intensif impliquant des rééchantillonnages massifs de type bootstrap ou des validations croisées par k-blocs, les objets complets de classe lm et glm présentent une volumétrie mémoire disproportionnée en raison du stockage par défaut de l’environnement parent, du cadre de modélisation intégral (model.frame) et des matrices de coordonnées. L’optimisation algorithmique impose alors d’élaguer les composants superflus des objets retournés en assignant la valeur NULL aux slots volumineux, ou en invoquant explicitement l’argument model = FALSE, y = FALSE, x = FALSE lors de l’appel initial, réduisant l’empreinte mémoire d’un facteur pouvant atteindre 95 % sans altérer la capacité de prédiction.

11. Démonstrations empiriques appliquées à la recherche en psychologie

11.1 Étude de cas 1 : Modélisation d’un score de régulation émotionnelle continu

Pour illustrer la convergence empirique absolue entre les moindres carrés ordinaires et le maximum de vraisemblance gaussien, considérons une recherche expérimentale en psychologie cognitive examinant les déterminants d’un score de régulation émotionnelle continue (mesuré via une échelle standardisée valide de 10 à 100 points, approximant une distribution continue symétrique). L’échantillon comporte 200 participants répartis entre une condition contrôle et un protocole d’entraînement attentionnel basé sur la pleine conscience, l’âge des participants étant intégré comme covariable continue.

L’estimation parallèle des deux formulations s’exprime syntaxiquement ainsi :
$$modele_lm = \text{lm}(Score \sim Age + Traitement, data = data_psy)$$
$$modele_glm = \text{glm}(Score \sim Age + Traitement, family = gaussian(link = « identity »), data = data_psy)$$
L’extraction comparative des matrices de coefficients fournit les résultats strictement superposables récapitulés dans le tableau suivant :

Paramètre estimé	Coefficient lm()	SE lm()	Coefficient glm(gaussian)	SE glm(gaussian)
Interception ($\beta_0$)	45,2134	2,1451	45,2134	2,1451
Âge du participant ($\beta_1$)	0,3142	0,0423	0,3142	0,0423
Traitement actif ($\beta_2$)	8,7651	1,2130	8,7651	1,2130

L’examen minutieux des décimales de calcul confirme l’égalité des estimations ponctuelles et des erreurs-types jusqu’au 12e ordre de grandeur après la virgule. L’inspection des diagnostics résiduels confirme l’homoscédasticité et la normalité approximative de la dispersion stochastique, validant empiriquement l’inférence selon laquelle le protocole d’entraînement élève en moyenne le score d’aptitude de près de 8,8 points à âge contrôlé, la significativité statistique ($p < 0,001$) étant rigoureusement identique dans les deux environnements de test.

11.2 Étude de cas 2 : Prédiction de l’abandon thérapeutique (attrition clinique binaire)

Considérons à présent une cohorte clinique en psychopathologie évaluant le phénomène d’abandon thérapeutique précoce (défini par une variable dichotomique : $1 =$ le patient interrompt prématurément le suivi psychothérapeutique ; $0 =$ le patient mène l’intervention à son terme). La cohorte regroupe 350 patients, et la variable prédictive majeure est le score initial d’alliance thérapeutique perçue lors du premier entretien d’admission.

Si un analyste ajuste naïvement un modèle de probabilité linéaire via :
$$fit_l\pm = \text{lm}(Attrition \sim Alliance, data = data_clinique)$$
l’observation des prédictions générées aux extrémités révèle la défaillance structurelle de l’outil : pour des scores d’alliance thérapeutique extrêmement élevés observés chez certains participants, le modèle produit des probabilités d’abandon négatives allant jusqu’à $\hat{Y} = -0,14$, une absurdité théorique et clinique évidente.

La bascule méthodologique vers une formulation adaptée s’opère par :
$$fit_logit = \text{glm}(Attrition \sim Alliance, family = binomial(link = « logit »), data = data_clinique)$$
Le coefficient estimé pour l’alliance thérapeutique ressort à $\hat{\beta} = -0,085$ ($z = -4,82, p < 0,0001$). En appliquant l'exponentiation, on obtient un rapport de cotes : $$\text{OR} = \exp(-0,085) = 0,918$$
L’interprétation clinique s’énonce ainsi : pour chaque incrément d’un point sur l’échelle d’alliance thérapeutique initiale, la cote (odds) d’abandon prématuré du traitement diminue de manière significative de 8,2 % ($1 – 0,918$). Plus important encore, les prédictions individuelles obtenues via predict(fit_logit, type = "response") demeurent harmonieusement asymptotiques, approchant de zéro sans jamais le franchir, respectant parfaitement la structure d’une distribution de probabilité conditionnelle.

11.3 Étude de cas 3 : Analyse des épisodes compulsifs quotidiens (données de comptage)

Le troisième scénario concerne une investigation en psychiatrie translationnelle analysant la fréquence quotidienne de manifestations compulsives chez des patients atteints de troubles obsessionnels compulsifs (TOC), mesurée sur une période de 30 jours consécutifs à l’issue de l’administration d’une molécule expérimentale versus un placebo. Les données sont des comptages entiers discrets présentant une très forte asymétrie positive, avec une moyenne d’épisodes compulsifs modérée mais certains individus affichant des valeurs très élevées.

L’estimation initiale par régression linéaire lm(Episodes ~ Molecule + Severite_Initiale) produit des résidus hautement hétéroscédastiques et non gaussiens (test de Shapiro-Wilk : $p < 10^{-6}$). De surcroît, le modèle linéaire prédit des nombres moyens négatifs de compulsions pour des patients présentant une sévérité initiale faible sous traitement expérimental.

L’application consécutive de la régression de Poisson via glm(Episodes ~ Molecule + Severite_Initiale, family = poisson) stabilise l’espace prédictif grâce à la fonction de lien logarithmique. Cependant, l’évaluation de la dispersion résiduelle livre une déviance résiduelle de 1420 pour 397 degrés de liberté, révélant un ratio de surdispersion massif :
$$\hat{\phi} = \frac{1420}{397} \approx 3,58$$
Les erreurs-types calculées sous le modèle de Poisson pur sont donc artificiellement réduites, gonflant indûment la significativité statistique apparente des traitements.

La correction définitive s’obtient en mobilisant le modèle quasi-poissonnien :
$$fit_quasi = \text{glm}(Episodes \sim Molecule + Severite_Initiale, family = quasipoisson)$$
L’estimation ponctuelle du paramètre de la molécule expérimentale demeure strictement identique ($\hat{\beta} = -0,412$), ce qui traduit un Ratio de Taux d’Incidence de :
$$\text{IRR} = \exp(-0,412) = 0,662$$
soit une réduction attendue de 33,8 % de la fréquence des compulsions quotidiennes par rapport au bras placebo. En revanche, l’erreur-type ajustée passe de $0,061$ sous le modèle de Poisson classique à $0,061 \times \sqrt{3,58} = 0,115$ sous le modèle quasi-poisson. Bien que l’effet thérapeutique demeure statistiquement significatif ($t = -3,58, p < 0,001$), l'élargissement de l'intervalle de confiance protège le chercheur contre les conclusions fallacieuses dictées par une hypothèse d'équidispersion invalidée par les données empiriques.

12. Guide décisionnel et recommandations méthodologiques pour la pratique de recherche

12.1 Algorithme de décision séquentiel selon la distribution de la variable cible

Afin de systématiser l’arbitrage méthodologique entre les fonctions lm() et glm() dans la pratique quotidienne de la modélisation statistique sous R, l’analyste peut déployer un protocole décisionnel séquentiel articulé autour de l’échelle de mesure et des propriétés distributionnelles de la variable cible. Le parcours logique s’organise selon les étapes critiques suivantes :

• Étape 1 : Qualification du support d’échelle de la variable dépendante $Y$
  • Si $Y$ est continue, théoriquement non bornée, et présente une distribution d’erreur symétrique sans contrainte de positivité absolue : privilégier l’usage canonique de lm().
  • Si $Y$ est de nature discrète, catégorielle, bornée ou strictement positive : rejeter lm() et s’orienter impérativement vers l’architecture glm().
• Étape 2 : Traitement des données binaires et proportionnelles
  • Pour $Y in {0, 1}$ (succès/échec, vivant/décédé) : mobiliser glm(..., family = binomial(link = "logit")). Si l’assise théorique postule un processus latent continu normalement distribué, basculer vers le lien link = "probit".
  • Pour des proportions agrégées issues de dénombrements ($k$ succès sur $n$ essais) : formuler la réponse sous forme matricielle à deux colonnes cbind(succes, echecs) et invoquer la famille binomiale.
• Étape 3 : Traitement des données de comptage et de temps d’attente
  • Pour des entiers positifs dénombrant des événements ($\mathbb{N}$) : initialiser avec glm(..., family = poisson(link = "log")).
  • Vérifier systématiquement le paramètre de dispersion. En cas de surdispersion constatée ($\hat{\phi} > 1,2$), ajuster via family = quasipoisson ou basculer vers un modèle binomial négatif (MASS::glm.nb).
  • Pour des temps d’attente continus positifs asymétriques : tester l’adéquation d’une distribution Gamma avec family = Gamma(link = "log").

Cet algorithme prévient les erreurs de spécification les plus courantes tout en garantissant que l’espace des valeurs prédites demeure en adéquation rigoureuse avec la géométrie de la variable empirique sous étude.

12.2 Bonnes pratiques de rapportage statistique selon les normes APA

La transcription académique des résultats issus de modélisations quantitatives doit répondre aux exigences de transparence méthodologique stipulées par la septième édition du manuel de l’American Psychological Association (APA). Lorsqu’un modèle lm() est déployé, il est impératif de rapporter systématiquement la statistique globale d’ajustement du modèle ($F$ avec ses deux degrés de liberté associés, la p-valeur correspondante, et le $R^2$ ou $R^2$ ajusté), ainsi que les coefficients bêtas standardisés ou non standardisés accompagnés de leurs intervalles de confiance à 95 % et de la statistique $t$ univariée.

Lorsque l’analyse fait appel à l’architecture glm(), les impératifs de description changent de nature et exigent une traçabilité accrue :

• Déclaration explicite de la structure stochastique : Mentionner sans ambiguïté la famille distributionnelle retenue ainsi que la fonction de lien utilisée (par exemple : « Régression logistique binaire via une fonction de lien logit »).
• Indicateurs d’ajustement global : Rapporter la déviance nulle et la déviance résiduelle en spécifiant les degrés de liberté résiduels correspondants, ou comparer les modèles concurrents à l’aide de l’AIC et des tests de rapport de vraisemblance ($\chi^2(df), p$).
• Rapportage des coefficients : Présenter les estimations non standardisées du prédicteur linéaire ($\beta$), leurs erreurs-types asymptotiques, la statistique de Wald ($z$), ainsi que l’indice de taille d’effet transformé adapté au contexte : Rapport de Cotes (OR) ou Rapport de Taux d’Incidence (IRR), systématiquement accompagnés de leurs intervalles de confiance à 95 % calculés sur l’échelle du lien puis retransformés.

Enfin, la politique contemporaine de science ouverte impose de déposer dans un répertoire public (tel que l’Open Science Framework) les scripts R intégraux documentant l’appel des fonctions, les éventuelles configurations manuelles d’optimisation via glm.control, et les fonctions de manipulation de contrastes, gage incontournable de la reproductibilité computationnelle des résultats publiés.

12.3 Transitions vers des modèles plus complexes : GAM, GLMM et modélisation multiniveau

Bien que les fonctions canoniques lm() et glm() constituent les pierres angulaires de l’analyse statistique sous R, leur cadre analytique postule formellement l’indépendance stricte des erreurs résiduelles d’échantillonnage. Dès lors que le devis de recherche intègre des mesures répétées dans le temps, des passations longitudinales, ou une structure hiérarchique en grappes (élèves emboîtés dans des classes, patients hébergés dans différents centres hospitaliers), ces deux fonctions deviennent conceptuellement inadaptées en raison de la violation du postulat d’indépendance.

Dans ce contexte, le prolongement direct de la régression linéaire lm() s’incarne dans les modèles linéaires mixtes (LMM), opérationnalisés dans R par la fonction lmer() du package lme4. Ces modèles décomposent la variance résiduelle en introduisant des effets aléatoires (interceptions et pentes aléatoires) capables de modéliser avec précision la covariance intra-sujet ou intra-grappe.

Symétriquement, l’extension naturelle des modèles linéaires généralisés de glm() vers les structures de données hiérarchiques se matérialise par les Modèles Linéaires Généralisés Mixtes (GLMM), mobilisables via la fonction glmer() de lme4. Les GLMM réalisent la synthèse conceptuelle ultime en combinant la flexibilité des fonctions de lien de la famille exponentielle avec la modélisation de l’hétérogénéité aléatoire multiniveau. Par ailleurs, lorsque les relations fonctionnelles liant les prédicteurs au terme de réponse défient toute spécification paramétrique simple, l’analyste peut franchir le pas vers les Modèles Additifs Généralisés (GAM) via le package mgcv, qui substitue aux coefficients linéaires rigides des fonctions de lissage non paramétriques (splines).

En maîtrisant la dialectique entre lm() et glm(), le chercheur acquiert bien plus qu’une compétence technique d’exécution de code sous R : il intègre les principes structurels fondamentaux régissant l’inférence statistique moderne, garantissant que ses choix méthodologiques soutiennent avec rigueur et fidélité l’avancement des connaissances empiriques.

Références

• Agresti, A. (2015). Foundations of linear and generalized linear models. John Wiley & Sons.
• American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000
• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
• Chambers, J. M., & Hastie, T. J. (Eds.). (1992). Statistical models in S. Wadsworth & Brooks/Cole Advanced Books & Software.
• Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). SAGE Publications.
• Hastie, T. J., & Tibshirani, R. J. (1990). Generalized additive models. Chapman & Hall/CRC.
• McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman & Hall/CRC. https://doi.org/10.1007/978-1-4899-3242-6
• Nelder, J. A., & Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A (General), 135(3), 370–384. https://doi.org/10.2307/2344614
• R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
• Robinson, D., Hayes, A., & Couch, S. (2023). broom: Convert statistical analysis objects into tidy tibbles (R package version 1.0.5). https://CRAN.R-project.org/package=broom
• Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S (4th ed.). Springer. https://doi.org/10.1007/978-0-387-21706-2
• Wood, S. N. (2017). Generalized additive models: An introduction with R (2nd ed.). Chapman & Hall/CRC. https://doi.org/10.1201/9781315370279
• Zeileis, A. (2004). Econometric computing with HC and HAC covariance matrix estimators. Journal of Statistical Software, 11(10), 1–17. https://doi.org/10.18637/jss.v011.i10