Les bases du langage

Introduction

Qu'est-ce que R ?

• R est un logiciel pour le traitement et la modélisation de données statistiques (Ihaka & Gentleman, 1996; Venables & Ripley, 2002).
• Dans cet environnement interactif, l'utilisateur envoie des commandes, R les interprète et renvoit un résultat.
• Il existe des interfaces graphiques plus ou moins avancées (Deducer, Rcmdr, etc.), mais on ne les utilisera pas.

R est gratuit, possède de nombreuses fonctionnalités pour la modélisation, le reporting et la visualisation de données. Il s'agit avant tout d'un langage de programmation inspiré du langage S développé dans les années 80 (Becker & Chambers, 1981; 1984). Le projet est maintenu par la R Foundation for Statistical Computing, et toutes les informations officielles sur le projet R se trouvent à cette adresse : http://www.r-project.org.

Représentation des données

Contrairement à SPSS ou Statistica, on ne travaille pas à partir d'une vue de type tableur des données, comme .

Un tableau Excel sera représenté sous forme d'un data.frame, c'est-à-dire un tableau avec des données de type mixte (chaînes de caractères, nombres, dates, etc.) et des lignes/colonnes nommées.

Des sources de données externes au format SPSS, Stata, SAS, ou des bases de données MySQL, peuvent être importées sans problème majeur.

Exemple de commandes R

Voici un exemple de régression simple à partir de données artificielles :

n <- 100
x <- runif(n)
y <- 0.9 + 0.65*x + rnorm(n)
## Fit a simple linear model
summary(lm(y ~ x))

Dans l'exemple ci-dessus, on voit que la syntaxe est proche de celle d'un langage de programmation : on travaille sur des variables à l'aide de fonctions.

Comme dans tout langage, il existe un vocabulaire et une syntaxe spécifique au langage R, ainsi que des règles sémantiques. Dans une optique introductive, on s'intéressera essentiellement aux fonctions de base qui permettent d'effectuer la plupart des opérations statistiques.

Environnement de travail

Pour travailler avec R, il est nécessaire d'avoir les outils adaptés, en particulier un environnement de développement "intelligent" qui ne gêne pas l'approche interactive, p.ex. RStudio : http://rstudio.org/.

Obtenir de l'aide

Il est indispensable de savoir où chercher l'information et de bien se documenter !

• sous R directement, help(topic, package=) ou help.search(pattern)
• la documentation R disponible sur CRAN (Manuals)
• la mailing-list R-help
• le moteur de recherche Google R-Seek

• les sites Stack Exchange (bien lire la FAQ) :

  • http://stackoverflow.com
  • http://crossvalidated.com

The R Inferno (Burns, 2011) est une lecture conseillée pour les personnes désirant approfondir leurs connaissances de R.

Ressources bibliographiques

En anglais :

• Introductory Statistics with R (Dalgaard, 2008)
• Biostatistical Design and Analysis Using R (Logan, 2010)
• Modern Applied Statistics with S (Venables & Ripley, 2002)

En français :

• R pour les débutants (Paradis, 2005)
• Comprendre et utiliser les statistiques dans les sciences de la vie (Falissard, 2005)

Voir aussi http://www.R-project.org/doc/bib/R-books.html.

Les objets R

L'élément de base est le vecteur (un scalaire est alors un vecteur de longueur 1) que l'on peut créer avec la commande c et dont les éléments sont de l'un des "types" (mode) suivants :

• numeric, des nombres (entiers ou flottants, integer/double)
• logical, des booléens à valeur dans {TRUE, FALSE}
• character, des caractères ASCII

Les expressions suivantes sont valides :

nombre <- 3.141593
v <- c(1,2,3,4)
b <- c(T,F,T,F)
s <- c("h","e","l","l","o")
a <- b <- rnorm(2)
u <- b  

Adresser les éléments d'un vecteur

On utilise [ pour indexer les \(n\) éléments d'un vecteur, sachant que l'index du premier élément vaut 1 (et non 0 comme dans certains langages) :

v[1]
v[c(1,4)]

On peut adresser les élements d'un vecteur en utilisant les valeurs contenues dans un autre vecteur (principe d'un dictionnaire), par exemple :

v[b]    ## valeurs de v telles que b vaut TRUE
s[v[3]] ## v[3]=3ème valeur de s

On peut encore adresser les éléments d'un vecteur par leur "nom", le cas échéant :

names(v) <- c("toto","titi","tata","tutu")
v["tata"]

Application

On peut également vouloir isoler des groupes de valeurs :

• Restriction d'un vecteur à certaines de ses valeurs (filtre),
• ventilation de scores par niveaux d'un facteur de classification.

Exemple :

n <- 100
age <- sample(25:45, n, replace=TRUE)
sex <- sample(c("M","F"), n, replace=TRUE)
length(age[sex=="F"])
length(sex[age >= 30 & age <= 40])

Une expression telle que age > 30 retourne un vecteur booléen ("les éléments du vecteur age sont-ils supérieurs (strictement) à 30"). Ce vecteur a la même longueur que age. Les indices de age pour lesquels cette condition est vérifiée peuvent être obtenus ainsi :

which(age > 30)

Recap'

Voici deux façons d'adresser spécifiquement un sous-ensemble des valeurs assignées à un vecteur x :

x <- c(1,4,10,3,1)
names(x) <- letters[1:length(x)]
idx <- c(1,3,4)
g <- c(T,F,T,T,F)
all(x[idx] == x[g])

avec les caractéristiques suivantes :

• length(idx) \(\le\) length(x)
• length(g) \(=\) length(x)

Propriétés des vecteurs

On peut connaître la taille d'un vecteur (length), son nombre d'éléments uniques (unique), le type de ses éléments (class ou str).

On peut également nommer les éléments d'un vecteur :

x <- rnorm(100)
res <- c(m=mean(x), s=sd(x))

On peut les "assembler" ou effectuer des opérations d'algèbre linéaire :

y <- rnorm(100)
rbind(x, y)
cbind(x, y)
t(x) %*% x * 1/(length(x)-1)

Générer des vecteurs

En dehors de l'opérateur c, on peut utiliser les commandes pré-définies dans R :

1:5
seq(1, 5)
rep(c(1,3), 5)
rep(c(1,3), each=5)  ## rep(c(1,3), c(5,5))

Comme on l'a vu, l'indexation d'un vecteur repose sur... un vecteur :

x <- c(0,1,1,0,1,0,1,0,0)
x[3]
x[length(x)]
x[2:4]
x[c(2,4)]
seq(1, 10)[seq(1, 10, by=2)]  ## nombres impairs
seq(1, 10)[seq(2, 10, by=2)]  ## nombres pairs

Opérer sur des vecteurs

On peut également effectuer des opérations de tri sur un vecteur :

x <- 1:10
xs <- sample(x)
sort(xs, decreasing=TRUE)
order(xs)

ou obtenir le rang de ses élements

rank(rev(xs))
rank(sample(LETTERS))

Attention au traitement des "ex-æquo" ; la commande

rank(c(1,3,5,5,2))

renvoit le vecteur 1 3 4.5 4.5 2, alors que

rank(c(1,3,5,5,2), ties.method= "first")

retourne 1 3 4 5 2.

Valeurs manquantes

Les valeurs manquantes sont représentées par NA, quel que soit le type des variables.

x <- 1:10
sum(x)

Certaines fonctions exigent de préciser ce que l'on doit faire en présence de valeurs manquantes.

x[c(1,3)] <- NA  ## ou is.na(x[c(1,3)]) <- TRUE
sum(x)
sum(x, na.rm=TRUE)
is.na(x)
sum(is.na(x))

D'autres non :

summary(x)

Les matrices

Une matrice est une structure composée de vecteurs de même taille et de même type.

x1 <- sample(1:100, 10)
x2 <- sample(1:100, 10)
x <- cbind(x1, x2)
x[,"x1"]

On peut créer explicitement des matrices comme suit :

x <- matrix(sample(1:100, 10), ncol=2)
colnames(x) <- paste("x", 1:2, sep="")
dim(x)
ncol(x)  ## ou nrow(x)
x[,1]    ## id. à x[,"x1"]
x[1:3,]
head(x)
str(x)

Opérations sur les matrices

L'adressage des éléments d'une matrice, ou l'assignation de valeurs, suit le même principe que pour les vecteurs, sauf que l'on a deux indices (ligne et colonne).

x <- matrix(1:10, nrow=2, byrow=TRUE)
x[1,2]

Lorsque l'on souhaite effectuer le même type d'opération sur les colonnes (ou les lignes) d'une matrice, au lieu d'une boucle

for (j in 1:ncol(x))
  print(mean(x[,j]))

on peut utiliser la fonction apply :

apply(x, 2, mean)  ## ou colMeans(x), plus efficace

L'argument ... (special variable length argument) permet de passer des arguments supplémentaires lors de l'appel à une fonction, p.ex. :

apply(x, 2, mean, na.rm=TRUE)

Les facteurs

Un facteur permet d'effectuer une partition des unités statistiques : variable catégorielle à valeurs discrètes (niveaux, modalités, strates).

La commande

factor(x, levels=, labels=, ordered=)

convertit x en une variable R de type factor.

Comme dans le cas des matrices, il existe des fonctions spéciales qui permettent d'opérer selon les niveaux d'un facteur :

y <- sample(1:100, 20)
x <- factor(rep(c("a","b"), 10))
tapply(y, x, mean)
by(y, x, mean)

On peut tabuler un facteur, comme n'importe quelle collection d'unités discrètes :

table(x) 

Accéder aux informations d'un facteur

Un facteur reste un vecteur :

x <- factor(rep(c("t1","t2","t3"), each=10), 
            ordered=TRUE)
x[c(1,3)]

avec des attributs supplémentaires qui peuvent être interrogés à l'aide de fonctions spécifiques :

levels(x)      ## unique(x)
nlevels(x)     ## length(unique(x))
is.factor(x)
summary(x)     ## id. à table(x)

R ne fait pas de distinction lexicale entre niveaux (variable ordinale) et modalités (variable nominale), mais l'argument ordered= permet d'imposer une relation d'ordre sur les modalités du facteur.

Créer des facteurs

Au lieu de factor, on peut utiliser gl ou utiliser/créer un vecteur et le convertir en factor.

x <- rep(letters[1:3], each=10)
x <- as.factor(x)
## ou de manière équivalente
gl(3, 10, labels=letters[1:3])

Les commandes interaction et expand.grid sont également utiles pour générer des combinaisons de niveaux de facteurs ("traitements").

Attention à ne pas confondre labels= et levels= :

factor(rep(1:2, 2), labels=2:1)
factor(rep(1:2, 2), labels=2:1, levels=2:1)

Modifier un facteur

On peut vouloir

• changer l'ordre des niveaux (p.ex., catégorie de référence pour une régression),
• supprimer un niveau ou regrouper des niveaux (p.ex., équilibrer les effectifs entre niveaux).

Exemple :

x <- gl(4, 5, 100)
relevel(x, ref=4)
levels(x) <- 4:1
levels(x)[1:2] <- "3.5"

Cela ne change en aucun cas le rang des éléments de x, seulement la modalité qui leur est associée.

Les erreurs fréquentes

En interne, les facteurs sont stockés sous forme numérique mais on leur associe des étiquettes (labels) :

as.numeric(x)

Les "niveaux" d'un facteur sont classés par ordre lexicographique, et la conversion d'un facteur en nombre peut être délicate parfois.

as.numeric(f <- factor(1:2, levels=2:1)[1])
as.numeric(levels(f))[f]

Il paut penser à garder en correspondance les niveaux observés et observables lorsque l'on souhaite modifier un facteur :

a <- gl(4, 5, 100)
b <- a[a != 1]
all(levels(a) == levels(b))
d <- droplevels(a[a != 1])

Il est recommendé de toujours vérifier que R code bien les données telles qu'on se les représente !

Structure de données avancée

Le data.frame est l'objet de base pour gérer des données de type variable (numeric, factor, etc.). Un data.frame est une collection de vecteurs de même taille, mais de type possiblement variable.

score <- c(1,5,3,1,5,2,2)
x1 <- c("M","M","F","F","F","F","M")
x2 <- c(92,96,97,93,101,100,97)
x3 <- c("C","A","C","C","C","B","B")
a <- data.frame(score, gender=x1, 
                IQ=x2, SES=x3)
dim(a)
summary(a)
str(a)
as.matrix(a)

Adresser les éléments d'un data.frame

Même principe d'adressage que pour de simple vecteurs ou des matrices, mais on peut accéder aux colonnes nommées avec l'opérateur $ :

a[,1:2]
a[1:3,c(2,4)]
a$score[1:3]
a$score[a$gender == "M"]
with(a, score[gender == "M"])
## ou de manière équivalente
subset(a, gender == "M")
a$score[a$SES %in% c("A","C")]

Les commandes with ou subset, comme beaucoup de commandes R, acceptent comme argument un data.frame ce qui simplifie l'accès aux variables contenues dans celui-ci. On peut également utiliser attach pour créer une copie dans l'espace de travail, mais ce n'est pas conseillé.

Propriétés d'un data.frame

Certaines commandes R traitent les data.frame de manière spéciale. Comme pour les matrices, on peut utiliser :

summary(a)
str(a)
dim(a)
colnames(a)

On peut ajouter des colonnes ou des lignes à un data.frame, comme dans le cas des matrices, p.ex.

cbind(a, 1:7)

mais

rbind(a, c(3, "M", 100, "C"))

ne produira pas le résultat escompté. Pourquoi ?

Opérations sur un data.frame

À ce stade, gender n'est pas considéré comme un facteur :

is.factor(gender)

mais il a été converti en facteur lors de la création du data.frame (cf. getOption("stringsAsFactors")):

is.factor(a$gender)

On peut convertir ou recoder les variables contenues dans a en facteurs comme suit :

within(a, SES <- factor(SES, labels=c("low","mid","high")))

Plus généralement,

within(data, expr, ...)

permet d'opérer sur les variables d'un data.frame et de rendre ces modifications permanentes, contrairement à with qui ne sert qu'à accéder aux composantes du data.frame (principe d'accès en "lecture seule").

Importer des données externes

Les fonctions read.table et read.csv permettent d'importer des données tabulées ou csv (séparateur de champs virgule ou point-virgule).

Exemple :

blood <- read.table("blood.txt", header=TRUE)
cereal <- read.csv("cereal.csv")

Bien vérifier le chemin d'accès au fichier (ou utiliser set.wd), la présence d'une ligne d'en-tête (header), le type de séparateur de champs (sep), ainsi que le codage des valeurs manquantes (na.strings).

Un data.frame est créé automatiquement lorsque l'on importe des données externes, p.ex. un fichier csv.

Ce qu'il faut retenir

• L'élément de base est le vecteur. Un vecteur a une taille définie et est constitué d'éléments du même type (nombre, caractères, etc.).
• Des structures plus complexes, comme les matrix ou les data.frame, sont constituées de vecteurs.
• On peut effectuer des opérations élémentaires d'indexation et de sélection sur ce type de structures de données.
• Les facteurs jouent un rôle particuliers en statistiques, et ils sont représentées sous forme d'une association nombre/étiquette sous R.
• L'importation de données externes, souvent sous forme de fichiers plats, résulte en un data.frame.

Analyse exploratoire des données

Objectifs

Au travers de l'analyse exploratoire des données, on cherche essentiellement à résumer la distribution de chaque variable (approche univariée) ainsi que les relations entre les variables (approche bivariée essentiellement), dont les caractéristiques pourraient suggérer un recodage ou une transformation des mesures (J. W. Tukey, 1977).

Plutôt que de modéliser directement les données, on s'attachera donc dans un premier temps à les décrire à l'aide de résumés numériques et graphiques. L'idée est de caractériser la forme d'une distribution et d'identifier les éventuelles valeurs influentes.

On profitera de cette approche pour présenter les principales fonctionnalités graphiques de R, en particulier l'interface lattice (Murrell, 2005; Sarkar, 2008).

Un jeu de données d'exemple

The low birth weight study

Il s'agit d'une étude prospective visant à identifier les facteurs de risque associés à la naissance de bébés dont le poids est inférieur à la norme (2,5 kg). Les données proviennent de 189 femmes, dont 59 ont accouché d'un enfant en sous-poids. Parmi les variables d'intérêt figurent l'âge de la mère, le poids de la mère lors des dernières menstruations, l'ethnicité de la mère et le nombre de visites médicales durant le premier trimestre de grossesse.

C'est une des études qui sert de base aux traitements statistiques présentés dans Hosmer & Lemeshow (1989).

Elle est disponible sous R dans le package MASS :

data(birthwt, package="MASS")

Traitement préalable

Les données fournies dans R nécessitent quelques recodages:

str(birthwt)
summary(birthwt)
birthwt <- within(birthwt, {
  low <- factor(low, labels=c("No","Yes"))
  race <- factor(race, labels=c("White","Black","Other"))
  smoke <- factor(smoke, labels=c("No","Yes"))
  ui <- factor(ui, labels=c("No","Yes"))
  ht <- factor(ht, labels=c("No","Yes"))
})

Il est intéressant d'ajouter les unités de mesure, pour s'en souvenir lorsqu'on en a besoin (le poids de la mère est en pounds, celui des bébés en kg !).

library(Hmisc)
birthwt <- within(birthwt, {
  units(age) <- "years"
  units(lwt) <- "pounds"
})

Valeurs extrêmes, atypiques ou outliers

Il n'y a pas vraiment de consensus sur la dénomination correcte des valeurs qui "ne ressemblent pas" à la majorité des valeurs observées. Une valeur extrême ou atypique est souvent assimilée à une valeur qui est située à plus de \(1.5\times\text{IQR}\) des quartiles supérieurs et inférieurs. Un outlier est une valeur susceptible d'influencer les résultats obtenus par un modèle statistique.

idx <- sapply(birthwt, is.numeric)
bwt <- apply(birthwt[,idx], 2, scale)
boxplot(bwt)

On peut aussi chercher des "patterns" multivariés (p.ex., des individus avec des valeurs sytématiquement élevées ou basses).

parallel(bwt, groups=birthwt$low, horiz=FALSE)
idx <- apply(bwt, 2, filter.perc, cutoff=c(.01,.99), 
             collate=TRUE)
my.col <- as.numeric(1:nrow(bwt) %in% unique(unlist(idx)))+1
splom(~ bwt, pch=19, col=my.col, alpha=.5, cex=.6)

Résumé de la structure de données

Un aperçu synthétique des données, stratifié par groupe de poids des bébés, peut être obtenu comme suit :

summary(low ~ ., data=birthwt[,-10], method="reverse")
library(latticeExtra)
marginal.plot(birthwt, data=birthwt, groups=low)    

Synthèse numérique et transformation

On peut aussi recoder certaines des variables discrètes (ftv et ptl) en variables binaires, pour la description ou pour la modélisation :

bwt.df <- transform(birthwt[,-10], 
                    ftv=factor(ftv>0, lab=c("No","Yes")),
                    ptl=factor(ptl>0, lab=c("None","1+")))
summary(low ~ ., data=bwt.df, method="reverse")

Il n'est pas nécessaire de stratifier pour résumer les données, mais dans le cas présent il est intéressant de vérifier la distribution des variables pour les deux groupes de bébés. En ajoutant l'option overall=TRUE dans l'expression ci-dessus, on obtient également le résumé numérique sur l'ensemble de l'échantillon.

Concernant les unités de mesure, on peut convertir "à la volée" lors de l'appel à des fonctions comme mean ou summary.formula (ci-dessus) :

mean(birthwt$lwt/2.2)  ## poids de la mère en kg
summary(lwt/2.2 ~ low + race, data=birthwt)

Synthèse numérique et transformation (2)

D'autres types de résumés numériques peuvent être produits avec summary.formula, en particulier des tableaux croisés ou des descriptions stratifiées. Pour plus d'informations, consulter l'excellent guide Hmisc : http://biostat.mc.vanderbilt.edu/wiki/Main/Hmisc.

Qu'est-ce qu'une distribution?

Considérons l'âge des mères, variable numérique souvent assimilée à une variable continue mais qui ici prend plutôt des valeurs discrètes. La plupart des valeurs prises par la variable age semble se concentrer autour de la tranche 20-25 ans.

stripplot(~ age, data=birthwt, jitter.data=TRUE, 
          amount=.3, aspect=.3, cex=.6)

Résumé numérique

Mesures de tendance centrale (moyenne, médiane) associées à des mesures de dispersion relative (écart-type, IQR).

## Tukey's five-point summary
summary(birthwt$age)
quantile(birthwt$age, probs=c(.1, .25, .5, .75, .9))
desc <- function(x, dig=2)
  round(c(ety=sd(x), iqr=IQR(x), "max-min"=diff(range(x))), 
        digits=dig)
desc(birthwt$age)

Autres possibilités : Estimateurs robustes

• Moyennée tronquée : mean(birthwt$age, trim=.025)
• Déviation médiane absolue : mad(birthwt$age)

Fonction de répartition

On peut représenter la distribution des effectifs cumulés à l'aide d'un diagramme quantile-quantile.

qqmath(~ age, data=birthwt, dist=qunif)

L'idée de base consistant à "trier" les données permet en outre d'identifier facilement n'importe quel quantile. Un algorithme pour calculer la valeur médiane est indiqué ci-dessous.

med <- function(x) {
  odd.even <- length(x) %% 2
  if (odd.even == 0) 
    (sort(x)[length(x)/2] + 
      sort(x)[1+length(x)/2]) / 2 
  else 
    sort(x)[ceiling(length(x)/2)]
}

On pourrait également comparer cette distribution empirique à une distribution théorique, par exemple une distribution gaussienne.

Approche graphique

On peut visualiser les cinq indicateurs renvoyés par la fonction summary à l'aide d'une "boite à moustaches" (Wickham & Stryjewski, sub.).

bwplot(~ age, data=birthwt)

Dans la figure ci-dessous, la moyenne est indiquée par un point, l'étendue des données est représentée par les "moustaches", et 50 % des observations sont comprises entre les 1er et 3ème quartiles figurés par la "boîte".

Histogramme

Un histogramme permet de représenter la répartition des valeurs d'une variable continue en classes. Ici, on voit que la distribution du poids des mères est asymétrique (asymétrie > 0, indépendant de l'unité de mesure).

histogram(~ lwt, data=birthwt, type="count")
library(e1071)
skewness(birthwt$lwt)

Histogramme (variations)

Le paramètre breaks permet de changer le nombre d'intervalles construits. Par défaut, la méthode utilisée est la méthode de Sturges (Sturges, 1926). L'ajout d'une loi de densité gaussienne dont les paramètres sont estimés* à partir de l'échantillon est également possible.

Densité empirique

Pour pallier à l'arbitraire du choix du nombre d'intervalles, on peut préférer représenter la fonction de densité empirique (Silverman, 1986; Venables & Ripley, 2002). Il reste toutefois à définir la largeur de la fenêtre de lissage associée à la fonction noyau de lissage (voir help(bw.nrd0)).

densityplot(~ lwt, data=birthwt)

Résumé numérique bivarié

Dans le cas de deux variables continues, on peut s'intéresser à la distribution jointe ou quantifier une certaine mesure de leur covariation, p.ex.

with(birthwt, cor(lwt, bwt))
with(birthwt, cor(lwt, bwt, method="spearman"))

Lorsque l'on croise une variable numérique avec une variable catégorielle, on peut produire des résumés numériques séparés pour chaque niveau de la variable de classification.

with(birthwt, by(age, low, summary))
with(birthwt, tapply(lwt, race, 
                     function(x) c(mean(x), sd(x))))

Des fonctionnalités étendues sont disponibles dans Hmisc. La commande ci-dessus se résume à

summary(race ~ lwt, data=birthwt, method="reverse")

ou encore, si l'on souhaite afficher moyenne et écart-type,

print(summary(race ~ lwt, birthwt, method="reverse"), 
      prmsd=TRUE)

Diagramme de dispersion

Pour deux variables numériques, un diagramme de dispersion permet de résumer rapidement la distribution conjointe, ici dans les mêmes unités* :

xyplot(lwt/2.2 ~ bwt/1000, data=birthwt, type=c("p","r"))

La concentration des points autour de la droite de régression est un bon indicateur de la "consistance" de l'hypothèse de linéarité entre x et y.

Scatterplot smoother

Supposer la linéarité n'est pas toujours la meilleure des options, et peut masquer certaines distorsions locales dans la relation entre les deux variables. Pour cette raison, on peut préférer utiliser une courbe plus flexible de type loess (Cleveland, 1979).

Pour cela, on remplacera avantageusement type=c("p","r") (points + droite de régression) par type=c("p","smooth") (points + lowess).

Distributions conditionnelles

Les diagrammes de type boîte à moustaches et densité empirique peuvent être conditionnés sur un facteur, ici l'âge de la mère en fonction de l'indicateur de poids à la naissance :

bwplot(low ~ age, data=birthwt, panel=Hmisc::panel.bpplot)

Distributions conditionnelles (2)

L'avantage des densités non-paramétriques est que l'on peut visualiser directement les variations dans l'allure des distributions conditionnelles (tendance centrale et forme de la distribution).

densityplot(~ age, data=birthwt, groups=low, 
            auto.key=list(columns=2), plot.points=FALSE)

Mesures aggrégées

Plutôt que des diagrammes en barres (barchart), on préférera les diagrammes de type dotplot (Cleveland, 1985) pour représenter les données aggrégées, telles que des moyennes ou des proportions*.

age.by.race <- aggregate(age ~ race, data=birthwt, FUN=mean)
dotplot(race ~ age, data=age.by.race)

Croiser plus de deux variables

Les graphiques en trellis (Becker, Cleveland, & Shyu, 1996; Cleveland, 1993) offrent une structure de graphique simple et efficace pour représenter des données multidimensionnelles. En particulier, ils introduisent la notion de facettes pour représenter des distributions conditionnelles.

La notation utilisée est la notation par formule :

y ~ x | a      ## y en fonction de x condit. à a
y ~ x | a + b  ## y en fonction de x condit. à a et b

Les variables continues peuvent être "catégorisées" à l'aide de shingles. En transformant les variables numériques en facteurs, il devient possible de représenter un plus grand nombre de croisement de variables, tout en tenant compte de la nature "continue" des données.

Quelques exemples :

xyplot(bwt/1000 ~ lwt/2.2 | smoke, data=birthwt, 
       groups=ptl>0)
bwplot(age ~ low | smoke + race, data=birthwt)

Illustrations

xyplot(bwt ~ lwt/2.2 | race, data=birthwt, layout=c(3,1),
       cex=sqrt(birthwt$ftv+.5), col=birthwt$smoke, 
       panel=function(...) {
         panel.xyplot(...)
         panel.abline(h=2500, lty=2)})

Illustrations (2)

Age <- equal.count(birthwt$age)
ftvc <- factor(birthwt$ftv>1, 
               labels=c("0 visite","1+ visite"))
xyplot(bwt/1000 ~ lwt/2.2 | Age + ftvc, data=birthwt, 
       groups=low, pch="+")

Ce qu'il faut retenir

• Il est important de vérifier le codage des variables, et de recoder en fonction des besoins de l'analyse ou de la visualisation.
• On caractérise d'abord les distributions univariées (toutes !) avant de passer aux visualisations/modèles multivariés. Cela permet de détecter les éventuelles valeurs aberrantes, la mauvaise représentation de certaines modalités d'une variable catégorielle, ou l'existence d'asymétrie dans les distributions.
• L'interface lattice utilise la même notation par formule que les fonctions R pour la modélisation statistique.
• La visualisation des données catégorielles n'a pas vraiment été abordée, mais de nombreux outils sont disponibles dans les packages vcd et vcdExtra (Friendly, 2011; Meyer, Zeilis, & Hornik, 2006).

Mesures et tests d'association

Objectifs

• Dans ce cours, on s'intéressera aux tests d'inférence simples visant à comparer des indices de tendance centrale calculés sur des données continues, des proportions, des mesures de concordance, estimés à partir de données observées à la suite d'un procédé d'échantillonnage.
• On considèrera deux échantillons ou plus, indépendants ou non.
• L'objectif est de démontrer le principe général de la démarche du test d'hypothèse à partir de quelques statistiques de test rencontrées fréquemment en recherche biomédicale.

Ce cours s'appuie en partie sur l'ouvrage de Everitt & Rabe-Hesketh (2001).

Hypothèse nulle, risque d'erreur en inférence fréquentiste

De manière générale, lorsque l'on s'intéresse à un effet particulier, l'idée est de postuler l'absence d'effet (hypothèse nulle, \(H_0\)) et de chercher à vérifier si les données observées sont compatibles ou non avec cette hypothèse. C'est le principe même de la démarche hypothético-déductive.

Pour réaliser un tel test, il est nécessaire d'avoir (ou de construire) une statistique de test, dont la distribution d'échantillonnage est connue (ou peut être approximée) sous \(H_0\), qui nous permettra de répondre à la question suivante : en supposant qu'il n'existe pas d'effet dans la population, quelle est la probabilité d'observer une statistique de test au moins aussi extrême que celle estimée à partir de l'échantillon choisi aléatoirement dans cette population ? Si cette probabilité se révèle "suffisamment petite", on concluera qu'il est vraisemblablement peu probable que le résultat observé soit dû simplement au hasard de l'échantillonnage.

Risque de première et deuxième espèce

Supposons que l'on est à prendre une décision concernant une hypothèse nulle. Le fameux "risque alpha" (type I) est le risque de conclure à tort à l'existence d'un effet alors qu'en réalité ce dernier n'existe pas.

À ce risque est typiquement associé, de manière asymétrique, le risque (type II) de ne pas rejeter \(H_0\) lorsque celle-ci est en réalité fausse ; le complémentaire de ce risque est appelée la puissance.

Ces risques sont effectivement asymétriques :

• dans un essai thérapeutique, si l'on doit décider si on un nouveau traitement est meilleur que le traitement courant, on cherche à minimiser le risque d'une mauvaise décision (\(\alpha=0.05\)). ;
• si à la fin de l'essai, on ne met en évidence aucune différences significatives, cela ne signifie pas que les traitements sont équivalents : il existe un risque \(\beta\) qu'il existe une réelle différence entre les deux.

En résumé

Comparer deux moyennes

Le test de Student est un test paramétrique permet de tester si deux moyennes de groupe (indépendants ou appariés) peuvent être considérées comme significativement différentes en considérant un seuil d'erreur \(\alpha\).

Exemples d'application : comparer un dosage biologique entre deux groupes de patients, comparer des mesures physiologiques avant et après traitement chez les mêmes patients.

Conditions d'application : normalité des distributions parentes, homogénéité des variances, indépendance (dans le cas du test t pour échantillons indépendants).

Application

Poids à la naissance : mères fumeur vs. non-fumeur. (Hosmer & Lemeshow, 1989)

data(birthwt, package=MASS)
birthwt$smoke <- factor(birthwt$smoke, 
                        labels=c("No","Yes"))
t.test(bwt ~ smoke, data=birthwt, var.equal=TRUE)

On conclut à l'existence d'une différence entre les poids moyens des enfants des deux groupes de mère, \(p<0.01\) (\(\Delta=+284\) g \([73;495]\), en faveur des non-fumeurs).

Vérification des conditions d'application

On doit essentiellement vérifier l'homogénéité des variances et la distribution des résidus.

bwplot(smoke ~ bwt, data=birthwt)
qqmath(~ bwt, data=birthwt, group=smoke)

Cas des données non indépendantes

Dans le cas où les échantillons ne sont pas indépendants, p.ex. sujets mesurés à deux reprises ou données appariées, le même type de test peut être utilisé, et la différence de moyennes est comparée à 0.

Le plus souvent ce sont les mêmes unités statistiques qui servent de support à la comparaison, mais n'importe quelle forme d'appariement, pourvu qu'elle fasse sens, peut justifier le recours à un test apparié.

Application

Effet de somnifères sur le temps de sommeil. (Student, 1908, p. 20)

t.test(extra ~ group, data=sleep, 
       var.equal=TRUE, paired=TRUE)

On constate une diminution de la durée de sommeil lorsque le médicament 1 est administré (\(\Delta=-1.6\) h \([-2.5;-0.7]\), \(p<0.01\)).

Analyse graphique

On peut utiliser un diagramme de dispersion, de Tukey (tmd) ou de Bland-Altman (Altman & Bland, 1983).

xyplot(extra[group==2] ~ extra[group==1], 
       data=sleep, type=c("p","g"))

Alternatives non-paramétriques

Dans le test t, on fait une hypothèse sur la nature de la distribution de la réponse mesurée dans la population parente. Si l'on relaxe cette hypothèse et que l'on exige simplement que les échantillons proviennent de populations ayant des distribution à peu près comparables en terme de forme, alors on peut utiliser des tests dits non-paramétriques. Dans le cadre de la comparaison de la tendance centrale de deux échantillons, l'alternative au test t est le test de Wilcoxon.

Ce type de test est souvent utilisé dans le cas des petits effectifs lorsque les données disponibles ne permettent pas réellement de vérifier la normalité des distributions parentes, comme dans le test t. On retiendra toutefois que le test de Wilcoxon a une puissance relative de 80 % par rapport au test de Student. Ce type de test présente en outre l'avantage d'être moins sensible à la présence de valeurs extrêmes.

Application

Hamilton depression scale scores. (Hollander & Wolfe, 1999, p. 29)

occ1 <- c(1.83,0.50,1.62,2.48,1.68,1.88,1.55,3.06,1.30)
occ2 <- c(0.878,0.647,0.598,2.05,1.06,1.29,1.06,3.14,1.29)
wilcox.test(occ1, occ2, paired=TRUE)

Le résultat du test suggère une amélioration de l'état (diminution du score Hamilton de mean(occ2-occ1)=-0.43 points).

Comparer plus de deux moyennes

L'ANOVA constitue une extension naturelle au cas où plus de deux moyennes de groupe sont à comparer. Attention, avec \(k\) échantillons, l'hypothèse nulle se lit : \[ H_0:\ \mu_1=\mu_2=\ldots=\mu_k, \] alors que l'alternative est l'existence d'au moins une paire de moyennes qui diffèrent (négation logique de \(H_0\)). Si l'on exprime chaque observation comme une déviation par rapport à sa propre moyenne de groupe, \(y_{ij}=\bar y_i+\varepsilon_{ij}\), on voit que la variabilité totale peut se décomposer comme suit : \[\underbrace{(y_{ij}-\bar y)}_{\text{totale}}=\underbrace{(\bar y_{i\phantom{j}}\hskip-.5ex-\bar y)}_{\text{groupe}} + \underbrace{(y_{ij}-\bar y_i)}_{\text{résiduelle}}.\]

Conditions d'application : normalité des distributions parentes pour chaque groupe, homogénéité des variances, indépendance des observations.

Illustration

Supposons trois groupes indépendants. Voici quelques scénarios imaginaires concernant la distribution des scores individuels :

Tableau d'ANOVA

On peut résumer les sources de variations observées au niveau de la variable réponse comme suit :

On rappelle qu'une variance se calcule comme la somme des écarts quadratiques à la moyenne, i.e. \(\sum_j (y_j-\bar y)^2/(n-1)\).

La \(SC_b\) représente la variance des moyennes de groupe (fluctuation autour de la moyenne générale), et la \(SC_w\) ("résiduelle") représente la moyenne des variances des mesures individuelles autour de la moyenne de groupe (dans les deux cas, les moyennes sont pondérées).

Application

Polymorphisme et gène du récepteur estrogène. (Dupont, 2009)

library(foreign)
polymsm <- read.dta("polymorphism.dta")[,-1]
with(polymsm, tapply(age, genotype, mean))
aov.res <- aov(age ~ genotype, data=polymsm)
summary(aov.res)
model.tables(aov.res)  
plot.design(age ~ genotype, data=polymsm)

Le test F indique que l'on peut rejeter l'hypothèse nulle d'absence de différence d'âge entre les trois génotypes (Df=3-1=2). Il est possible d'utiliser le rapport 2316/(2316+8246)=0.219 (\(\eta^2\)) pour quantifier la variance expliquée. On pourrait poser la question de savoir quelles paires de moyennes diffèrent significativement (comparaisons multiples post-hoc).

Distributions conditionnelles

Ci-dessous un diagramme en forme de boites à moustaches pour les trois génotypes. Le trait horizontal indique la moyenne globale.

bwplot(age ~ genotype, data=polymsm)

Diagnostic du modèle

Essentiellement, il faut vérifier la distribution des résidus

• à l'aide d'un QQ plot (normalité des résidus),
• selon les valeurs prédites (homoscédasticité).

On peut utiliser plot(aov.res) ou rfs (ci-dessous).

Diagramme en barres pour publication

Les diagrammes en barres ne sont pas très économes en termes de "data-ink ratio". On peut leur préferer des diagrammes en points (Cleveland, 1993), voir dotplot.

Quelques remarques

• Dans le cas où on compare deux groupes (indépendants), un test t ou un test F donneront le même résultat (\(t^2 = F\)) :

  summary(aov(bwt ~ smoke, data=birthwt))
  t.test(bwt ~ smoke, data=birthwt, var.equal=TRUE)

• Pour quoi ne pas faire une simple série de tests t au lieu de l'ANOVA ? Comparer les \(m=\complement_2^{k(k-1)}\) paires de moyennes entraîne une inflation du risque de première espèce puisque le risque global vaut alors \(1-(1-\alpha)^m\) (en supposant les tests indépendants), soit 40 % de chances de conclure à tort à une différence lorsqu'en réalité il n'y en a pas (\(H_0\) vraie) sur 5 groupes.

• On peut voir l'ANOVA comme une simple comparaison de modèles, i.e. modèle complet (avec prédicteur) vs. modèle nul (résiduelle) :

  aov0 <- aov(age ~ 1, data=polymsm)
  aov1 <- aov(age ~ genotype, data=polymsm)     
  anova(aov0, aov1)

Alternative non-paramétrique

La comparaison de k groupes peut se faire à l'aide d'un test non-paramétrique, reposant sur des sommes de rangs. L'ANOVA de Kruskal-Wallis, qui est une extension du test de Wilcoxon-Mann-Whitney, permet de tester si la distribution de la variable réponse peut être considérée comme identique entre les k groupes. On suppose toujours les groupes indépendants.

kruskal.test(age ~ genotype, data=polymsm)

Mesures de corrélation

Le coefficient de corrélation de Bravais-Pearson permet d'évaluer le degré d'association linéaire entre deux variables continues. Les variables jouent un rôle symétrique. En pratique, il est calculé à partir de la covariance entre les deux variables, et normalisé pour prendre des valeurs entre -1 (corrélation négative parfaite) et +1 (corrélation positive parfaite).

Il est possible de tester si un coefficient de corrélation est nul (\(H_0:\,\rho=0\)) ou égal à une certaine valeur.

Illustration

Dans chacun des cas ci-dessous, le coefficient de corrélation de Bravais-Pearson vaut 0.816 ! D'où l'importance de prendre en considération de possibles déviations par rapport à l'hypothèse de linéarité (b), la présence de valeurs atypiques ou extrêmes (c et d).

Application

Health survey of paint sprayers. (Everitt & Rabe-Hesketh, 2001, p. 158)

paint <- read.table("PAINT.DAT", header=TRUE)
with(paint, cor.test(HAEMO, PCV))

On conclut à l'existence d'une corrélation significative entre le taux hématocrite et la concentration en hémoglobine (\(p<0.001\)).

Approche graphique

Un diagramme de dispersion est nécessaire pour vérifier la présence éventuelle de valeurs aberrantes ou extrêmes. Il est possible de recalculer la corrélation sans les points définissant l'enveloppe du nuage de points ("convex hull").

xyplot(PCV ~ HAEMO, data=paint, type=c("p","g","r"))

Alternative non-paramétrique

Le coefficient de corrélation de Spearman est calculé à partir des rangs des observations et par conséquent permet d'estimer une mesure de la monotonicité de la relation entre deux variables continues. Il est évidemment moins sensible aux valeurs extrêmes. Le calcul de la statistique de test est comparable au cas de la corrélation linéaire de Pearson, excepté que l'on remplace les valeurs observées par leurs rangs.

Le coefficient de Spearman peut également être utilisé sur des variables ordinales, en particulier lorsque l'on souhaite s'affranchir des hypothèses sur la distribution continue des réponses mesurées.

library(foreign)
anorex <- read.spss("anorectic.sav", to.data.frame=TRUE)
with(subset(anorex, time==1), 
     cor.test(binge, purge, method="spearman"))

Test d'association pour tableau à deux entrées

Lorsque les données sont qualitatives, ou catégorielles, un tableau de contingence permet de résumer la distribution des effectifs pour chaque croisement des modalités des variables. On se limitera ici au cas le plus simple avec deux variables. Notons que l'on a besoin de connaître les marges du tableau (totaux-lignes et colonnes) donc on a besoin des effectifs, et pas seulement de données aggrégées comme des proportions.

La question de l'association entre ces deux variables revient le plus souvent à formuler une hypothèse nulle d'indépendence entre les deux variables étudiées. La notion d'association entre variables catégorielles peut être envisagée comme une certaine mesure de corrélation telle que celle discutée précédemment.

Une discussion approfondie de ces mesures d'association est disponible dans Bishop, Fienberg, & Holland (2007).

Illustration

Dans le cas général, un tableau des effectifs (ou des fréquences relatives) peut être construit comme indiqué ci-dessous. Les \(n_{i\cdot}\) (\(i=1,\dots,I\)) et les \(n_{\cdot j}\) (\(j=1,\dots,J\)) sont les distributions marginales des variables A et B, respectivement.

Les modalités de A et B peuvent être nominales et/ou ordinales, les totaux lignes/colonnes peuvent être fixés ou non, il peut y avoir des cellules vides, les variables A et B peuvent être liées (pre/post, inter-juges), etc.

Extensions possibles

Plusieurs cas de figure sont envisageables : tableau 2x2 (e.g., exposition/maladie), marge(s) fixe(s), tableau I lignes x J colonnes stratifié sur deux groupes, tableau croisé avec variables en commun, tableau I lignes x J colonnes répétés dans le temps.

Écart à l'indépendance

Rappelons que, dans le cas général d'un tableau I x J, la statistique de Pearson consiste à comparer les données observées aux données attendues sous une hypothèse d'indépendance ("produit des marges"), soit, à une constante près, la somme des écarts quadratiques entre les effectifs observée et les effectifs théoriques ("résidus"), normalisés par ces derniers, soit \[ \Phi^2=\sum_{i=1}^I\sum_{j=1}^J\frac{(p_{ij}-p_{i\cdot}p_{\cdot j})^2}{p_{i\cdot}p_{\cdot j}}. \]

Le "chi-deux" de Pearson est simplement \(N\Phi^2\), et cette statistique de test suit une loi \(\chi^2\) à \((I-1)(J-1)\) degrés de liberté. Dans le cas 2x2, on a

\[ \Phi^2=\sum_{i=1}^2\sum_{j=1}^2\frac{(p_{ij}-p_{i\cdot}p_{\cdot j})^2}{p_{i\cdot}p_{\cdot j}}=\frac{(p_{11}p_{22}-p_{21}p_{12})^2}{p_{1\cdot}p_{2\cdot}p_{\cdot 1}p_{\cdot 2}},\]

Approche générale

On peut regrouper les mesures d'association pour les tableaux 2 x 2 en deux classes, celles qui reposent sur une fonction du produit croisé des cellules et celles qui se basent sur le coefficient de corrélation de Pearson.

Dans le premier cas, on retrouve l'odds-ratio, \[\alpha=\frac{p_{11}p_{22}}{p_{12}p_{21}},\] encore exprimé sous la forme \(\alpha=\frac{p_{11}/p_{12}}{p_{21}/p_{22}}\) lorsque l'on considère les totaux-lignes fixés (e.g., exposition) : \(p_{11}/p_{12}\) est alors la probabilité d'être malade (vs. non-malade) sachant que l'on a été exposé (vs. non-exposé).

Dans le second cas, si l'on assimile les modalités des variables à deux scores 0/1, alors la corrélation se définit naturellement comme

\[ \rho=\frac{p_{11}p_{22}-p_{21}p_{12}}{\sqrt{p_{1\cdot}p_{2\cdot}p_{\cdot 1}p_{\cdot 2}}}=\frac{p_{22}-p_{2\cdot}p_{\cdot 2}}{\sqrt{p_{1\cdot}p_{2\cdot}p_{\cdot 1}p_{\cdot 2}}},\quad \text{d'où}\, r^2=\chi^2/N. \]

Quelques remarques

• Il y a de nombreux "tests de chi-deux". Ne pas confondre la loi que suit une statistique de test avec l'usage consacré du \(\chi^2\) pour les tableaux de contingence !
• Les tests du \(\chi^2\) d'indépendance dans les tableaux I x J permettent de répondre globalement à la question : y'a-t-il ou non association entre les variables étudiées ? La nature précise de cette association doit être étudiée plus en détails (résidus, décomposition du \(\chi^2\), direction et taille d'effet, etc.).
• Outre le fait qu'ils nécessitent généralement de grands échantillons, ce type de tests ne reposent que sur les sommes marginales, sans considération de la nature des variables (nominale, ordinale). Des tests plus adaptés (et économes en termes de degrés de liberté) existent dans le cas où l'une des variables ou les deux sont ordonnées.
• Bien que l'on ne les présente pas, il existe de nombreuses fonctions graphiques pour la visualisation des données catégorielles dans les packages vcd et vcdExtra (Friendly, 2000).

Application 1

Swedish study on aspirin use and myocardial infarction. (Agresti, 2002, p. 72)

aspirin <- matrix(c(28,18,656,658), nrow=2)
dimnames(aspirin) <- list(c("Placebo","Aspirin"), 
                          c("Yes","No"))
library(vcd)
asp.or <- oddsratio(aspirin, log=FALSE)
print(list(or=asp.or, conf.int=confint(asp.or)))
summary(oddsratio(aspirin))

L'estimé du "vrai" OR est plutôt imprécis (l'IC à 95% est large et contient la valeur 1), et le test sur le log(OR) n'est pas significatif (\(p=0.071\)). On concluera que rien ne permet d'affirmer que la probabilité de décès dûe à un infarctus diffère selon le facteur d'exposition.

Application 2

Consommation de cafféine et statut marital. (Dalgaard, 2008, p. 84)

coffee <- matrix(c(652,1537,598,242,36,46,38,21,
                   218,327,106,67), nrow=3, byrow=TRUE)
dimnames(coffee) <- list("marital status"=c("Married",
  "Prev.married","Single"), consumption=c("0","1-150",
  "151-300",">300"))
prop.table(coffee, 1)
(chsq <- chisq.test(coffee))
chsq$residuals

On conclut à l'existence d'une association forte entre le statut marital et la consommation de caffeine (\(p<0.001\)). Les célibataires sont beaucoup plus nombreux à ne pas consommer de café en comparaison de ce que l'on attendrait en cas d'indépendance.

Illustration

Comme on l'a dit, des graphiques de type dotplot sont largement préférables à la plupart des alternatives de type diagramme en barres.

dotplot(consumption ~ value, data=melt(coffee), 
        groups=marital.status)

Application 3

Données fictives pré/post traitement.

x <- cbind(c(50,10),c(20,40))
dimnames(x) <- list(c("+","-"),
                    c("+","-"))
margin.table(x, 1)
mcnemar.test(x)
mcnemar.test(x)
binom.test(apply(x, 1, sum), .5)

Le résultat non-significatif du test (\(p=0.100\)) suggère que la répartition entre cas positifs et négatifs n'a pas évolué avec le traitement.

Pour plus d'informations, se référer à Campbell (2007).

Application 4

Fisher's Tea Drinker. (Agresti, 2002, p. 92)

TeaTasting <- matrix(c(3,1,1,3), nrow=2,
       dimnames=list(Guess=c("Milk","Tea"),
                     Truth=c("Milk","Tea")))
fisher.test(TeaTasting, alternative="greater")

Rien ne permet d'établir qu'il existe une association entre l'ordre d'ajout du thé et du lait et le choix du répondant sur la base de ces résultats (\(p=0.2429\)).

Concordance

Lorsque deux ou plusieurs évaluations catégorielles indépendantes sont effectuées sur la même unité d'analyse, et que l'on s'intéresse à l'association ou la corrélation entre les séries d'évaluations, on parle de concordance ("agreement" dans la littérature anglo-saxonne).

Dans le cas de données numériques, on utiliserait plutôt le coefficient de corrélation intra-classe (Shrout & Fleiss, 1979).

Application

Diagnoses of carcinoma. (Agresti, 2002, p. 432)

data(pathologist.dat, package="exactLoglinTest")
patho <- xtabs(y ~ A + B, data=pathologist.dat)
sum(diag(patho))/sum(patho)
library(reshape)
patho.expand <- untable(pathologist.dat[,2:3], 
                        pathologist.dat$y)
library(irr)
kappa2(patho.expand)

La concordance entre les deux médecins apparaît relativement bonne puisque la différence avec ce que l'on attendrait sous l'hypothèse d'indépendance est presque de 50 %. On pourrait y associer un IC à 95 # calculé par bootstrap.

Modèle linéaire et applications

Objectifs

• Dans ce cours, on s'intéressera à l'approche de modélisation par régression linéaire, dans un cadre uni- et multivarié, en montrant le lien avec l'analyse de variance et de covariance.
• Le cas des données corrélées (mesures répétées) sera traité séparément.
• On insistera en particulier sur l'estimation ponctuelle et par intervalles des paramètres des modèles et les procédures de diagnostic des modèles.

Lecture conseillée : Vittinghoff, Glidden, Shiboski, & McCulloch (2005)

Un même objectif

Expliquer les variations observées au niveau d'une variable réponse ("dépendante") numérique, \(y\), en fonction de variables prédictrices ("indépendantes"), \(x_j\), pouvant être de nature qualitative ou quantitative.

Illustration

L'idée revient toujours à considérer qu'il existe une part systématique et une part aléatoire (résidus) dans ces variations. Le modèle linéaire permet de formaliser la relation entre \(y\) et les \(x_j\), en séparant ces deux sources afin d'estimer la contribution relative des \(x_j\) dans les fluctuations de \(y\).

\[ \text{réponse} = \text{effet prédicteur} + \hskip-11ex \underbrace{\text{bruit}}_{\text{erreur de mesure, temps d'observation, etc.}}\]

sachant que le modèle théorique relie fonctionnellement la réponse au(x) prédicteur(s) de manière additive : \(\mathbb{E}(y \mid x) = f(x; \beta)\).

Régression linéaire simple

La régression linéaire permet de modéliser la relation linéaire entre une variable réponse continue et un prédicteur d'intérêt.

Contrairement à l'approche corrélationnelle, dans ce cas les deux variables jouent un rôle asymétrique : la variable explicative ou prédictrice (x) est supposée expliquer une partie des variations observées au niveau de la variable réponse y. On peut vouloir quantifier cette part de variance expliquée, ou encore estimer la contribution de x dans les variations de y.

La linéarité de la relation et l'influence des observations sont deux aspects critiques de la validité des résultats.

Alternatives possibles :

• méthodes résistantes (Huber) ou
  robustes (LAD, quantile)
• splines ou restricted cubic splines

Illustration

On simule 10 paires d'observations indépendantes, liées par la relation \(y=5.1+1.8\times x\), à laquelle on rajoute des aléas gaussiens N(0,1).

n <- 10
x <- runif(n, 0, 10)
y <- 5.1 + 1.8 * x + rnorm(n)
summary(lm(y ~ x))

Illustration (2)

Les valeurs prédites \(\tilde y\) sont estimées conditionnellement aux valeurs prises par \(x\). L'hypothèse de normalité (\(\varepsilon_i\sim\mathcal{N}(0,\sigma^2)\)) n'est pas nécessaire pour estimer la pente par MCO, mais seulement pour l'inférence sur les paramètres du modèle.

Application

Health survey of paint sprayers. (Everitt & Rabe-Hesketh, 2001, p. 158)

paint <- read.table("PAINT.DAT", header=TRUE)
xyplot(PCV ~ HAEMO, data=paint, type=c("p","r"))

Estimation des paramètres du modèle de régression

lm.fit <- lm(PCV ~ HAEMO, data=paint)
summary(lm.fit)
confint(lm.fit)
anova(lm.fit)

Le test pour \(H_0: \beta_1=0\) est significatif (\(p<0.001\)). Lorsque la concentration en hémoglobine (HAEMO) varie de une unité, le taux hématocrite varie de 2.3 %.

Diagnostic du modèle

L'essentiel de l'activité de diagnostic du modèle (mauvaise spécification et points influents) repose sur l'analyse des résidus, ceux-ci étant accessibles à partir de la commande resid.

xyplot(resid(lm.fit) ~ HAEMO, data=paint)
xyplot(resid(lm.fit) ~ fitted(lm.fit))

Avec la commande influence, on dispose également d'un ensemble d'indices statistiques permettant d'évaluer la qualité d'ajustement du modèle ; voir aussi help(influence.measures).

influence.measures(lm.fit)

Observations extrêmes et influentes

On peut distinguer différents types de points influents :

• "outlier" univarié (x ou y)
• point présentant un effet levier ("leverage") : valeur extrêmes sur x
• point ayant un large résidu (en valeur absolue) : valeur extrêmes sur y
• point exerçant un effet sur les coefficients de régression (sensibilité)

La fonction resid permet d'obtenir les résidus "simples", \(e_i=y_i-\hat y_i\), de variance non-constante, mais on montre que celle-ci est fonction de \(\sigma_{\varepsilon}^2\) : \[\text{Var}(e_i)=\sigma_{\varepsilon}^2(1-h_i),\] où \(h_i\) ("hat value") est l'effet levier de la \(i\)ème observation sur l'ensemble des valeurs ajustées via le modèle de régression.

Les points exerçant un fort effet levier tendent donc à avoir de plus faibles résidus (puisqu'ils "tirent" la droite de régression vers eux).

Mesures d'influence et critères d'interprétation

• effet levier : large si \(h_i>2(k+1)/n\)
• résidus standardisés/studentisés : \(e'_i=e_i/s_e\sqrt{1-h_i}\) et \(e^*_i=e_i/s_{e^{(-i)}}\sqrt{1-h_i}\), où \(s_{e^{(-i)}}\) est une estimation de \(\sigma_{\varepsilon}\) sans l'observation \(i\). Comme \(e^*_i\sim t(n-k-2)\), on s'attend à ce que 5 % des valeurs au plus soient situées en dehors de l'intervalle [-2;+2].
• distance de Cook : \(D_i=\frac{h_i}{1-h_i}\times\frac{e_i^{'2}}{k+1}\), considérée large si \(D_i>4/(n-k-1)\).
• DFFITS\(=\frac{\hat y_i-\hat y_i^{(-i)}}{s_{e^{(-i)}}}\sqrt{h_i^{(-i)}}\), considéré large si \(>2\sqrt{p/n}\).
• DFBETAS\(=\frac{\hat\beta_j-\hat\beta_j^{(-i)}}{s_{e^{(-i)}}\sqrt{X'X}_{jj}}\), pour le \(j\)ème coefficient du modèle (incluant l'intercept), considéré large si \(>2/\sqrt{n}\).

La plupart de ces indices sont estimés en enlevant l'observation en question, et donc permettent d'apprécier son poids dans la qualité de la prédiction ou la stabilité des coefficients de régression. Voir http://bit.ly/JdaJM3.

Illustration

Les graphiques suivants visent essentiellement à mettre en évidence de potentielles valeurs extrêmes ou influentes ("outliers"). Ils reposent tous sur les données extraites de influence.measures. Voir aussi le package car (Fox, 2011) pour d'autres fonctionnalités graphiques.

Prédiction ponctuelle et par intervalle

Les valeurs prédites sont obtenues avec fitted ou predict (pour de nouvelles observations). Par exemple,

fitted(lm.fit)
predict(lm.fit, data.frame(HAEMO=seq(13, 18, by=1)))

Il existe deux types de prédictions par intervalle (95 %), de la forme \[\hat y\pm t_{n-p-1,1-\alpha/2}\text{SE}.\]

• Prédiction de \(y\) : predict(..., interval="prediction")
  \(\text{SE}(\hat y) = s\sqrt{1+\frac{1}{n}+\frac{(x-\bar x)^2}{\sum(x_i-\bar x)^2}}\)

• Prédiction de \(\mathbb{E}(y\mid x)\) : predict(..., interval="confidence")
  \(\text{SE}(\mathbb{E}(y\mid x))=s\sqrt{\frac{1}{n}+\frac{(x-\bar x)^2}{\sum(x_i-\bar x)^2}}\)

Illustration

Les intervalles de confiance de type prédiction sont plus larges que les intervalles associés à la prédiction de valeurs moyennes. Dans les deux cas, leur demi-largeur est toujours plus petite autour du point moyen \((\bar x, \bar y)\).

Voir aussi xYplot dans le package Hmisc.

Régression vs. ANOVA

La matrice de dessin (design) utilisée dans les deux cas est identique :

x <- gl(5, 1, 10, labels=letters[1:5])
model.matrix(rnorm(10) ~ x)

Le modèle s'écrit, pour chaque observation, \[y_i=\beta_0+\sum_{j=1}^k\beta_jx_i\] Soit les valeurs prédites, \[ \begin{array}{cc} \tilde y_{\cdot} &= b_0+b_1\times 0 + b_2\times 0 +b_3\times 0 + b_4\times 0\quad (x=a) \\ \tilde y_{\cdot} &= b_0+b_1\times 1 + b_2\times 0 +b_3\times 0 + b_4\times 0\quad (x=b) \\ \vdots & \vdots\\ \tilde y_{\cdot} &= b_0+b_1\times 0 + b_2\times 0 +b_3\times 0 + b_4\times 1\quad (x=e) \\ \end{array}\] c'est-à-dire les moyennes de groupes en ANOVA.

Application

Supposons que seules 10 classes soient disponibles pour la variable HAEMO.

haemo.dec <- cut2(paint$HAEMO, g=10)
fm <- PCV ~ haemo.dec
summary(aov.fit <- aov(fm, data=paint))
summary(lm.fit <- lm(fm, data=paint))
.Last.value$sigma^2
grp.means <- tapply(paint$PCV, 
                    haemo.dec, 
                    mean)   
grp.means[2:10] - grp.means[1]
coef(lm.fit)

La manière dont R code les contrastes (voir options("contrasts")) est importante !

contr.treatment(10)
contr.sum(10)
contr.helmert(10)

Usage des formules sous R

R suit les conventions de notation proposées par Wilkinson & Rogers (1973), discutées dans Chambers & Hastie (1992), pp. 24--30.

Notons x, y, ... des variables continues, a, b, ... des variables catégorielles. On utilise "~" pour dénoter une relation fonctionnelle entre une réponse, y, et

• x, régression linéaire simple
• x + 0 ou x -1, idem avec suppression de l'intercept
• a + b, deux effets principaux
• a * b, équivalent à 1 + a + b + a:b, idem avec interaction

Dans le cas des interactions, on a les équivalences suivantes :

• factor:factor \(\leftrightarrow \gamma_{ij}\) (croisement des niveaux des facteurs)
• factor:numeric \(\leftrightarrow \beta_jx\) (pente variable selon les niveaux du facteur)
• numeric:numeric \(\leftrightarrow \beta xz\) (produit élément par élément)

Usage des formules sous R (2)

• a + b + a:d, effets principaux pour A et B, et interaction AxD
• a * b * d - a:b:d, tous les effets sauf l'interaction AxBxD (inclut les interactions AxB, AxD et BxD)
• a / b, équivalent à 1 + a + b %in% a (relation d'emboitement)
• a / (b * d), tous les termes BxD pour chaque niveau de A
• .^2, tous les effets principaux et les interactions de second-ordre

Il est également possible de mettre à jour un modèle existant (en ajoutant ou supprimant des termes) :

fm <- y ~ x + a * b
mod1 <- lm(fm, data=dat)
update(mod1, . ~ . - a:b)  ## supprime interaction AxB

Analyse de variance à deux facteurs

Dans l'analyse de variance à deux facteurs considérés comme des effets fixes, on considère deux variables explicatives catégorielles et une variable réponse continue. À la différence de l'ANOVA à un facteur, l'inclusion d'un second facteur pose la question de la co-relation entre les deux facteurs dans la prédiction de la variable réponse. Cet effet d'interaction peut être ou non l'objet de l'étude, mais le plus souvent il est important de le quantifier et de le tester afin de pouvoir discuter les effets principaux des facteurs d'intérêt.

Traditionnellement, ce type d'analyse se retrouve dans les plans d'expérience (psychologie, agronomie, industrie, etc.) mais reste applicable aux études cliniques. Les conditions d'application de l'ANOVA sont les mêmes que dans le cas à un facteur, en particulier celles concernant les résidus.

La régression logistique ordinale (Armstrong & Sloan, 1989), des techniques de permutation (Anderson & Ter Braak, 2003), ou l'analyse median polish (Mosteller & Tukey, 1977) complémentent cette approche paramétrique.

Application

The effect of Vitamin C on tooth growth in Guinea Pigs. (Bliss, 1952)

Quelques statistiques descriptives (hors aggregate) :

data(ToothGrowth)
ToothGrowth$dose <- factor(ToothGrowth$dose)
fm <- len ~ supp * dose
replications(fm, data=ToothGrowth)
library(Hmisc)
f <- function(x) apply(x, 2, function(x) 
                       c(mean=mean(x), sd=sd(x)))
summary(fm, data=ToothGrowth, fun=f)

Notons que nous avons délibéremment ignoré le fait que le facteur dose est ordonné. Autre méthode de calcul des moyennes conditionnelles et marginales (ajouter margins=TRUE) :

library(reshape2)
m <- acast(ToothGrowth, supp ~ dose, mean, value.var="len")

Graphique d'interaction

L'inspection visuelle de l'éventuel effet d'interaction peut se faire avec un simple diagramme en lignes. Le non-parallélisme des deux droites VC et OJ suggère que l'effet dose n'est pas le même selon le type de supplément.

xyplot(len ~ dose, data=ToothGrowth, groups=supp,
       type=c("p","a"))

Analyse du modèle complet

aov.fit <- aov(fm, data=ToothGrowth)
summary(aov.fit)
model.tables(aov.fit, type="means", se=TRUE, 
             cterms="supp:dose")

Le test F correspondant à l'interaction supp:dose est significatif à 5 %, suggérant que l'on ne peut se limiter à l'interprétation seule des effets principaux : l'effet dose n'est pas le même selon le type de supplément (VC ou OJ). Pour quantifier les différences moyennes entre traitements :

apply(m, 2, diff)

Voir aussi le package effects.

Vérification des conditions d'application

Concernant la normalité, on peut inspecter la distribution des résidus avec un histogramme ou des quantiles. L'homogénéité des variances peut être vérifiée avec des boîtes à moustaches (et un test statistique).

qqmath(~ resid(aov.fit))
bwplot(len ~ interaction(supp, dose), data=ToothGrowth)
bartlett.test(len ~ interaction(supp,dose),data=ToothGrowth)

Analyse de covariance

L'analyse de covariance consiste à tester différents niveaux d'un facteur en présence d'un ou plusieurs co-facteurs continus. La variable réponse et ces co-facteurs sont supposées reliés, et l'objectif est d'obtenir une estimation des réponses corrigée pour les éventuelles différences entre groupes (au niveau des cofacteurs).

Ce type d'analyse est fréquemment utilisé dans le cas des données pré/post avec des mesures continues et un facteur de groupe, et reste préférable à une simple analyse des scores de différences (Miller & Chapman, 2001; Senn, 2006).

Il existe également des alternatives non-paramétriques (Young & Bowman, 1995), disponibles dans le package sm (sm.ancova).

Application

Weight change data for young female anorexia patients. (Hand, Daly, McConway, & Ostrowski, 1993)

data(anorexia)
anorexia$Treat <- relevel(anorexia$Treat, ref="Cont")
anorex.aov0 <- aov(Postwt ~ Prewt + Treat, data=anorexia)
anorex.aov1 <- aov(Postwt ~ Prewt * Treat, data=anorexia)
summary(anorex.aov0)

En tenant compte du poids initial, le modèle de base (anorex.aov0) suggère bien qu'il existe des différences de poids moyens entre les trois groupes. Ce modèle suppose que la relation entre poids avant/après traitement est indépendante du traitement.

Vérification graphique

Un diagramme de dispersion faisant apparaître les différents groupes permet de vérifier visuellement si l'hypothèse de parallélisme est réaliste. Ici, clairement cette hypothèse ne tient pas.

xyplot(Postwt ~ Prewt, data=anorexia, groups=Treat, 
       aspect="iso", type=c("p","r"))

Test de parallélisme

Comparaison avec un modèle incluant l'interaction Prewt:Treat :

anova(anorex.aov0, anorex.aov1)

Le test pour le terme d'interaction indique que celle-ci est significative ; par conséquent, on ne peut pas considérer que la relation entre poids avant et après traitement est la même dans les trois groupes. Il faudra donc décrire les résultats groupe par groupe.

summary.lm(anorex.aov1)

On peut vouloir plutôt tester \(H_0:\beta=1\), d'où

lm(Postwt ~ Prewt + Treat + offset(Prewt), data=anorexia)

Interprétation des coefficients du modèle

Le modèle sans interaction (coef(anorex.aov0)) s'écrit \[\tilde y_i = 45.67 + 0.43\times\text{Prewt}_i +4.10\times\mathbb{I}(\text{Treat}_i=\text{CBT}) +8.66\times\mathbb{I}(\text{Treat}_i=\text{FT}).\] Pour les patientes du groupe contrôle, \(\tilde y_i = 45.67 + 0.43\times\text{Prewt}_i\), alors que pour celles du groupe FT, \(\tilde y_i = 45.67 + 0.43\times\text{Prewt}_i{\color{myred2}+8.66}\). Ceci correspond bien à l'idée que l'effet de Prewt est le même pour toutes les patientes et que le facteur de groupe induit simplement un changement moyen (+4.10 ou +8.66) par rapport au groupe contrôle.

Pour le modèle avec interaction avec Prewt centré, on a
\[\begin{array}{rl}\tilde y_i =& 80.99 - 0.13\times\text{Prewt}_i\\ &\phantom{80.99}+4.46\times\mathbb{I}(\text{Treat}_i=\text{CBT})\\ &\phantom{80.99}+8.75\times\mathbb{I}(\text{Treat}_i=\text{FT})\\&\phantom{80.99}+0.98\times\text{Prewt}_i\times\mathbb{I}(\text{Treat}_i=\text{CBT})\\& \phantom{80.99}+1.04\times\text{Prewt}_i\times\mathbb{I}(\text{Treat}_i=\text{FT}).\end{array}\]

Régression linéaire multiple

L'inclusion de plusieurs prédicteurs conduit à généraliser le modèle de régression simple à la régression multiple et à s'intéresser à l'effet de chaque prédicteur en tenant compte des autres variables dans le modèle : on parlera de l'effet d'un prédicteur à niveau constant des autres prédicteurs. Les coefficients de régression reflète toujours le poids de chaque variable explicative, via le coefficient de corrélation partiel.

La régression multiple pose la question de la sélection des prédicteurs dans le modèle final, et de l'évaluation des effets partiels de chaque prédicteur. L'analyse des résidus du modèle doit tenir compte de l'ensemble des prédicteurs en présence dans le modèle.

Fox (2011) décrit en détails les différentes étapes de construction et de validation d'un modèle de régression multiple. Voir aussi http://tinyurl.com/carbook.

Cas des données corrélées

Tous les modèles discutés précédemment supposaient les observations indépendantes les unes des autres. Dans le cas où la même unité statistique est utilisée pour collecter différentes mesures (mesures répétées), il est important de prendre en compte la corrélation intra-unité induite par ce type de recueil de données. Cela présente l'avantage d'offir un gain de puissance statistique et la possibilité de modéliser la structure de covariance.

L'ANOVA dite "à mesures répétées" permet, en incorporant un effet sujet aléatoire, d'incorporer la corrélation intra-unité et de mieux estimer la résiduelle. On fait une hypothèse de symétrie composée selon laquelle la covariance intra-unité est la même pour toutes les unités statistiques. Des résultats équivalents sont obtenus par un modèle linéaire mixte à intercept aléatoire.

Application

Blood cholesterol concentrations. (Zar, 1998)

chs <- read.table("cholesterol.txt", header=TRUE)
chs$Subject <- factor(chs$Subject)  ## important
chs <- melt(chs, id.vars="Subject")
aov1 <- aov(value ~ variable + Error(Subject), data=chs)
summary(aov1)

Le terme d'erreur pour tester l'effet intra-sujet est calculé en enlevant la variabilité inter-sujets (SS=18731, soit 96.4 % de la résiduelle d'un modèle à un facteur). Le résultat significatif du test suggère que le niveau de cholestérol varie bien selon le traitement.

Application (2)

Un modèle à effet aléatoire donnera sensiblement le même résultat, tout en permettant d'estimer la corrélation intraclasse :

library(nlme)
lme1 <- lme(value ~ variable, data=chs, 
            random= ~ 1 | Subject)
summary(lme1)
anova(lme1)
31.96^2/(31.96^2+7.61^2)  ## ICC
intervals(lme1)

Illustration

On vérifie aisément dans le graphique des valeurs prédites que la corrélation intra-sujet est identique pour tous les sujets et que seuls les niveaux moyens varient entre les sujets.

Régression multiple et extensions

Objectifs

• Dans ce cours, on approfondira les notions abordées pour le modèle de régression multiple, en particulier l'interprétation des coefficients et la comparaison de modèles.
• En présence d'un grand nombre de variables explicatives, la question du choix des variables (sélection de variables) se pose naturellement et différentes techniques peuvent être utilisées.
• Le traitement des données manquantes sera également traité séparément.

Quelques éléments de contexte

On distingue essentiellement deux types d'approche :

• explicative : construire un modèle permettant de rendre compte des relations causale entre des prédicteurs et une variable réponse
• prédictive : construire un modèle permettant de prédire avec un minimum d'erreur (régression ou classification) les réponses observées ou futures ; pas nécessairement de relations de causalité

Concernant les modèles multivariés prédictifs, il peut s'agir d'études diagnostiques ou pronostiques. Dans tous les cas, il existe des recommendations concernant la sélection de variables ou l'évaluation des performances du modèle (Bouwmeester et al., 2012; Harrell, 2001; Steyerberg, 2009).

Recommendations et checklists : CONSORT, STARD, STROBE, REMARK.

Rappel sur le modèle de régression multiple

Dans le cas où la variable réponse est continue, le modèle de régression multiple s'écrit simplement comme une combinaison linéaire de différents termes représentant les effets de chaque régresseur : \[\mathbb{E}(y\mid x_1,x_2,\dots) = \beta_0 + \sum_{j=1}^p\beta_jx_p.\]

Les coefficients de régression représentent l'effet partiel de chaque prédicteur, c'est-à-dire l'effet de l'augmentation d'une unité de \(x_j\) sur \(y\) lorsque tous les autres prédicteurs sont maintenus constants. On peut tester l'égalité de l'ensemble des coefficients (test F), la nullité d'un coefficient (test t) ou d'un ensemble de coefficients (test F partiel).

Interprétation des coefficients

Par exemple, si \(\hat y=\hat\beta_0+\hat\beta_1x_1+\hat\beta_2x_2\), avec \(x_1\) l'âge et \(x_2\) le sexe (0=M, 1=F), \[\begin{array}{ll}\mathbb{E}(y|x_1,x_2)=\hat\beta_0+\hat\beta_1x_1,& \text{pour les hommes,}\\ \mathbb{E}(y|x_1,x_2)=(\hat\beta_0+\hat\beta_2)+\hat\beta_1x_1,& \text{pour les femmes.}\end{array}\]

d'où \(\beta_1\) qui représente l'effet lorsqu'on augmente \(X_1\) de une unité, à niveau constant des autres facteurs (ici, \(X_2\)).

Rappel sur la corrélation partielle : \(r_{xy.z}=\frac{r_{xy}-r_{xz}r_{yz}}{\sqrt{(1-r_{xz}^2) (1-r_{yz}^2)}},\) où \(r_{xz}^2\) s'interprète comme la part de variance expliquée par \(X\) lorsque \(Z\) est maintenu constant.

Analyse des résidus

Vérifier :

• homoskédascticité : \(\sigma^2\) constante, indépendament des \(x_j\)
• linéarité des effets : l'effet de chaque prédicteur est linéaire, lorsque les autres prédicteurs sont maintenus constants
• présence d'outliers ou de points influents
• normalité et indépendance des erreurs (pour l'inférence)

En plus des techniques d'analyse des résidus vues dans le cas de la régression linéaire simple s'ajoutent les graphiques des résidus partiels : residus lorsqu'on prédit \(y\) à partir de tous les \(x\) sauf \(x_j\) vs. résidus lorsqu'on prédit \(x_j\) à partir tous les \(x\) sauf \(x_j\), ce qui revient à se demander comment ce que l'on ne peut prédire de \(y\) en ignorant \(x_j\) pourrait dépendre de ce que l'on ne peut pas prédire de \(x_j\) à partir des autres \(x\).

Transformations

text

Co-facteurs généralement d'intérêt

• "patient outcome studies" et facteurs de risque (Iezzoni, 1994) :
  age, sexe, stabilité clinique, diagnostique, séverité du diagnostique, comorbidités (nombre et sévérité), état physique fonctionnel, status psychologique, cognitive et psycho-social, données culturelles, ethniques et socio-économiques, santé générale et qualité de vie, préférences patient

Qualité d'ajustement

text

Choix des prédicteurs dans un modèle

On peut se demander

• préalablement ou durant la construction d'un modèle multivarié, quelles sont les variables à inclure dans le modèle, par exemple lorsqu'on cherche à isoler des prédicteurs d'intérêt dans le cas où \(n\ll p\) : il s'agit d'un problème de sélection de variables ;
• après avoir construit un premier modèle, quel sous-ensemble des prédicteurs d'origine fournit un sous-modèle essentiellement "équivalent" en termes de qualité de prédiction : il s'agit d'un problème de sélection de modèles.

Procédures de sélection automatique de prédicteurs

• Tous les sous-ensembles possibles (Leaps and Bounds) :
  exploration de tous les sous-modèles possibles (\(2^p-1\)), et sélection du meilleur modèle sur la base du \(C_p\) de Mallows (\(C_p=(\text{RSS}_kp/s^2)-(n-2p)\), pour un modèle à \(p-1\) prédicteurs).
• Méthodes pas-à-pas
  • ascendante : à partir d'un modèle de base contenant l'intercept, ajouter des variables tant que celles-ci sont compatibles avec le critère d'inclusion (généralement basé sur une diminution de la résiduelle, ou une p-valeur type p=0.1)
  • descendante : même chose, mais en partant du modèle complet (tous les prédicteurs)
  • mixte : principe de l'approche ascendante (descendante), avec possibilité d'enlever (ajouter) les variables précédemment incluses (omises), et utilisation de critères d'inclusion/exclusion potentiellement différents (p=0.1 pour l'entrée, p=0.2 pour la sortie).

Problèmes des méthodes pas-à-pas

Les méthodes pas-à-pas posent plusieurs problèmes : sur-estimation des \(R^2\), les statistiques usuelles (\(F\) ou \(\chi^2\)) ne suivent plus les distributions théoriques, sous-estimation des erreurs standard des coefficients de régression (eux-mêmes biasés vers les valeurs élevées), \(p\)-valeurs trop libérales (comparaisons multiples), pas de solution à la multi-colinéarité (Harrell, 2001, pp. 56–60).

Par ailleurs, on ne peut utiliser que les cas complets. Les procédures ascendante et mixte ne permettent pas de traiter le cas des effets de confusion "négatifs" (lorsque ceux-ci amènent à sous estimer l'effet d'un facteur d'exposition sur la maladie, par exemple).

Enfin, on a toujours le problème de la taille de l'échantillon : avec de petits échantillons, on manque de puissance et de précision et on risque de ne pas inclure des prédicteurs importants ; avec de grands échantillons, des prédicteurs de peu d'intérêt peuvent se retrouver inclus.

Application

Framingham Heart Study. (Dupont, 2009, §3.17; Levy, 1999)

fhs <- read.csv("Framingham.csv")
fhs <- subset(fhs, fhs$id != 9999 & complete.cases(fhs))
fhs$sex <- fhs$sex-1
fm <- log(sbp) ~ log(bmi) + age + log(scl) + sex + 
               sex:log(bmi) + sex:age + sex:log(scl)
mod0 <- lm(log(sbp) ~ 1, data=fhs)
mod1 <- lm(fm, data=fhs)
summary(mod1)

Application (2)

On obtient des résultats différents par ajout ou suppression de termes reposant sur le critère AIC (Venables & Ripley, 2002, pp. 172–177).

addterm(mod0, mod1)    ## ou add1()
drop1(mod1, test="F")

En revanche, les méthodes pas à pas donnent les mêmes résultats :

st.fwbw <- step(mod1)  ## backward + forward
preds <- attr(terms(fm), "term.labels")
preds <- paste("~", paste(preds, collapse="+"))
step(mod0, scope=as.formula(preds), direction="forward")

Le cas des données manquantes

Les données manquantes réduisent la puissance des tests (réduction d'effectif dans les analyses cas complets) et peuvent induire des biais au niveau de l'interprétation dans certains cas (p. ex., s'il y a beaucoup de valeurs dans un seul groupe detraitement)

On distingue trois types de VM :

Régression logistique

Objectifs

• Dans ce cours, on aborde la régression logistique (simple et multiple) qui sert à modéliser la relation entre des prédicteurs catégoriels ou continus et une variable réponse binaire.
• Sans rentrer dans les détails des modèles linéaires généralisés, on insistera sur la correspondance entre la régression logistique et la régression linéaire sur le plan de la modélisation, tout en soulignant les différences les plus notables : interprétation des coefficients, qualité d'ajustement du modèle, analyse des résidus.

Lectures conseillées : Vittinghoff et al. (2005), Hosmer & Lemeshow (1989).

Régression logistique

La régression logistique permet de traiter le cas où la variable réponse est de type binaire (oui/non, malade/pas malade, etc.), et non pas continu comme dans le modèle de régression linéaire. Tout en relaxant certaines des hypothèses du modèle de régression multiple, on maintient quand même l'idée d'une relation linéaire entre la réponse et les prédicteurs.

Dans le cas d'une variable binaire, sa "moyenne" correspond à la proportion d'individus possédant la caractéristique étudiée ou répondant positivement à l'événement, d'où l'idée de modéliser la probabilité de succès, comprise entre 0 et 1, en fonction d'un certain nombre de prédicteurs.

Dans les enquêtes épidémiologiques cas-témoin (avec ou sans matching) où l'incidence de la maladie n'est pas connue, la régression logistique fournit une estimation de l'odds-ratio ajusté sur les co-facteurs d'intérêt (âge, sexe, etc.). D'autre part, lorsque la prévalence de la maladie est faible, l'OR fournit une bonne approximation du risque relatif.

Illustration

L'exemple ci-dessous montre les réponses observées (0/1) et les probabilités individuelles correspondantes, telles qu'estimées à partir d'une régression logistique. En considérant \(\tilde y_i=\mathbb{I}(P(x)\geq 0.5)\), on dénombre 8 mal classés (10 %).

Parallèle avec la régression linéaire

Comme dans la régression linéaire, on cherche la meilleure combinaison linéaire des données d'entrée pour modéliser la réponse, à ceci près que c'est une transformation de cette combinaison (on parle d'une fonction de lien) qui est utilisée en sortie.

Différences avec le modèle linéaire

L'analyse des résidus du modèle permet de vérifier si celui-ci est satisfaisant (en terme de spécification et de qualité d'ajustement).

Les différences principales avec la régression linéaire sont les suivantes :

• On ne parle plus de sommes de carrés (OLS, résidus, variance) mais de déviance (dans le cas gaussien, elle est équivalente à la somme de carrés de la résiduelle), mais cette dernière reflète toujours l'écart entre les données et le modèle.
• En raison de la nature binaire de la variable réponse, l'analyse classique des résidus en fonction des valeurs prédites ou la notion d'hétéroskédasticité ne font plus sens ; en revanche, on s'intéresse toujours à la qualité d'ajustement du modèle, et à la comparaison de modèles emboîtés qui permettent d'évaluer l'apport d'un ou plusieurs prédicteurs par rapport à un modèle de base. Pour cela, on utilise des tests de rapport de vraisemblance.

Interprétation des coefficients

Le terme d'intercept s'interprète comme un odds, et les coefficients de régression comme des odds-ratio : lorsque \(X_j\) augmente de \(d=1\) unité, l'odds de \(y=1\) augmente de \(\exp(\beta_jd)\) (de manière équivalente, le log-odds augmente de \(\beta_jd\)). Dans le cas où l'on a un seul prédicteur, binaire, on peut vérifier à partir de la relation \(\frac{P(x)}{1-P(x)}=\exp(\beta_0+\beta_1x)\) que
\(\frac{P(1)/[1-P(1)]}{P(0)/[1-P(0)]}=\exp(\beta_1).\)

\[\begin{array}{l} \log\left[\frac{P(56)}{1-P(56)}\right]-\log\left[\frac{P(55)}{1-P(55)}\right]\\ \phantom{==}=(-5.940+0.074\times 56)-(-5.940+0.074\times 55)=0.074. \end{array}\] D'où un odds-ratio de \(\exp(0.074)=1.077\) associé à une augmentation d'âge d'un an sur le risque d'infarctus, c'est-à-dire une augmentation du risque de 8 %. Pour une variation de 10 ans, l'OR est de \(\exp(0.074\times 10)=2.105\).

Données groupées vs. non groupées

Dans le cas où on dispose des données individuelles, soit n réponses binaires codées sous forme numérique 0/1 ou à l'aide d'un facteur dont le premier niveau code l'échec (0), le modèle de régression logistique s'écrit :

glm(y ~ x, family=binomial)

Dans le cas des données "groupées" (aggrégées par niveaux d'un ou plusieurs facteurs), on utilise une matrice à deux colonnes représentant les effectifs 0/1 (n0=n-n1) :

glm(cbind(n1, n0) ~ x, family=binomial)

ou directement la fréquence empirique et le nombre d'observations :

glm(n1/n ~ x, family=binomial, weights=n)

On peut toujours passer d'un format à l'autre avec aggregate ou untable (package reshape).

Qualité d'ajustement du modèle et diagnostics

La déviance et le \(\chi^2\) généralisé de Pearson sont les deux outils principaux pour l'évaluation et la comparaison de modèles logistiques. La déviance se définit comme deux fois la différence entre les log vraisemblances de deux modèles, tandis que \(\chi^2=\sum_i(y_i-\hat y_i)^2/V(\hat y)\). (On peut aussi définir des résidus basés sur la déviance.)

Le test de Hosmer-Lemeshow consiste à évaluer la concordance entre les valeurs prédites et observées des observations regroupées en quantiles, typiquement des déciles. Ce test dépend du nombre de groupes fixés a priori, et il est peu puissant en cas de mauvaise spécification. Les techniques de lissage non-paramétrique ou d'estimations stratifiées sont utiles pour identifier des déviations locales ou globales (Harrell, 2001).

Il existe des mesures de type pseudo-\(R^2\) (Cox-Snell, Nagelkerke) qui permettent d'évaluer la qualité du modèle global (par rapport au modèle nul). Concernant la capacité prédictive du modèle, on utilise généralement le score de Brier ou l'index \(C\) (probabilité de concordance), en lien avec la courbe ROC.

Qualité d'ajustement du modèle et diagnostics (2)

Comme dans le cas du modèle linéaire, on s'intéressera aux points influents et aux "outliers", tous deux dérivés à partir des résidus standardisés de Pearson. Le problème est que dans ce cas, ces résidus appartiennent à deux classes ce qui rend délicat, voire inutile, les représentations graphiques de type résidus vs. valeurs prédites. On se concentrera donc plutôt sur les mesures d'influence, comme les DFBETAS. On peut représenter ces dernières en fonction des probabilités prédites. Voir aussi residual.plots (package alr3). Il est également possible d'examiner les variations dans le \(\chi^2\) (ou la déviance) en fonction du leverage ou des probabilités prédites, voir par exemple http://bit.ly/JMtaho.

Concernant la linéarité de la relation entre les prédicteurs et le log-odds, on peut utiliser de simple diagramme de dispersion, ou utiliser ce que l'on appelle les "marginal model plot" (par ailleurs plus utile dans le cas des données non groupées), et on peut juger de la linéarité de la relation à partir d'une courbe lowess. Voir mmps (package car) ainsi que le package marginalmodelplots.

Sélection de variables, modèles prédictifs

La construction et la validation de modèles prédictifs dans le domaine clinique ont fait l'objet de nombreux ouvrages, en particulier Harrell (2001) et Steyerberg (2009) (voir aussi Steyerberg et al. (2001)). Le site RMS fournit de nombreuses ressources sur ce sujet : http://biostat.mc.vanderbilt.edu/wiki/Main/RmS

Pour un tour d'horizon rapide des enjeux de la modélisation à partir d'enquêtes épidémiologiques, voir Greenberg & Kleinbaum (1985). De nombreux autres articles fournissent des recommendations pour le reporting des résultats (Bagley, White, & Golomb, 2001; Bouwmeester et al., 2012; K. J. Ottenbacher, Ottenbacher, Tooth, & Ostir, 2004).

Application 1

Heart disease and blood pressure. (Everitt & Rabe-Hesketh, 2001, p. 208)

bp <- read.table("hdis.dat", header=TRUE)
blab <- c("<117","117-126","127-136","137-146",
          "147-156","157-166","167-186",">186")
clab <- c("<200","200-209","210-219","220-244",
          "245-259","260-284",">284")
bp <- within(bp, {
  bpress <- factor(bpress, labels=blab)
  chol <- factor(chol, labels=clab)
})
round(xtabs(hdis/total ~ bpress + chol, data=bp), 2)

Ici, les données sont groupées et on cherche à modéliser la probabilité de survenue d'une crise cardiaque en fonction de la tension artérielle, que l'on peut considérer comme une variable ordinale, ou numérique en considérant les valeurs centrales des catégories.

Résultats

midpoint <- function(x) {
  x <- as.numeric(unlist(strsplit(x, "-")))
  return(sum(x)/2)
}
val <- sapply(levels(bp$bpress)[-c(1,8)], midpoint)
dfrm <- aggregate(bp[,3:4], list(bpress=bp[,1]), sum)
dfrm$bpress <- c(val[1]-10, val, val[6]+15)
mod1 <- glm(cbind(hdis, total-hdis) ~ bpress, 
            data=dfrm, family=binomial)
summary(mod1)

Illustration

Ci-dessous sont représentées les valeurs observées et prédites, sous forme de proportion ou de logit. Dans le deuxième cas (droite), on voit clairement la linéarité entre le log odds, \(\hat g(x)\), et \(x\).

16/271 = 0.059
log(0.059/(1-0.059)) = -2.769
1/(1+exp(-(-6.082+0.024*141.5))) = 0.067

Application 2

The low birth weight study. (Hosmer & Lemeshow, 1989)

data(birthwt, package="MASS")
ethn <- c("White","Black","Other")
birthwt$race <- factor(birthwt$race, labels=ethn)
fm <- low ~ age + lwt + race + ftv
glm1 <- glm(fm, data=birthwt, family=binomial)
summary(glm1)
confint(glm1)

L'âge de la mère ne semble pas être un prédicteur intéressant, mais son poids oui, et on retrouve des différences entre les mères White vs. Black.

Interprétation des coefficients

Le logit estimé s'écrit (Hosmer & Lemeshow, 1989 p. 36) \[\begin{array}{rl}\hat g(x)=&1.295-0.024\times\text{age}-0.014\times\text{lwt}\\ &\phantom{1.295}+1.004\times\mathbb{I}(\text{race}=\text{Black})+0.433\times\mathbb{I}(\text{race}=\text{Other})\\&\phantom{1.295}-0.049\times\text{ftv},\end{array}\] ce qui permet d'obtenir les valeurs prédites ("probabilité d'un enfant avec un poids inférieur à 2,3 kg à la naissance") facilement. Pour un enfant dont la mère possède les caractéristiques moyennes ou modales de l'échantillon (âge 23,2 ans, poids 59 kg, White, aucune visite chez le gynécologue durant le premier trimestre de grossesse), on obtient : \[\begin{array}{rl}P(\text{low}=1)&=1/\big(1+\exp(-(1.30-0.02\times 23.24-0.01\times 129.82))\big)\\&=0.253.\end{array}\]

log.odds <- predict(glm1, data.frame(age=mean(birthwt$age), 
                                     lwt=mean(birthwt$lwt),
                                     race="White", ftv=0))
exp(log.odds)/(1+exp(log.odds))  ## 0.2483566

Significativité des termes du modèle

Les facteurs age et ftv apparaissent non significatifs. On peut a priori simplifier le modèle comme suit :

glm2 <- update(glm1, . ~ . - age - ftv)

Ces deux modèles sont bien emboîtés (glm2 opère juste une restriction dans l'espace des prédicteurs).

À partir de l'analyse des déviances résiduelles, on pourrait confirmer que ce modèle n'est pas "moins bon" que le précédent (\(\chi^2(2)=0.686\), \(p=0.710\)) :

anova(glm2, glm1, test="Chisq")

Significativité des termes du modèle (2)

En raison des multiples degrés de liberté, l'évaluation des tests de significativité (Wald) associés aux coefficients de régression des variables catégorielles est délicat.

Pour évaluer la significativité du facteur race, on peut procéder par comparaison de modèles également et le test du rapport de vraisemblance s'obtient comme suit :

anova(update(glm2, . ~ . - race), glm2)

D'où un test non significatif à 5 % :
(2 degrés de liberté)

pchisq(5.4316, 2, lower.tail=FALSE)

Toutefois, ce facteur étant connu pour être cliniquement important, on peut préférer le conserver dans le modèle final même si à l'évidence il ne change pas significativement la qualité explicative ou prédictive du modèle.

Vérification des hypothèses du modèle

Pour la linéarité, on peut examiner les prédictions (probabilités) en fonction de \(g(x)\), ou des variables continues.

library(car)
mmps(glm2, terms=~lwt)

Vérification des hypothèses du modèle (2)

Les résidus (Pearson ou déviance) peuvent être obtenus à partir de la fonction rstandard. Quant aux mesures d'influence, il s'agit du même principe que pour le modèle linéaire :

influence.measures(glm2)

Fonctionnalités avancées

La fonction lrm du package rms (Harrell, 2001) est beaucoup plus informative que glm, en particulier en ce qui concerne la qualité d'ajustement du modèle et son pouvoir discriminant.

library(rms)
ddist <- datadist(birthwt)
options(datadist="ddist")
glm2b <- lrm(low ~ lwt + race, data=birthwt)
glm2b
exp(coef(glm2b))

Prédictions

Le package rms fournit également tout un ensemble de fontions additionnelles pour les prédictions et leur représentation graphique (à l'image du package effects).

Predict(glm2b, lwt, race, fun=plogis, conf.type="mean")

Régression logistique et diagnostique médical

Le modèle de régression logistique est idéal pour fournir des probabilités (au niveau individuel, ou en moyenne conditionnellement à certains co-facteurs). Dans un contexte de classification, on réduit ces prédictions à deux classes, 0/1, à partir d'un "cut-off" optimisé sur la base d'un compromis entre sensibilité et spécificité, sans prendre en considération un coût différentiel pour les faux-positifs et les faux-négatifs. D'autre part, transformer une probabilité (continue dans [0;1]) en une variable binaire, sans autoriser de zone d'incertitude, est risqué.

Pour plus d'informations sur la présentation des résultats d'un modèle de régression logistique dans les études diagnostiques, voir Steyerberg (2009).

There is a reason that the speedometer in your car doesn't just read "slow" and "fast". Frank Harrell on R-help, February 2011

Voir aussi Direct Measures of Diagnostic Utility Based on Diagnostic Risk Models, http://1.usa.gov/Mba7xK.

Comparaison avec d'autres méthodes (classification)

Ci-dessous figurent quelques résultats concernant l'utilisation de différents modèles de classification et la comparaison de leur capacité prédictive. Les taux de classification correcte sont estimés par validation croisée (25 x 10 fold), avec le package caret (Kuhn, 2008).

Sur les données d'origine, le taux de classification correcte est de 68.8 % :

table(predict(glm1, type="resp")>=.5, birthwt$low)

Attention, il s'agit dans tous les cas d'un taux de classification "apparent" !

Agresti, A. (2002). Categorical Data Analysis. John Wiley & Sons.

Altman, D. G., & Bland, J. M. (1983). Measurement in medicine: the analysis of method comparison studies. Statistician, 32, 307–317.

Anderson, M. J., & Ter Braak, C. J. F. (2003). Permutation tests for multi-factorial analysis of variance. Journal of Statistical Computation and Simulation, 73, 85–113.

Armstrong, B. G., & Sloan, M. (1989). Ordinal regression models for epidemiological data. American Journal of Epidemiology, 129, 191–204.

Bagley, S. C., White, H., & Golomb, B. A. (2001). Logistic regression in the medical literature: Standards for use and reporting, with particular attention to one medical domain. Journal of Clinical Epidemiology, 54, 979–985.

Becker, R. A., & Chambers, J. M. (1981). S: A Language and System for Data Analysis.

Becker, R. A., & Chambers, J. M. (1984). S: An Interactive Environment for Data Analysis and Graphics. Wadsworth.

Becker, R. A., Cleveland, W. S., & Shyu, M. J. (1996). The Visual Design and Control of Trellis Display. Journal of Computational and Statistical Graphics, 5, 123–155.

Bishop, Y. M., Fienberg, S. E., & Holland, P. W. (2007). Discrete Multivariate Analysis. Springer.

Bliss, C. I. (1952). The Statistics of Bioassay. Academic Press.

Bouwmeester, W., Zuithoff, N. P. A., Mallett, S., Geerlings, M. I., Vergouwe, Y., Steyerberg, E. W., Altman, D. G., et al. (2012). Reporting and Methods in Clinical Prediction Research: A Systematic Review. PLoS Medicine, 9, 1001221.

Burns, P. (2011, April). The R Inferno. Retrieved from http://www.burns-stat.com

Campbell, I. (2007). Chi-squared and Fisher-Irwin tests of two-by-two tables with small sample recommendations. Statistics in Medicine, 26, 3661–3675.

Chambers, J. M., & Hastie, T. J. (Eds.). (1992). Statistical Models in S. Wadsworth & Brooks.

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74, 829–836.

Cleveland, W. S. (1985). The Elements of Graphing Data. Monterey, CA: Wadsworth.

Cleveland, W. S. (1993). Visualizing Data. Hobart Press.

Dalgaard, P. D. (2008). Introductory Statistics with R. Springer.

Dupont, W. D. (2009). Statistical Modeling for Biomedical Researchers. Cambridge University Press.

Everitt, B., & Rabe-Hesketh, S. (2001). Analyzing Medical Data Using S-PLUS. Springer.

Falissard, B. (2005). Comprendre et utiliser les statistiques dans les sciences de la vie. Masson.

Fox, J. (2011). An R Companion to Applied Regression. Sage Publications.

Friendly, M. (2000). Visualizing Categorical Data. SAS Publishing.

Friendly, M. (2011). Tutorial: Working with categorical data with R and the vcd package. Retrieved from http://www.datavis.ca/courses/VCD

Greenberg, R. S., & Kleinbaum, D. G. (1985). Mathematical modeling strategies for the analysis of epidemiological research. Annual Review of Public Health, 6, 223–245.

Hand, D. J., Daly, F., McConway, K., & Ostrowski, E. (Eds.). (1993). A Handbook of Small Data Sets. Chapman & Hall.

Harrell, F. E. (2001). Regression modeling strategies with applications to linear models, logistic regression and survival analysis. Springer.

Hollander, M., & Wolfe, D. A. (1999). Nonparametric Statistical Methods. Wiley.

Hosmer, D., & Lemeshow, S. (1989). Applied Logistic Regression. New York: Wiley.

Iezzoni, L. I. (1994). Using risk-adjusted outcomes to assess clinical practice: An overview of issues pertaining to risk adjustment. Annals of Thoracic Surgery, 58, 1822–1826.

Ihaka, R., & Gentleman, R. (1996). R: A language for data analysis and graphics. Journal of Computational and Graphical Statistics, 5, 299–314.

Kuhn, M. (2008). Building Predictive Models in R Using the caret Package. Journal of Statistical Software, 28.

Levy, D. (1999). 50 Years of Discovery: Medical Milestones from the National Heart, Lung, and Blood Institute’s Framingham Heart Study.

Logan, M. (2010). Biostatistical Design and Analysis using R. Wiley.

Meyer, D., Zeilis, A., & Hornik, K. (2006). The Strucplot Framework: Visualizing Multi-way Contingency Tables with vcd. Journal of Statistical Software, 17.

Miller, G. A., & Chapman, J. P. (2001). Misunderstanding Analysis of Covariance. Journal of Abnormal Psychology, 110, 40–48.

Mosteller, F., & Tukey, J. (1977). Data Analysis and Regression. Reading, MA: Addison-Wesley.

Murrell, P. (2005). R Graphics. Chapman & Hall/CRC.

Ottenbacher, K. J., Ottenbacher, H. R., Tooth, L., & Ostir, G. V. (2004). A review of two journals found that articles using multivariable logistic regression frequently did not report commonly recommended assumptions. Journal of Clinical Epidemiology, 57, 1147–1152.

Paradis, E. (2005, Septembre). R pour les débutants. Retrieved from http://cran.r-project.org

Sarkar, D. (2008). Lattice, Multivariate Data Visualization with R. Springer.

Senn, S. (2006). Change from baseline and analysis of covariance revisited. Statistics in Medicine, 25, 4334–4344.

Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlation: Uses in assessing rater reliability. Psychological Bulletin, 86, 420–428.

Silverman, B. W. (1986). Density estimation. Chapman and Hall.

Steyerberg, E. W. (2009). Clinical prediction models: a practical approach to development, validation, and updating. Springer.

Steyerberg, E. W., Harrell, F. E., Borsboom, G. J. J. M., Eijkemans, M. J. C., Vergouwe, Y., & Habbema, J. D. F. (2001). Internal validation of predictive models: Efficiency of some procedures for logistic regression analysis. Journal of Clinical Epidemiology, 54, 774–781.

Student. (1908). The probable error of a mean. Biometrika, 6, 1–25.

Sturges, H. A. (1926). The choice of a class interval. Journal of the American Statistical Association, 65–66.

Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.

Venables, W. N., & Ripley, B. D. (2002). Modern Applied Statistics with S. Springer.

Vittinghoff, E., Glidden, D. V., Shiboski, S. C., & McCulloch. (2005). Regression Methods in Biostatistics. Linear, Logistic, Survival, and Repeated Measures Models. Springer.

Wickham, H., & Stryjewski, L. (sub.). 40 years of boxplots. The American Statistican.

Wilkinson, G. N., & Rogers, C. E. (1973). Symbolic description of factorial models for analysis of variance. Applied Statistics, 22, 392–399.

Young, S. G., & Bowman, A. W. (1995). Nonparametric analysis of covariance. Biometrics, 51, 920–931.

Zar, J. H. (1998). Biostatistical Analysis. Pearson, Prentice Hall.