Ecrire un programme en R

On a vu à plusieurs reprises, dans le cadre des exercices de travaux pratiques, qu'il était parfois intéressant d'écrire des (petits) programmes pour effectuer de manière automatique certaines tâches répétitives. R offre également cette possibilité, qui est en réalité un autre atout majeur du logiciel. Nous allons illustrer cet aspect au travers d'exemples qui sont utiles dans le cadre des travaux pratiques.

I. Mélanger un tableau de données

A plusieurs reprises dans le cours, nous avons montré qu'il était possible d'obtenir la distribution sous l'hypothèse nulle d'une statistique d'intérêt. Ainsi, par exemple, si on effectue une régression linéaire, on calculera la valeur de b, estimateur du coefficient de régression linéaire β, puis on désirera tester l'hypothèse nulle H₀: β=0. Comme, en général, la valeur calculée de b (notée b_R) sera non nulle, se posera alors la question de savoir si obtenir une valeur aussi éloignée de 0 que celle obtenue arrive souvent par hasard (dans ce cas, on accepterait H₀) ou bien est un événement peu probable (au quel cas on rejetterait H₀). Si la distribution de b est connue, on peut répondre de manière directe à cette question en calculant sur cette distribution P(b > |b_R|), la probabilité de dépasser b_R (en valeur absolue). Une alternative est de générer la distribution de b à partir des données: si les valeurs de Y étaient mélangées, et que les couples (X,Y) utilisés dans la régression étaient reformés de manière aléatoire, détruisant ainsi l'éventuelle relation existant en réalité entre X et Y, les couples ainsi reformés pourraient être vus comme des couples générés à partir d'une situation où β=0, c'est-à-dire générés sous l'hypothèse nulle. Le b, noté b_P, calculé sur ces données permutées est donc un échantillon de la distribution de b recherchée. L'intérêt de cette approche est évidemment qu'il existe, si le nombre de couples (X,Y) est raisonnablement grand, un très grand nombre de possibilités de créer ces listes de couples mélangés (le nombre de possibilités est N!, où N est le nombre de couples. Par exemple, pour 10 couples, il y a 10! = 3628800 possibilités !). L'idée serait donc de mélanger un grand nombre de fois n le tableau des Y (par exemple), de générer ainsi n valeurs de b_P, puis de confronter b_R à la distribution de b_P obtenue pour vérifier si une valeur aussi grande (en valeur absolue) que celle obtenue est fréquente (P > α) ou pas (P < α). Nous commençons donc à décrire cette approche en fournissant une fonction permettant de mélanger un tableau, fonction qui pourra être invoquée pour la suite de la procédure décrite ci-dessus.

A. Le principe

Pour mélanger de manière aléatoire un tableau T, une manière simple de procéder est la suivante:

On tire au hasard un nombre (entier) entre 1 et N (la taille du tableau). Appelons x ce nombre.
On permute T[1] et T[x].
On tire au hasard un nombre (entier) entre 2 et N (la taille du tableau). Appelons y ce nombre.
On permute T[2] et T[y].
On continue ainsi jusqu'à l'élément (N-1) du tableau

Comme on souhaite invoquer cet algorithme de mélange de manière répétée, le plus simple est de créer une fonction R, qu'on pourrait appeler 'melange'. Si on veut mélanger le tableau Y, il suffira alors de taper dans R : Ym<-melange(Y), qui renverra dans Ym une version mélangée du tableau Y.

B. Le code

On veut donc écrire une fonction R qui implémente l'algorithme suivant:

pour chaque élément i entre 1 et (N-1)
tirer un nombre r au hasard entre i et N
échanger les éléments r et i

Le code est donc le suivant:

melange<-function(Y){
  Z<-Y
  n<-length(Z)
  for (i in 1:(n-1)) {
    r<-round(runif(1,min=i,max=n),digits=0)
    temp<-Z[r];Z[r]<-Z[i];Z[i]<-temp
  }
Z
  }

On peut ensuite tester le code avec l'exemple simple suivant:

Y<-c(1:10)
Ym<-melange(Y)
Ym

##  [1]  3  2  8  1 10  7  4  5  9  6

II. Effectuer une régression linéaire:

Nous allons à présent développer une autre fonction, que nous appelerons 'regression', qui recevra en arguments les vecteurs X et Y, et qui retournera les résultats de la régression sous une forme que nous allons décrire. Dans la version simplifiée que nous allons développer, les résultats seront les estimateurs a et b de α et β de l'équation de régression E[Y|X] = α + β * X, la valeur de la statistique F, la probabilité théorique (c'est-à-dire obtenue en utilisant la distribution théorique de F) de dépasser F et la probabilité empirique (c'est-à-dire obtenue en permutant le vecteur Y) de dépasser F. Si les données respectent les hypothèses de la régression, ces deux valeurs de probabilité devraient être proches. On pourrait bien entendu ajouter d'autres résultats à la fonction, mais dans cet exemple, on se limitera à ces 5 valeurs.

Une première question est: comment la fonction va-t-elle renvoyer les 5 valeurs ? Plusieurs options sont possibles, comme renvoyer ces 5 valeurs dans un tableau (par exemple, un tableau nommé RL) et d'accéder aux résultats en énumérant les éléments du tableau (RL[1] contient a, RL[2] contient b, etc...). Dans l'exemple, nous allons travailler avec une liste, ce qui rend plus explicite l'accés aux éléments: si la liste s'apelle RL, l'accès aux résultats se fait via 'RL$a', RL$b', 'RL$F', RL$pth', 'RL$pemp'. L'algorithme aura cette fois la structure suivante:

Calculer la moyenne pour X (Xm) et pour Y (Ym)
Calculer les écarts à la moyenne x = X - Xm et y = Y - Ym
Calculer les coefficients b = Σxy/Σx² et a = Ym - b * Xm
Calculer les sommes de carrés SCM = b*Σxy, SCT = Σy² et SCE = SCT - SCM
Calculer F = (SCM/1)/(SCE/n-2) = (n-2)*SCM/SCE
Calculer p(théorique) = p(>F)
Calculer p(empirique) par permutations. On répète 1000 fois:
- mélanger y
- on recalcule b(permuté) avec le nouveau vecteur d'y
- on comptabilise les valeurs de b(permuté) qui sont supérieures à b(réel).
On retourne la liste qui contient tous les résultats.

Voici le code qui implémente cet algorithme:

regression<-function(X,Y){
  #definition de la liste
  l<-list(a=0,b=0,F=0,pth=1,pemp=1)
  #calcul des moyennes
  mX<-mean(X)
  mY<-mean(Y)
  #calcul des écarts à la moyenne
  x<-X-mX
  y<-Y-mY
  #calcul de a et de b
  l$b<-sum(x*y)/sum(x*x)
  l$a<-mY-l$b*mX
  #calcul des sommes de carrés
  scm<-l$b*sum(x*y)
  sct<-sum(y*y)
  sce<-sct-scm
  #calcul de F et de p(>F) théorique
  l$F<-scm*(length(X)-2)/sce
  l$pth<-1-pf(l$F,1,length(X)-2)
  #calcul de p(>F) empirique
  l$pemp=0.0
  for (i in 1:1000){
    ys<-melange(y)
    bs<-sum(x*ys)/sum(x*x)
    if (abs(bs)>=abs(l$b))
    {l$pemp<-l$pemp+0.001}
  }
  l
}

On peut illustrer l'utilisation de cette fonction en générant 20 couples de la manière suivante:

X<-c(1:20)
Y<-rnorm(20, mean=2*X+3,sd=20)

Les valeurs de X varient donc de 1 à 20, et les valeurs de Y sont obtenues en tirant au hasard dans une distribution normale de moyenne 2*X+3 et de déviation standard égale à 20 (ces valeurs de Y sont prises comme exemple, pour simuler une situation où une relation existerait entre X et Y, soit Y = 2*X + 3, mais serait obscurcie par un bruit aléatoire gaussien qui écarte les points de la droite théorique). L'appel à la fonction de régression est alors très simple: soit on appelle directement la fonction et les 5 résultats sont affichés par R, soit on mémorise les résultats dans une liste, et on les affiche à la demande, comme montré ci-dessous:

RL<-regression(X,Y)
RL$a

## [1] -10.79611

RL$b

## [1] 3.291102

RL$F

## [1] 26.6051

RL$pth

## [1] 6.603348e-05

RL$pemp

## [1] 0

On remarquera que, dans cet exemple, malgré le bruit gaussien important, l'estimation de b n'est pas trop mauvaise (on obtient b = 3.291102 alors que β = 2). La régression est déclarée comme significative (p(>F) = 6.603348e-05). De plus, la valeur obtenue de manière empirique n'est pas très éloignée de la valeur théorique (on est, dans cet exemple, dans les conditions d'utilisation de la distribution théorique de F). Enfin, signalons que les résultats que vous obtiendriez en recopiant ce code diffèreront certainement légèrement de ce qui est montré ici, en raison du caractère aléatoire de la génération des valeurs de Y.

III. Effectuer une régression logistique

La régression logistique partage certaines caractéristiques avec la régression linéaire mais est compliquée par le fait qu'il n'existe pas de formule analytique pour obtenir les estimateurs de α et β, comme c'était le cas pour la régression linéaire. A la place, on est amené à utiliser un algorithme itératif qui, partant de valeurs de départ, notées a₀ et b₀, pour a et b, va faire évoluer ces valeurs progressivement vers des valeurs qui rendent la vraisemblance maximale (voir le cours théorique). Si on veut créer une fonction 'regression.logistique' de manière similaire à ce qui a été fait ci-dessus, on devra donc introduire l'algorithme en question dans la fonction.

Schématiquement:

On fournit des valeurs pour a₀ et b₀, et on calcule la valeur l₀ de la log-vraisemblance sur base de ces paramètres.
On obtient, via l'algorithme décrit plus bas, de nouvelles valeurs a₁ et b₁ telles que l₁>l₀.
On obtient, via l'algorithme décrit plus bas, de nouvelles valeurs a₂ et b₂ telles que l₂>l₁.
On répète cette procédure jusqu'à ce que, soit on ait atteint un nombre d'itérations fixés au préalable, soit il n'y ait plus de changement appréciable de la log-vraisemblance (l_i+1 ≈ l_i)

Un algorithme, nommé 'Newton-Raphson', a été proposé au cours théorique, pour passer des valeurs courantes de a et b (a_i et b_i) aux valeurs suivantes (a_i+1 et b_i+1):

Calculer: $\pi_i(X) = \frac{exp(a_i+b_i*X)}{1+exp(a_i+b_i*X)}$
Calculer: $p = \sum n_i*\pi_i*(1-\pi_i)$
Calculer: $q = \sum n_i*\pi_i*X_i*(1-\pi_i)$
Calculer: $r = \sum n_i*\pi_i*{X_i}^2*(1-\pi_i)$
Calculer: $s = \sum (n_i-r_i*\pi)$
Calculer: $t = \sum (n_i-r_i*\pi)*X_i$
Calculer: $a_{i+1} = a_i + (r*s-q*t)/(p*r-q^2)$
Calculer: $b_{i+1} = b_i + (p*t-q*s)/(p*r-q^2)$

Rappelons également que, à une constante près qui ne dépend ni de a ni de b, la log-vraisemblance se calcule via:

$l = \sum r_i*ln(\pi_i)+ \sum (n_i-r_i)*ln(1-\pi_i)$

Avec toutes les pièces du puzzle maintenant rassemblées, on peut passer à l'écriture de la fonction 'regression.logistique'. Cette fonction recevra en entrée un vecteur X et un vecteur Y, ce dernier ne contenant que des 0 et des 1. La fonction retournera une liste avec les valeurs estimées de a et b et un indicateur qui vaudra 0 si tout s'est déroulé correctement et 1 dans le cas contraire:

regression.logistique<-function(X,Y){
  #defintion de la liste
  l<-list(a=0,b=0,err=1)
  #valeur initiale de a et b
  a<-0
  b<-0
  #Boucle de calcul
  #On itère jusqu'à ce que: soit on ait atteint le nombre maximal d'itérations
  #soit la variation de vraisemblance est inférieure à un seuil donné
  maxiter<-100 #Nombre maximal d'itérations
  minseuil<-0.00001 #Seuil de variation de la vraisemblance
  ll<-1.0 #Initialisation du log-likelihood
  iter<-0 #Initialisation du nombre d'itérations
  seuil<-1.0 #Initialisation du seuil
  #Boucle while
  while ((iter<maxiter) && (seuil>minseuil)){
    pi<-exp(a+b*X)/(1+exp(a+b*X))
    p<-sum(pi*(1-pi))
    q<-sum(pi*(1-pi)*X)
    r<-sum(pi*(1-pi)*X*X)
    s<-sum(Y-pi)
    t<-sum((Y-pi)*X)
    newa<-a+(r*s-q*t)/(p*r-q*q)
    newb<-b+(p*t-q*s)/(p*r-q*q)
    newl<-sum(X*log(pi)+(1-X)*log(1-pi))
    iter<-iter+1
    seuil<-abs(ll-newl)
    a<-newa
    b<-newb
    ll<-newl
  }
  #Fin de la boucle
  #Si une solution a été trouvée, iter<maxiter
  #Sinon, une erreur est reportée
  l$a<-a
  l$b<-b
  l$err<-0
  if(iter>maxiter){ l$err<-1}
  #On retourne la liste
  l
}

On peut tester cette procédure en utilisant les données qui ont été proposées lors du cours théorique. C'est ce qui est fait dans le code qui suit:

X<-c(rep(34.5,8),rep(37.5,28),rep(40.5,41),rep(43.5,17),rep(46.5,4))
Y<-c(rep(1,1),rep(0,7),rep(1,6),rep(0,22),rep(1,18),rep(0,23),rep(1,9),rep(0,8),rep(1,3),rep(0,1))
regression.logistique(X,Y)

## $a
## [1] -10.38889
## 
## $b
## [1] 0.2463824
## 
## $err
## [1] 0

On obtient donc bien les mêmes résultats que ceux obtenus au cours théorique. Si on veut tester la signification de b, on peut écrire ce qui suit:

rlog<-regression.logistique(X,Y)
pemp<-0.0
np<-0
for (i in 1:1000){
  Ys<-melange(Y)
  temp<-regression.logistique(X,Ys)
  if (temp$err==0){
    np<-np+1
    if (abs(temp$b)>abs(rlog$b)) {pemp<-pemp+1.0}
  }
}
if (np>0){pemp<-pemp/np}
pemp

## [1] 0.007

On obtient une régression logistique très significative (p = 0.007): l'hypothèse nulle (β=0) est donc rejetée.