S4_programmation.Rnw

\chapter{Introduction à la programmation\label{s:progra}}

Dans les séances précédentes, nous avons utilisé des fichiers \texttt{.R} pour sauvegarder des listes d'instructions.
Nous avons aussi chargé et manipulé des jeux de données.
Il est souvent nécessaire d'automatiser tout ou partie de ce processus, ce qui implique de faire appel à de la programmation.

L'objectif de cette séance est de se familiariser avec les principaux concepts qu'on utilise pour concevoir un programme.
La première partie couvre les bases en algorithmie, c,-à-d. les boucles et les tests.
La deuxième partie concerne les fonctions, leur définition dans R, et leur utilisation.
\`A l'issue de cette séance, vous serez en mesure de vous attaquer a pratiquement tous les problèmes nécéssitant de manipuler des données.
Tout programme que vous aurez a écrire ne sera qu'une combinaison plus ou moins complexe des éléments abordés dans les séances précédentes. 

\section{Boucles}

Les boucles permettent parcourir une liste, ou de répéter une série d'instructions, dans des conditions bien définies; c'est une des structures de base de l'algorithmique.
R propose deux types de boucles, les boucles \texttt{for} et les boucles \texttt{while}.
En français, on peut les résumer par <<pour chaque>> et <<tant que>>. 

\subsection{Boucles de type \emph{for}}

Une boucle \texttt{for} permet de répéter un bloc d'instructions un nombre prédéfini de fois, ou d'éxécuter des commandes pour chaque élément d'un tableau de données.
La syntaxe de base est la suivante:

<<>>=
for (step in c(1:10)) cat(step)
@

En clair, pour chaque valeur entre 1 et 10, qu'on nomme \texttt{step} (mais qu'on aurait pû nommer n'importe comment), on affiche (\texttt{cat}) la valeur de step.
On peut bien sûr spécifier plusieurs instructions qui doivent être éxécutées à chaque \emph{itération} (étapes de la boucle) en utilisant les accolades:

<<>>=
for (step in c(1:3))
{
	cat(step)
	print(summary(rnorm(100,mean=step)))
}
@

Cette commande affiche le numéro de l'itération en cours (\texttt{cat(step)}), puis affiche les information de base (\texttt{summary}) sur une distribution normale (\texttt{rnorm}) centrée sur \texttt{step}.
Les boucles \texttt{for} peuvent contenir des instructions aussi longues que souhaité.

Une autre application des boucles \texttt{for} est de parcourir un objet.
Par exemple, on peut souhaiter, pour chaque élément d'un objet, afficher sa valeur.
R permet de réaliser ce genre d'opérations, avec la syntaxe suivante:

<<>>=
vect =  c('a','b','c','d')
for (val in vect) cat(val)
@ 

Pour chaque élément du vecteur \texttt{vect}, que l'on nomme \texttt{val} pour pouvoir y accéder pendant les itérations, R va afficher la valeur que l'élément contient.

Les boucles \texttt{for} sont donc une façon explicite de faire la même chose que ce qui a été traité dans la séance précédente, avec un inconvénient majeur:
elles sont extrèmement demandantes en temps de calcul.
Dans la majorité des cas, il est préférable d'avoir recours autant que possible aux fonctions \texttt{*apply}, qui sont bien plus optimisées que des boucles. 

\subsection{Boucles de type \emph{while}}

Les boucles de type \texttt{while}, littéralement \emph{pendant}, permettent de répéter une série d'instructions tant qu'une condition n'a pas été atteinte.
Pour cette raison, il faut bien prendre en compte le fait que mal utilisées, ces boucles peuvent ne jamais stopper.
Il faut donc faire particulièrement attention à la condition qui est évaluée à chaque itération.

Un exemple simple d'utilisation d'une boucle \texttt{while} est le calcul d'une factorielle.
On veut calculer $n!$, ce qui se fait simplement en multipliant l'ensemble des $1\leq k \leq n$.

<<>>=
n = 5
k = n
while(k > 1){
	k = k - 1
	n = n*k
}
print(n)
@ 

On remarquera que dans la parenthèse après \texttt{while} se trouve un test logique; les tests sont abordés dans la partie suivante.

\subsection{Sortir d'une boucle et sauter des étapes}

Lors de l'éxécution d'une boucle, on peut ne pas vouloir éxécuter toutes les instructions.
R possède des structures de contrôle pour effectuer ce type d'opérations.
Les principales qu'il faut connaître sont \texttt{break} et \texttt{next}.
L'instruction \texttt{break} permet de stopper l'éxécution de la boucle, c'est-à-dire de sortir de la boucle comme si la condition de sortie était remplie, ou le nombre maximal d'itérations atteint.
La boucle suivante va afficher les valeurs de \texttt{i} tant qu'elles sont inférieures à 3, et afficher \texttt{bye!} puis stopper la boucle sinon.

<<>>=
for(i in c(1:5)){
	if(i < 3){
		cat(i)
	} else {
		cat(' bye!')
		break
	}
}
@  

Une autre structure de contrôle intéréssant est \texttt{next}, qui permet de sauter une itération si la valeur de l'itérateur ne nous plaît pas.
Cette structure est particulièrement utile quand on réalise un nombre important d'opérations à chaque itération, et qu'on ne veut pas perdre de temps a traiter des valeurs qui ne nous intéréssent pas.
On peut utiliser \texttt{next} pour avoir, par exemple, une liste de tous les nombres pairs entre 1 et 10:

<<>>=
is.even <- function(x) x %% 2 == 0 
for(i in c(1:10)){
	if(is.even(i)){
		cat(i)
		cat(' is even\n')
	} else {
		next
	}
}
@

\section{Tests}

Cette partie est consacrée à l'utilisation des tests.
On a vu dans les séances précédentes l'existence de variables de type booléen, qui prennent les valeurs \texttt{TRUE} ou \texttt{FALSE}.
Des variables de ce type sont utilisées dans le cadre d'expressions conditionelles, c'est-à-dire quand on souhaite effectuer différentes instructions en fonction de la valeur d'une condition.

\subsection{Expressions conditionnelles}

La structure de base d'une expression conditionnelle est la suivante:

<<eval=FALSE>>=
if (condition){
	instruction(1)
} else {
	instruction(2)
}
@

Une notation comme

<<eval=FALSE>>=
if (condition) instruction
@

est aussi acceptable, de même que

<<eval=FALSE>>=
val = ifelse(condition,valeur_if,valeur_else)
@

Cette dernière notation permet de gagner du temps quand on veut que la valeur d'une variable dépende d'une condition.
Par exemple, $(n \geq 2)\lor(n = 0)$ s'écrit \texttt{(n <= 2) or (n == 0)}, ou encore \texttt{(n <= 2) | (n == 0)}.
Les différents opérateurs logiques sont regroupés dans le tableau~\ref{t:p:logical}.

\begin{table}[tb]
   \centering
   \begin{tabularx}{\textwidth}{lXlr}
\toprule
   \textbf{Opérateur} & \textbf{Signification} & \textbf{Version 1} & \textbf{Version 2}\\\midrule
$\lor$ & ou (au moins une des deux) & \texttt{or} & \texttt{|}, \texttt{||}\\
$\land$ & et (les deux) & \texttt{and} & \texttt{\&}, \texttt{\&\&}\\
$\lnot$ & non & \texttt{not} & \texttt{!}\\
$\oplus$ & ou conditionnel (seulement une des deux) & \texttt{xor(a, b)} & \\\midrule
$\in$ & est compris dans & \texttt{\%in\%} & \\\midrule
$=$ & égalité &  & \texttt{==} \\
$\leq$ & inférieur ou égal &  & \texttt{<=} \\
$\geq$ & supérieur ou égal &  & \texttt{>=} \\
$>$ & supérieur &  & \texttt{>} \\
$<$ & inférieur &  & \texttt{<} \\\bottomrule
\end{tabularx}
   \caption{Différents opérateurs logiques disponibbles dans R.}
   \label{t:p:logical}
\end{table}

L'argument \texttt{condition} prend la forme d'un test logique, qui peut être d'une complexité aussi grande que l'on veut.
Pour certains des opérateurs, il existe deux variantes (\texttt{|} et \texttt{||}, \texttt{\&} et \texttt{\&\&}).
On peut comprendre pourquoi avec les exemples suivants:

<<>>=
a = c(1,2,3,4,5,6)
b = c(1,2,3,4,4,6)
d = c(3,2,3,4,5,4)
(a == b)
(a == d)
(b == d)
(a == d) & (a == b)
(a == d) && (a == b)
(a == d) | (a == b)
(a == d) || (a == b)
all.equal(a,b)
@

L'opérateur répété une seule fois travaille par élément: les éléments des vecteurs sont comparés deux à deux (et on applique du recyclage).
L'opérateur répété deux fois ne prend en compte que les premières valeurs des vecteurs (un message d'avis sera émis si le vecteur a une taille supérieure à 1).

\subsection{Manipulation des valeurs booléennes}

Comme tous les autres types de données de R, les valeurs booléennes peuvent être manipulées dans des additions.

<<manipbool>>=
TRUE + TRUE
FALSE + TRUE
TRUE * TRUE
@

Cette possibilité de travailler sur les valeurs booléennes comme si \texttt{TRUE = 1} et \texttt{FALSE = 0} permet des raccourcis intéréssants.
Par exemple, on peut facilement connaître le nombre d'éléments d'un vecteur qui répondent à une condition particulière, ici $ x\in \left(-0.5; 1\right]$:

<<utilManipBool>>=
a = rnorm(100)
sum(a[(a<1)&(a>=0.5)])
sum((a<1)&(a>=0.5))
@

Dans l'exemple précédemt, la première ligne fait la somme de \texttt{a} aux indices satisfaisant la condition fixée.
La deuxième ligne fait la somme de ces indices, c.-à-d. d'un vecteur de booléens.

\section{Fonctions}

\subsection{Généralités}

L'utilisation des fonctions va permettre de gagner du temps dans la programmation.
Comprendre le principe des fonctions dépasse de beaucoup le cadre de R, et mérite qu'on s'y arrête.
Qu'est-ce qu'une fonction? Une série d'instructions qui vont, à partir d'arguments, renvoyer un résultat.
En quoi est-ce différent des scripts que nous avons utilisé jusqu'ici?
Écrire une fonction revient en quelque sorte a `expliquer' le code une fois, et R se charge ensuite de redonner la bonne valeur aux arguments.

Le parallèle le plus évident est celui des fonctions mathématiques: si $f(x) = x^2$, on peut calculer $f(x)$ pour tout $x$, parce qu'on sait quoi faire.
Ça devient donc très avantageux si on doit calculer $f(x)$ un grand nombre de fois.
En écrivant uniquement des scripts, pour calculer la valeur de beaucoup de $x^2$, on aurait écrit:

<<eval=FALSE>>=
1^2
2^2
3^2
4^2
5^2
@

ou encore

<<eval=FALSE>>=
for(i in c(1:5)) i^2
@

Si il faut revenir sur ce code plus tard, et transformer tous les \texttt{\^{}2} en \texttt{\^{}3}, la première solution implique de tout corriger manuellement.
En utilisant une fonction, la logique est différente:

<<>>=
f = function(x) x^2
f(2)
@ 

Peut importe le nombre de fois ou on devra effectuer l'opération contenue dans \texttt{f}, si on veut la modifier, elle sera toujours stockée au même endroit.
Dans la pratique, la majorité des instructions utilisées jusau'à présent sont des fonctions.
Par exemple, \texttt{print}, \texttt{read.table}, et \texttt{apply} sont des fonctions rendues disponibles par R.
Écrire une instruction comme \texttt{read.table('file.txt', h=TRUE)}, c'est faire appel à une fonction, en lui spécificiant certains \emph{arguments}.
L'objectif des parties suivantes est de comprendre comment les fonctions fonctionnent (!), afin de pouvoir en écrire de nouvelles. 

Les fonctions existent comme un espace <<à part>> dans R: ce qui se déroule dans une fonction, reste dans une fonction.
Prennons le cas du code suivant:

<<>>=
a = 2
b = 3
f = function(x){
	y = x+a
	z = y
	return(z)
}
f(b)
@

La fonction \texttt{f} peut aller chercher la valeur de \texttt{a} dans l'environnement global, mais tout ce qui est défini au sein de \texttt{f} est innaccessible.
D'ailleurs, tout ce qui est créé dans la fonction est détruit -- retiré de la mémoire -- une fois que la dernière instruction est éxécutée.
Voyez la section sur la fonction \texttt{return} pour plus de détails. Cette notion de quel objet est accessible est extrèmement importante à maîtriser.

Dans R, les objets existent dans deux <<mondes>>, l'environnement global, et l'intérieur de chaque fonction.
L'intérieur de chaque fonction peut avoir accès aux objets et aux variables de l'environnement global, même si cette pratique est à éviter pour différentes raisons (une variable globale peut être modifiée entre deux appels à la fonction, notamment).
En revanche, l'environnement global n'a pas accés aux objets crées ou modifiés à l'intérieur d'une fonction.
La communication avec les fonctions se fait  ...

Le fait qu'on puisse passer des objets d'un environnement à l'autre peut entraîner un comportement assez surprenant.
Dans l'exemple suivant:

<<>>=
a = 2
print(a)
f = function(a)
{
	a = a+1
	return(a)
}
print(f(a))
print(a)
@

l'appel à la fonction \texttt{f} devrait modifier \texttt{a}, puisque la seule instruction de cette fonction est \texttt{a = a+1}.
Or, quand on appelle cette fonction, puis qu'on affiche la valeur de \texttt{a}, elle n'a pas changé.
Ce comportement vient du fait que l'objet \texttt{a} qui existe dans la fonction n'est pas celui qui existe dans l'environnement global.
Par conséquent, le \texttt{a} de la fonction peut être modifiè, sans que cela n'affecte le \texttt{a} de l'environnement de travail.

Cet example illustre aussi pourquoi le nom des variables est important. On sait que R est capable d'aller chercher, depuis une fonction, des variables de l'environnement global.
Quelle version de \texttt{a} faut-il aller chercher?
Pour éviter les erreurs liées au fait que plusieurs variables aient le même nom, on essaie de donner un nom unique à toutes les variables.
C'est sans doute plus long à écrire -- même si cet argument n'est pas valable avec un éditeur qui auto-complète le code --, mais ça évite surtout les erreurs à l'éxécution.  

\subsection{Déclarer une fonction}

Comme illustré dans les exemples précédents, la déclaration d'une fonction se fait par

<<eval=FALSE>>=
nom_de_la_fonction = function(argument1,argument2)
{
	instructions
	return(sortie)
}
@

Les éléments les plus importants sont les arguments et l'instruction \texttt{return}.

\subsubsection{La commande \emph{return}}

La commande \texttt{return} est en général la dernière ligne d'une fonction: ce qui se passe après est ignoré.
Cette commande va renvoyer ce qui se trouve entre ses parenthèse dans l'environnement global de R: c'est un des points de communication entre l'environnement global et l'environnement de la fonction.
La commande \texttt{return} ne peut prendre qu'un seul argument, il faut donc regrouper les variables sous forme de liste ou autres si vous en avez plusieurs.

Prennons un exemple simple:

<<explFuncReturn>>=
somme_1 = function(a,b){
	S = a + b
}
somme_2 = function(a,b){
	S = a + b
	return(S)
}

a = somme_1(2,3)
a
b = somme_2(2,3)
b
@   

Quelques précisions.
On peut tout a fait nommer des variables \texttt{a} et \texttt{b} dans l'environnement global, même si les fonctions \texttt{somme\_1} et \texttt{somme\_2} utilisent des variables avec ces noms,
grâce au fait que le \emph{scope} de ces variables n'est pas le même.
Il est toutefois déconseillé de le faire.
Dans l'exemple précédent, la variable \texttt{S} n'existe pas, à aucun moment, dans l'environnement global de R.
Le seul moyen de récupérer ce qui a été renvoyé par les fonctions est de les assigner à une variable de l'environnement global.

Lorsque que la fonction ne se compose que d'une seule ligne, le résultat de cette dernière ligne est renvoyé et il n'y a pas besoin de spécifier \texttt{return}:

<<returnImplicit>>=
multipl = function(a,b) a*b
multipl(2,3)
@  

\subsubsection{Arguments}

Les arguments sont des variables que l'on peut échanger entre R et l'intérieur d'une fonction.
Dans l'exemple précédent, pour faire une somme, il faut additionner deux nombres: $f(a,b) = a + b$ s'écrit \texttt{f = function(a,b) a + b}.
Les éléments entre parenthèse sont les arguments de la fonction.
Typiquement, un argument est noté par \texttt{nom = valeur}, où \texttt{nom} est l'identifiant de cet argument \emph{dans la fonction}, et \texttt{valeur} est la valeur par défaut.

Les exemples suivants permettent de comprendre le fonctionnement des arguments et des valeurs par défaut.

<<argDefault>>=
somme = function(a=1,b=1) a + b
somme()
somme(a = 2)
somme(b = 3)
somme(1,2)
@ 


Il existe un argument particulièrement intéressant dans R: \texttt{...}.
Cet argument contient tout objet passé à une fonction, pour lequel il n'y pas pas de nom d'argument correspondant.
Par exemple:

<<ellipsisArg>>=
print_message = function(msg='Hello',...){
	print(msg)
	print(...)
}
print_message(msg='Hello','world!')
@ 

L'utilisation de cet argument \emph{catch all} demande un peu de réflexion et beaucoup de pratique,
mais s'avère particulièrement utile quand le format exact des arguments qui seront donnés a une fonction n'est pas connu d'avance.

\section{Mise en application}

\subsection{Test par permutation}

Lorsque les données devient de la normalité, on peut préférer réaliser un test paramétrétique avec des permutations plutôt qu'un test non paramétrique.
La plupart des programmes de statistique n'offrent pas cette possibilité qui demande pourtant très peu d'efforts pour être implémentée dans R.
Dans cette mise en application, on veut effectuer un test t, pour comparer deux distributions, disponibles dans un fichier \texttt{s4-data.txt}.

Le principe d'un test par permutations est simple.
La première étape est d'effectuer le test sur l'échantillon non permuté, pour obtenir la valeur de la statistique ($T$).
Dans le cas du test t, R propose la fonction \texttt{t.test}, et un rapide survol de \texttt{?t.test} vous donnera les arguments nécéssaires et la manière de récupérer la statistique. 
Commencent ensuite les permutations a proprement parler. Pour un nombre $n$ d'itérations choisies (en général 9999), on mélange l'ensemble des valeurs des deux distributions.
On reconstruit ensuite, en tirant au hasard dans le \emph{pool} de valeurs ainsi formées, deux distributions de taille égale.
Cette étape peut, par example, prendre la forme d'une fonction \texttt{resample}, qui prendrait une \texttt{data.frame} avec deux colonnes (la valeur, et le groupe d'origine) en argument.
R propose la fonction \texttt{sample}, qui permettra de mélanger la colonne correspondant au groupe (ce qui recréera automatiquement les deux distributions -- économisons nous!).
Une fois les deux distributions reconstruites, on calcule la nouvelle statistique $T'$.
Si la valeur de $T'$ est inférieure ou égale à la valeur de $T$, on incrémente une variable $N$ de 1.
Sinon, la valeur de $N$ reste la même.

Le calcul de la \emph{p-value} se fait de la manière suivante:

\begin{equation}
	p = \frac{N}{n+1}
\end{equation}

À partir de ces informations, et des informations données dans l'introduction de cette séance, vous devez être en mesure de programmer sans difficultés une fonction \texttt{t.test.permut}, qui permet de réaliser un test t par permutations.
En bonus, vous pouvez ajouter des arguments qui permettent de contrôler le nombre de permutations qui doivent être réalisées. 

\clearpage
\section{Solution des mises en application}

\subsection{Test par permutation}

On commence par écrire une fonction \texttt{resample}, qui mélange l'attribution des valeurs à un des deux groupes.

<<>>=
resample = function(df)
{
	df$group = sample(df$group)
	return(df)
}
@

Puis on écrit la fonction \texttt{t.test.permut}, qui effectue les permutations à proprement parler

<<>>=
t.test.permut = function(df,n=9999)
{
	baseStat = t.test(value~group,df)$statistic
	N = 0
	for(repl in c(1:n))
		if(t.test(value~group,resample(df))$statistic<=baseStat)
			N = N + 1
	return(N/(n+1))
}
@

On peut générer un jeu de données de test:

<<>>=
value = c(rnorm(100,mean=0),rnorm(100,mean=1))
group = c(rep('a',100),rep('b',100))
test_df = as.data.frame(cbind(value = value, group = group))
test_df$value = as.numeric(as.vector(test_df$value))
@

puis vérifier que la \emph{p-value} est inférieure à 0.05:

<<>>=
print(t.test.permut(test_df, n = 9))
@

Pour éviter les problèmes de lenteur des boucles \texttt{for}, on peut écrire la même fonction basée sur la fonction \texttt{replicate}:

<<>>=
t.test.permut = function(df,n=9999)
{
	baseStat = t.test(value~group,df)$statistic
	return(sum(replicate(n,
		t.test(value~group,resample(df))$statistic<=baseStat))/(n+1))
}
@

La différence sur 10000 réplicats est assez minime (3 secondes environ).
Mais dans un contexte ou il faut analyser plusieurs fois des jeux de données, ces petits écarts finissent par faire une différence importante.