Uncertainty and collective science

22 mai 2011

There are some things that you, the reader of this preface, know to be true, and others that you know to be false; yet, despite this extensive knowledge that you have, there remain many things whose truth or falsity is not known to you. We say that you are uncertain about them. You are uncertain, to varying degrees, about everything in the future; much of the past is hidden from you; and there is a lot of the present about which you do not have full information. Uncertainty is everywhere and you cannot escape from it. Truth and falsity are the subjects of logic, which has a long history going back at least to classical Greece. The object of this book is to tell you about work that has been done in the twentieth century about uncertainty.

[...]

Research is carried out by individuals and often the best research is the product of one person thinking deeply on their own. For example, relativity is essentially the result of Einstein’s thoughts. Yet, in a sense, the person is irrelevant, for most scientists feel that if he had not discovered relativity, then someone else would; that relativity is somehow ‘‘out there’’ waiting to be revealed, the revelation necessarily being made by human beings but not necessarily by that human being. This may not be true in the arts so that, for example, if Shakespeare had not written his plays it would not follow that someone else would have produced equivalent writing. Science is a collective activity, much more so than art, and although some scientists stand out from the rest, the character of science depends to only a very small extent on individuals and what little effect they have disappears over time as their work is absorbed into the work of others.

Dennis Lindley, Understanding uncertainty

Je ne sais pas pour vous, mais moi, quand je lis ça, j’ai envie de m’installer dans un fauteuil et lire toute la journée.

2 Commentaires | Citations, Science | Tagué: , , | Permalien
Publié par walrus

Le modèle linéaire (1) – présentation

11 octobre 2010

Motivation

Lorsque l’on étudie un phénomène naturel, on se retrouve souvent à mesurer quelque chose, la variable d’intérêt, par exemple la taille du géranium sur le balcon du voisin, et on aimerait bien savoir si cette chose (ici la taille du géranium) change en fonction d’une autre chose, par exemple selon qu’il est à la lumière ou à l’ombre (facteur qualitatif), ou bien selon la quantité d’eau que le voisin lui verse tous les jours (facteur quantitatif), etc.

Bien sur, en faisant une expérience et en regardant bien tous les jours, on aura peut-être la réponse à notre question. Mais pour être plus confiant dans nos conclusions, surtout quand le phénomène semble compliqué, c’est mieux de construire un modèle statistique et de tester une hypothèse, par exemple dans notre cas “le fait que le géranium soit au soleil ou à l’ombre influe sur sa taille”.

Dans l’exemple ci-dessus, notre variable d’intérêt est la taille du géranium et elle dépend de l’expérience: le géranium sera plus ou moins grand s’il y a eu plus ou moins de soleil le mois où on l’a planté. Initialement, on ne sait pas quelle taille le géranium aura au bout d’une semaine, deux semaines, etc.: la mesure de la taille peut prendre plusieurs valeurs selon le temps qu’il fait, c’est-à-dire selon le résultat de l’expérience. On parle alors de variable aléatoire. Soit Y la variable aléatoire correspondant à la taille du géranium, mesurée en centimètre. Elle suit une loi de probabilité, que l’on ne connaît pas a priori, mais à chaque évènement est associé une probabilité. A l’évènement “le géranium mesure 15 cm” est associé la probabilité P( Y = 15 ), idem pour toutes les tailles possibles.

Comme Y est une variable aléatoire, on peut représenter sa loi en représentant la distribution de ses valeurs, c’est-à-dire en traçant un histogramme. On peut aussi calculer son espérance (valeur moyenne de sa distribution) et sa variance (mesure de la dispersion de sa distribution).

Définition

Le modèle linéaire est le modèle statistique le plus simple cherchant à expliquer une variable observée à l’aide de variables explicatives. Il s’écrit de la manière suivante:

y_{tj} = m_t + e_{tj}

  • les y_{tj} sont les valeurs prises par notre variable aléatoire Y lors de n expériences, ce sont les observations;
  • t est l’indice d’un traitement, un traitement étant une combinaison de niveaux de facteurs qualitatifs et quantitatifs utilisés lors des n expériences (un traitement est donc fixe, non aléatoire);
  • j est un indice de répétition;
  • m_t est l’espérance de y_{tj};
  • les e_{tj} sont les résidus entre les valeurs observées et ce qui est expliqué par m_t, correspondant à la variabilité du matériel expérimental, celle due aux facteurs non contrôlés, aux erreurs de mesure et à la randomisation.

Le modèle linéaire repose sur deux grandes hypothèses:

  • la partie fixe du modèle, m_t, est linéaire;
  • la partie aléatoire du modèle, e_{tj}, suit une loi normale d’espérance 0 et de variance \sigma^2.

Ceci suppose que, si l’on modélise un phénomène naturel avec ce modèle-là, il faut vérifier que les hypothèses sont vérifiées, sinon on ne peut pas tirer de conclusion.

Le modèle linéaire a une écriture plus générale sous forme matricielle:

y = X \Theta + e

  • y est un vecteur de dimension n (matrice colonne à n lignes) dans lequel on range les observations (l’ordre est arbitraire);
  • e est aussi un vecteur de dimension n dans lequel on range les variables résiduelles du modèle (même ordre que pour y bien sûr…);
  • \Theta est une matrice à p lignes, contenant les p paramètres;
  • X est une matrice à n lignes et p colonnes contenant les valeurs des facteurs qualitatifs et/ou quantitatifs, c’est la matrice du plan d’expérience.

Avant d’aller plus loin, voyons ce que signifie “plan d’expérience” et “randomisation“.

Prenons un cas bateau où y correspond au rendement d’une variété de blé que l’on cherche à expliquer à l’aide de k facteurs: y = f( x_1, x_2, ..., x_k). Selon la méthode classique, on fixe tous les facteurs sauf un que l’on fait varier, et ainsi de suite pour chaque facteur. Si k = 7 et que l’on veuille 5 répétitions de l’expérience afin d’avoir une bonne précision de mesure, on a 5^7 = 78125 expériences à faire: beaucoup trop… Grâce à la méthode du plan d’expérience, on fait varier les niveaux de tous les facteurs à la fois à chaque expérience, ce qui permet d’explorer une bonne partie des combinaisons sans pour autant réaliser toutes les expériences possibles.

Supposons maintenant que l’on ait un variable V_1 à deux niveaux A, B et une variable V_2 à deux niveaux C, D. Dans quel ordre doit on faire l’expérience ? V_1 = A puis V_1 = B puis V_2 = C puis V_2 = D ? Si on fait ça on risque une distribution non aléatoire des erreurs; pour l’éviter, on tire au sort. La randomisation consiste à déterminer aléatoirement l’ordre des expériences lors d’une étude statistique avec plan d’expérience. Cela augmente un peu l’erreur expérimentale mais diminue les biais.

Exemples

Voici maintenant quelques exemples pour illustrer ces notions.

Exemple 1: plan complètement randomisé à un facteur qualitatif

y: rendement d’une variété de blé

facteur qualitatif: mode de semis (par exemple labour ou semis direct)

randomisation: trois répétitions sur deux parcelles chacune, avec tirage au sort

Pour un exemple de ce type, on utilise le modèle dit d’analyse de la variance (ANOVA):

y_{ij} = \mu_{i} + e_{ij} ou bien y_{ij} = \mu + \alpha_{i} + e_{ij}

avec i = 1, ..., I=2 et j = 1, ..., J=3

Cela revient à écrire E(y_{ij}) = \mu_{i} avec le premier paramétrage et E(y_{ij}) = \mu + \alpha_{i} avec le deuxième.

On peut remarquer que le modèle est indéterminé dans le deuxième paramétrage puisque pour une valeur \lambda quelconque:

\mu + \alpha_i = ( \mu + \lambda) + (\alpha_i - \lambda) = \mu' + \alpha_i'

Ainsi, pour déterminer les paramètres et leur donner un sens, il faudra rajouter une contrainte. Mais par contre, la différence \alpha_i - \alpha_i' ne dépend pas de la contrainte choisie.

En écriture matricielle, dans le cas du premier paramétrage ça donne:

y = X_{(1)} \Theta_{(1)} + e

\begin{pmatrix}y_{11}\\ y_{12}\\ y_{13}\\ y_{21}\\ y_{22}\\ y_{23}\end{pmatrix} = \begin{pmatrix}1&0\\1&0\\1&0\\ 0&1\\ 0&1\\ 0&1\end{pmatrix} \begin{pmatrix}\mu_1\\\mu_2\end{pmatrix} + \begin{pmatrix}e_{11}\\e_{12}\\e_{13}\\e_{21}\\e_{22}\\e_{23}\end{pmatrix}

ce qui équivaut au système d’équations y_{11} = \mu_1 + e_{11}, y_{12} = \mu_1 + e_{12}, et ainsi de suite

Et dans le cas du deuxième paramétrage:

y = X_{(2)} \Theta_{(2)} + e

\begin{pmatrix}y_{11}\\ y_{12}\\ y_{13}\\ y_{21}\\ y_{22}\\ y_{23}\end{pmatrix} = \begin{pmatrix}1&1&0\\1&1&0\\1&1&0\\1&0&1\\1&0&1\\1&0&1\end{pmatrix} \begin{pmatrix}\mu\\\alpha_1\\\alpha_2\end{pmatrix} + \begin{pmatrix}e_{11}\\e_{12}\\e_{13}\\e_{21}\\e_{22}\\e_{23}\end{pmatrix}

Dans les équations concises ci-dessus, il est nécessaire de distinguer X_{(1)} et X_{(2)} ainsi que \Theta_{(1)} et \Theta_{(2)} pour montrer qu’elles sont différentes, mais pas e puisque cette matrice reste inchangée.

Exemple 2: plan complètement randomisé à un facteur quantitatif

y: nombre d’épis par pied d’une certaine variété de blé

facteur quantitatif: dose d’engrais azoté, en kg

randomisation: on attribue au hasard chacune des I valeurs x_i de doses choisies à chacune des I parcelles expérimentales

On peut choisir un modèle de régression simple:

y_i = a + bx_i + e_{i(1)} avec i = 1, .., I = 5

ou bien un modèle de régression polynomiale par exemple de degré 2:

y_i = a_0 + a_1x_i + a_2x_i^2+ e_{i(2)} avec i = 1, .., I = 5

Dans le premier cas, E(y_i) = a + bx_i et dans le deuxième E(y_i) = a_0 + a_1x_1 + a_2x_2.

En écriture matricielle pour le premier modèle:

y = X_{(1)} \Theta_{(1)} + e_{(1)}

\begin{pmatrix}y_1\\y_2\\y_3\\y_4\\y_5\end{pmatrix} = \begin{pmatrix}1&x_1\\1&x_2\\1&x_3\\1&x_4\\1&x_5\end{pmatrix} \begin{pmatrix}a\\b\end{pmatrix} + \begin{pmatrix}e_{1(1)}\\e_{2(1)}\\e_{3(1)}\\e_{4(1)}\\e_{5(1)}\end{pmatrix}

Et pour le deuxième modèle:

y = X_{(2)} \Theta_{(2)} + e_{(2)}

\begin{pmatrix}y_1\\y_2\\y_3\\y_4\\y_5\end{pmatrix} = \begin{pmatrix}1&x_1&x_1^2\\1&x_2&x_2^2\\1&x_3&x_3^2\\1&x_4&x_4^2\\1&x_5&x_5^2\end{pmatrix} \begin{pmatrix}a_0\\a_1\\a_2\end{pmatrix} + \begin{pmatrix}e_{1(2)}\\e_{2(2)}\\e_{3(2)}\\e_{4(2)}\\e_{5(2)}\end{pmatrix}

Ici on distingue aussi e_{(1)} et e_{(2)} car on change de modèle et pas seulement de paramétrage. On peut d’ailleurs s’attendre à ce que la variance de e_{(2)} soit un plus faible que celle de e_{(1)} car on a rajouté dans l’espérance E(y_i) l’effet du carré de x_i.

Exemple 3: plan en blocs complets randomisés à un facteur qualitatif.

à faire

Exemple 4: plan en blocs incomplets randomisés à un facteur qualitatif.

à faire

Exemple 5: plan en blocs complets randomisés à un facteur quantitatif.

à faire

Exemple 6: plan factoriel à deux facteurs qualitatifs.

à faire

Exemple 7: plan factoriel à deux facteurs quantitatifs.

à faire

Exemple 8: plan factoriel à un facteurs qualitatif et un facteur quantitatif.

à faire

Exemple 9: plan factoriel à deux facteurs hiérarchisés.

à faire

Suite: “Le modèle linéaire – estimation des paramètres (2)”

Source: polycopié “Le Modèle Linéaire” de Camille Duby, mai 2003, INA P-G + quelques ajouts

1 Commentaire | Science | Tagué: , | Permalien
Publié par walrus

Qu’est-ce que l’indépendance conditionnelle ?

14 septembre 2010

L’autre jour, je lisais un bouquin de stats et je suis tombé sur l’exemple ci-dessous expliquant ce qu’est l’indépendance conditionnelle (par rapport à l’indépendance tout court bien sûr). Je l’ai donc recopié ici afin de ne pas l’oublier.

Le problème. Considérons deux usines, A et B, fabriquant des montres. L’usine A fabrique en moyenne une montre cassée sur 100, et la B une sur 200. Maintenant supposons qu’un revendeur reçoive une livraison de montres sans savoir de quelle usine elles proviennent. Il regarde la première montre: elle marche. Quelle est la probabilité que la seconde marche aussi ?

La formalisation. Soit X_n l’état de la n-ième montre livrée, avec X_n = 1 si elle marche et X_n = 0 sinon. Soit Y l’usine d’origine. A priori on n’a aucun indice sur l’usine de provenance, donc on pose l’hypothèse suivante:

P( Y = A ) = P( Y = B ) = \frac{1}{2}

On suppose aussi que, sachant Y = A, les états des montres testées successivement sont indépendants (idem pour Y = B):

P( X_1 = 1, X_2 = 0 | Y = A ) = P( X_1 = 1 | Y = A ) \times \\ P( X_2 = 0 | Y = A )

Par contre, on sait quand même que P( X_n = 0 | Y = A ) = 0.01 et que P( X_n = 0 | Y = B ) = 0.005.

La reformulation de la question. On veut ici calculer P( X_2 = 1 | X_1 = 1 ).

La solution. D’après la formule de Bayes, on peut écrire:

P( X_2 = 1 | X_1 = 1 ) = P( X_1 = 1, X_2 = 1 ) / P( X_1 = 1 )

Or on peut aussi dire que le numérateur vaut:

P( X_1 = 1, X_2 = 1 ) = P( X_1 = 1, X_2 = 1 | Y = A ) \times \\ P( Y = A ) + P( X_1 = 1, X_2 = 1 | Y = B ) \times P( Y = B )

En utilisant le fait que X_1 et X_2 sont indépendants conditionnellement à Y:

P( X_1 = 1, X_2 = 1 ) = P( X_1 = 1 | Y = A ) \times P( X_2 = 1 | Y = A ) \times P( Y = A ) + P( X_1 = 1 | Y = B ) \times P( X_2 = 1 | Y = B ) \times P( Y = B )

Puis en remplaçant par les valeurs numériques:

P( X_1 = 1, X_2 = 1 ) = \frac{99}{100} \times \frac{99}{100} \times \frac{1}{2} + \frac{199}{200} \times \frac{199}{200} \times \frac{1}{2}

On aboutit à:

P( X_1 = 1, X_2 = 1 ) = \frac{1}{2} \times ( \frac{99}{100}^2 + \frac{199}{200}^2 )

De même, le dénominateur vaut:

P( X_1 = 1 ) = P( X_1 = 1 | Y = A ) \times P( Y = A ) + \\ P( X_1 = 1 | Y = B ) \times P( Y = B )

P( X_1 = 1 ) = \frac{1}{2} \times ( \frac{99}{100} + \frac{199}{200} )

En mettant les deux ensembles, on obtient:

P( X_2 = 1 | X_1 = 1 ) = \frac{ \frac{99}{100}^2 + \frac{199}{200}^2 }{ \frac{99}{100} + \frac{199}{200} }

P( X_2 = 1 | X_1 = 1 ) = 0.9925063

La différence entre “indépendant” et “conditionnellement indépendant”. Il faut bien noter que les états des deux montres ne sont pas indépendants “tout court”, mais indépendants conditionnellement à leur usine de provenance… Et oui, si les états des deux montres étaient inconditionnellement indépendants, on aurait:

P( X_2 = 1 | X_1 = 1 ) = P( X_2 = 1 )

Et de là:

P( X_2 = 1 ) = P( X_2 = 1 | Y = A ) \times P( Y = A ) + P( X_2 = 1 | Y = B ) \times P( Y = B )

Soit encore:

P( X_2 = 1 ) = \frac{1}{2} \times ( \frac{99}{100} + \frac{199}{200} )

P( X_2 = 1 ) = 0.9925

Alors oui c’est proche, mais ce n’est pas égal…

Laisser un Commentaire » | Science | Tagué: , | Permalien
Publié par walrus

Utiliser R pour placer des lycées sur la carte d’Ile-de-France

12 septembre 2010

Prenons un exemple au hasard, tiens, une liste de lycées d’Ile-de-France. Je veux apprendre à utiliser R pour tracer une carte de l’Ile-de-France et y placer dessus les lycées de ma liste.

Pourquoi R ? parce que c’est un environnement informatique de calcul statistique disposant de fonctionnalités graphiques, et puis c’est libre et gratuit, et ça fonctionne sur Windows, Linux et Mac OS. Donc tout lecteur de ce billet peut tester sur son propre ordinateur si je ne raconte pas n’importe quoi.

Tout d’abord, il me faut récupérer les données géographiques de ma carte. Pour cela, je vais sur le site GADM (global administrative areas) qui met gracieusement à disposition une base de données contenant les coordonnées spatiales des régions administratives du monde entier. Puis je télécharge les données correspondant à la France, niveau 2, au format “R data”.

Avant d’aller plus loin, je dois installer plusieurs librairies R spécialisées dans les coordonnées spatiales: “sp”, “maptools”, gpclib” et l’excellent “ggplot2″. La procédure d’installation est simple: ouvrez R et utilisez la commande install.packages( “sp” ). Pour ceux qui utilisent Ubuntu, l’installation de “maptools” directement dans R utilise trop de mémoire (voir ce thread). Il faut donc l’installer via la commande suivante:

$ R CMD INSTALL --no-latex package.tar.gz

Maintenant, je dois apprendre à manipuler les coordonnées spatiales dans R afin de tracer la carte en question. Pour cela, j’ouvre R, je charge les librairies dont j’ai besoin, je charge les données spatiales, je les convertis et je les trace:

library(ggplot2)
library(sp)
library(maptools)
library(gpclib)
gpclibPermit()
load("FRA_adm2.RData")
gadm_idf <- gadm[ gadm$NAME_1 == "\xcele-de-France", ]
idf <- fortify( gadm_idf, region="ID_2" )
qplot( long, lat, group=group, data=idf, geom="polygon", colour=I("white") )

Et voilà, on obtient le graphique ci-dessous:

Voilà une bonne chose de faite. A présent, je dois récupérer les coordonnées géographiques exactes des lycées qui m’intéressent (longitude et latitude). Heureusement Google Maps permet de faire ça mais il faut au préalable activer cette fonctionnalité (en savoir plus ici). Je reconnais que c’est un peu long et fastidieux, mais à la fin vous obtenez quelque chose comme ça:

lycée Jean Rostand,Mantes-la-Jolie,48.998052,1.696765

lycée Romain Rolland,Ivry-sur-Seine,48.800683,2.392895

etc…

Enfin, il ne me reste plus qu’à ajouter un point sur la carte à l’emplacement de chaque lycée. Pour cela rien de plus facile, commençons par ajouter le lycée Jean Rostand dont les coordonnées sont ci-dessus:


qplot( long, lat, group=group, data=idf, geom="polygon", colour=I("white") )
last_plot() + geom_point( aes(x, y, group=1), colour="red", data=data.frame(x=1.696765, y=48.998052) )

Pour ajouter tous les lycées d’un coup et placer des points de taille proportionnelle au nombre d’élèves dans ces lycées, voici la marche à suivre. Les données sont maintenant prises au hasard, à titre d’exemple. La fonction “theme_whiteMap()” n’est utilisée qu’à titre esthétique.


n <- 80
lycees 0, 1, 0))
lycees$bac[lycees$bac < 0.3] <- NA
lycees <- lycees[ lycees$lat1.5, ]
theme_whiteMap <- function(base_size = 12)
{
structure(list(axis.line = theme_blank(),
axis.text.x = theme_blank(), axis.text.y = theme_blank(), axis.ticks = theme_blank(),
axis.title.x = theme_blank(), axis.title.y = theme_blank(),
axis.ticks.length = unit(0.3, "lines"), axis.ticks.margin = unit(0.5, "lines"),
legend.background = theme_rect(colour = NA),
legend.key = theme_rect(colour = "grey80"), legend.key.size = unit(1.2, "lines"),
legend.text = theme_text(size = base_size * 0.8),
legend.title = theme_text(size = base_size * 0.8, face = "bold", hjust = 0),
legend.position = "right",
panel.background = theme_rect(fill = "white", colour = NA),
panel.border = theme_blank(),
panel.grid.major = theme_blank(),
panel.grid.minor = theme_blank(),
panel.margin = unit(0.25, "lines"),
strip.background = theme_blank(), strip.text.x = theme_blank(), strip.text.y = theme_blank(),
plot.background = theme_blank(),
plot.title = theme_text(size = base_size * 1.2, face = "bold"), plot.margin = unit(c(1,
1, 0.5, 0.5), "lines")), class = "options")
}
png( "ex_lycees_idf.png", width=600, height=500 )
pl <- qplot(long, lat, data = lycees, geom = "point")
pl <- pl + geom_polygon(aes(group=group), data=idf, colour=I("lightcyan4"), fill=I("transparent"))
pl <- pl + geom_point(aes(colour = bac, shape = factor(internat), size = eleves))
pl <- pl + theme_whiteMap()
pl <- pl + opts(title = "Exemples de lycées en Ile-de-France")
pl <- pl + labs(colour = "Réussite au bac", shape = "Internat", size = "Nombre d'élèves")
print(pl)
dev.off()

Au final, on obtient une carte avec plein d’infos dessus:

Merci à TJ sans qui j’étais bien incapable de faire tout ça…

Les observateurs auront remarqué que j’ai publié ce billet dans les catégories “Politique” et “Science”. La première parce que les données sont fortement inspirées des activités d’une association “œuvrant pour la cité”. La deuxième parce que l’acquisition de données, leur analyse et leur visualisation sont au cœur du travail quotidien des scientifiques.

Laisser un Commentaire » | Politique, Science | Tagué: , , | Permalien
Publié par walrus

Modéliser l’évolution des séquences d’ADN

28 mars 2010

Bien que je m’empêche de réagir à chaud, l’actualité médiatique est une source permanente d’inspiration. Alors quand les mathématiciens avertissent du desintérêt croissant pour leur discipline et que l’apparition de virus HINI mutants crée des frayeurs, l’envie me prend de dire quelques mots mêlant modèle mathématique et mutations dans des séquences d’ADN.

Attention le billet est long mais à la fin puisse votre persévérance être récompensée. Non seulement vous connaîtrez une formule estimant la distance entre deux séquences d’ADN, mais vous serez aussi passés par toutes les étapes du calcul (ce qui est généralement négligé partout ailleurs sur internet). Et puis le fait de rentrer dans le coeur d’un modèle, suivant pas à pas la démarche du modélisateur, pourrait réconcilier certains lecteurs avec les maths ? En tout cas moi ça m’a bien servi d’écrire ce billet.

On l’a vu précédemment, l’ADN est la molécule support de l’information génétique. Elle est formée de quatre “briques” appelées nucléotides, chacune symbolisée par une lettre: A, T, G et C. Chaque cellule possède une molécule d’ADN qui est copiée dans son intégralité au moment de la division cellulaire. Or, lors de cette réplication, des erreurs peuvent être commises, par exemple un A remplacé par un G. On parle alors de mutation (de substitution pour être plus précis). Par simplicité on ne traite pas le cas des insertions/délétions.

Depuis les années 1970, avec la possibilité de séquencer de l’ADN, on peut observer un fragment du génome d’un organisme et le comparer avec le même fragment mais provenant d’un autre organisme. En alignant ces deux séquences, on peut avoir une idée du nombre de match (:) et mismatch (*) les séparant, comme le montre l’image ci-dessous:

figure 1: un alignement global entre deux séquences d’ADN

En voyant cela, on se rappelle qu’une séquence d’ADN est une chaîne de caractères de longueur finie. Chaque position est appelée un site. Chaque site est dans un état particulier parmi quatre possibles (A, T, G ou C). Lorsque plusieurs séquences ont un ancêtre commun, on parle de séquence homologues. Sur la figure 1 la séquence S_1 et la séquence S_2 ont la même longueur L, 20 nucléotides, l’alignement résultant a donc 20 sites.

Pour mesurer la distance entre deux séquences, le plus simples est de calculer la proportion de sites différents. Sur la figure 1 une distance de 0.1 sépare la séquence S_1 de la séquence S_2 (2 mismatches sur 20 sites). Le problème c’est qu’en faisant cela on néglige les mutations cachées. Par exemple au site n°2 supposons que l’ancêtre commun ait été le nucléotide A. La séquence S_1 n’a pas eu de mutation mais la séquence S_2 aurait pu en avoir deux, d’abord de A vers C puis de C vers A. Or avec la distance décrite ci-dessus on ne compte pas de mutation puisque les deux séquences sont dans l’état A au 2e site. Dans la même veine, au 8e site supposons que la séquence ancestrale ait été dans l’état G. Aujourd’hui on observe que la séquence S_1 est toujours dans l’état G mais la séquence S_2 est dans l’état A, on compte donc une mutation, mais qui nous dit qu’il n’y en a pas eu plusieurs, par exemple de G vers T puis de T vers A ? Alors allons-y, écrivons un modèle mathématique prenant cela en compte !

On suppose que tous les sites suivent la même distribution de probabilités et que chaque site évolue de façon indépendante des autres. Ainsi la probabilité de passer de la séquence S_1 à la séquence S_2 (toutes les deux de longueur L) est donnée par:

\mathbb{P}( S_1 \rightarrow S_2 ) = \displaystyle \prod_{i = 1} ^L \mathbb{P}( S_1[i] \rightarrow S_2[i] )

Dire que deux probabilités sont indépendantes revient à les multiplier. Comme on suppose que tous les suites suivent la même loi de probabilité, on peut se concentrer sur l’évolution d’un seul site. Grâce à la formule ci-dessus, si on arrive à modéliser \mathbb{P}( S_1[i] \rightarrow S_2[i] ) on arrivera à modéliser \mathbb{P}( S_1 \rightarrow S_2 ).

Supposons maintenant que le temps avance en “tic”, comme les aiguilles d’une horloge, et qu’à chaque “tic” une mutation peut subvenir ou non. Pour modéliser cela on va utiliser une chaîne de Markov. Une chaine de Markov est un processus stochastique (synonyme de probabiliste). Par exemple, si on note X l’évènement “le dé jeté affiche la valeur x“, on dit que X est une variable aléatoire prenant les valeurs 1, 2 ... 6. Et bien, si on étudie plusieurs lancés de dé au cours du temps, on se retrouve à étudier un processus probabiliste: pas compliqué…

Si maintenant on considère que la probabilité au temps t_{n+1} ne dépend que de l’état présent, c’est-à-dire du temps t_{n}, et non des états passés, les temps t_{n-1}, t_{n-2}, etc, on dit que le processus possède la propriété de Markov. Résumé d’une autre façon: “le futur ne dépend du passé qu’au travers de l’instant présent”. Et c’est tout naturellement ce qui arrive en génétique: une mutation à la génération des petit-enfants ne va pas dépendre du nucléotide en question à la générations des grand-parents mais uniquement du nucléotide en question à la génération des parents.

Une chaîne de Markov est caractérisée par sa matrice de transition P. Quand on modélise l’évolution d’une séquence d’ADN, la chaîne a 4 états (pour A, T, G et C), et la matrice aura 4 lignes et 4 colonnes. La valeur au croisement de la ligne i et de la colonne j est la probabilité p_{ij} d’être dans l’état i et de passer dans l’état j.

En 1969, Jukes et Cantor proposaient de modéliser l’évolution d’une séquence via une chaîne de Markov sous l’hypothèse que la probabilité f(t) de passer d’un nucléotide à un autre pendant la durée t était constante au cours du temps. Voici la matrice de transition correspondante:

P = \begin{pmatrix}1 - 3f(t) & f(t) & f(t) & f(t) \\ f(t) & 1 - 3f(t) & f(t) & f(t) \\ f(t) & f(t) & 1 - 3f(t) & f(t) \\ f(t) & f(t) & f(t) & 1 - 3f(t) \end{pmatrix}

Je précise de manière arbitraire que les 1e ligne et colonne correspondent au nucléotide A, les 2e au nucléotide T, les 3e à G et les 4e à C: p_{23} à l’intersection de la 2e ligne et de la 3e colonne est la probabilité que le site soit dans l’état T et mute vers l’état G.

i \ne j: \mathbb{P}( i \rightarrow j ) = p_{ij}(t) =f(t)

i = j: \mathbb{P}( i \rightarrow j ) = p_{ii}(t) = 1 - 3f(t)

Si on est à la génération t et que le nucléotide que l’on est en train d’analyser est un “A”, deux scénarios sont possibles:

  • il y a une mutation, de A vers T, G ou C, chaque évènement ayant une probability f(t) d’arriver;
  • il n’y a pas de mutation, ceci avec une probabilité 1 - 3f(t) (la somme de tous les évènements doit faire 1).

Maintenant calculons f(t) = p_{ij}(t) et pour ceci commençons par différencier la matrice de transition P, c’est-à-dire regardons ce que vaut cette matrice P à l’instant t + \Delta t (c’est-à-dire très peu de temps après l’instant t) on a:

P( t + \Delta t ) = P(t) \times P(\Delta t)

Donc si maintenant on fait tendre \Delta t vers 0 (ça vous rappelle la définition de la dérivée n’est-ce pas ?):

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P( t + \Delta t ) - P(t)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P(t) \times P(\Delta t) - P(t + 0)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P(t) \times P(\Delta t) - P(t) \times P(0)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} P(t) \frac{P(\Delta t) - P(0)}{\Delta t}

P'(t) = P(t) \lim_{\Delta t \rightarrow 0} \frac{P(\Delta t) - P(0)}{\Delta t}

Ainsi: P'(t) = P(t) \times P'(0)

P'(0) = \begin{pmatrix} - 3f'(0) & f'(0) & f'(0) & f'(0) \\ f'(0) & - 3f'(0) & f'(0) & f'(0) \\ f'(0) & f'(0) & - 3f'(0) & f'(0) \\ f'(0) & f'(0) & f'(0) & - 3f'(0) \end{pmatrix}

Posons f'(0) = \alpha, on a alors:

P'(0) = \begin{pmatrix} - 3\alpha & \alpha & \alpha & \alpha \\ \alpha & - 3\alpha & \alpha & \alpha \\ \alpha & \alpha & - 3\alpha & \alpha \\ \alpha & \alpha & \alpha & - 3\alpha \end{pmatrix}

En multipliant, par exemple, la 1e ligne de P(t) avec la 2e colonne de P'(0) on obtient:

p'_{12}(t) = \alpha - 3\alpha f(t) - 3\alpha f(t) + \alpha f(t) + \alpha f(t)

p'_{12}(t) = \alpha - 4 \alpha f(t)

Or on sait aussi que si i \ne j, on a p'_{12}(t) = p'_{ij}(t) = f'(t).

C’est-à-dire: f'(t) = \frac{df}{dt} = \alpha - 4 \alpha f(t)

Et maintenant on intègre cette équation différentielle:

\int \frac{df}{\alpha - 4 \alpha f(t)} = \int dt

\Rightarrow \frac{ln(\alpha - 4 \alpha f(t))}{-4\alpha} = t + c

\Rightarrow ln( \alpha - 4 \alpha f(t) ) = - 4 \alpha t + c

\Rightarrow \alpha - 4 \alpha f(t) = c \exp^{- 4 \alpha t}

\Rightarrow f(t) = \frac{1}{4} - \frac{c}{4 \alpha} \exp^{- 4 \alpha t}

Il nous faut maintenant calculer c qui est la constance d’intégration. Pour cela supposons que f(0) = 0 ce qui signifie qu’au temps t=0 on commence dans un état constant. Par exemple, si à t=0 on est dans l’état A alors la probabilité d’avoir une substitution de ce “A” à t=0 vaut 0. Ainsi:

f(0)=0 \Rightarrow f(0) = \frac{1}{4} - \frac{c}{4\alpha} = 0 \Rightarrow c = \alpha

On a donc: f(t) = \frac{1}{4} - \frac{\exp^{- 4 \alpha t}}{4}

Pour calculer \alpha on se place à t=0 puisque \alpha = f'(0), et on appelle \Pi_i la probabilité d’être dans l’état i au temps t=0 (\Pi_i = 1/4):

\sum_{i} \sum_{j \ne i} \Pi_{i} P'(0)_{ij} = 1

\Pi_{A} P'(0)_{AT} + \Pi_{A} P'(0)_{AG} + \Pi_{A} P'(0)_{AC} + \Pi_{T} P'(0)_{TA} + ... = 1

12 \frac{1}{4} \alpha = 1

\alpha = \frac{1}{3}

Finalement:

f(t) = \frac{1}{4} - \frac{exp(-4t/3)}{4}

C’est bien gentil vous allez me dire, on connaît maintenant f(t) mais ça ne résout pas notre problème initial qui était de prendre en compte les mutations cachées… En fait si mais il reste encore un peu de calcul à faire. Pour cela on doit estimer la valeur de la variable t, c’est-à-dire la distance qui sépare nos deux séquences S_1 et S_2. Grâce à notre modèle décrit ci-dessus cette distance est bien sûr reliée aux nombres de mutations observées entre les deux séquences tout en prenant en compte le fait que certaines mutations soient arrivées sans qu’on puisse les voir.

Afin d’estimer t on va appliquer la méthode du maximum de vraisemblance. Je rappelle que la vraisemblance (notée L pour likelihood) est la probabilité d’observer les données sachant le modèle: P(data/model). Dans notre cas on veut calculer la probabilité que la séquence S_1 (de taille n) ait pu évoluer en S_2, ce qui s’écrit:

L = P(S_{1}[1]) P(S_{1}[1] -> S_{2}[1]/t) P(S_{1}[2]) P(S_{1}[2] -> S_{2}[2]/t) ... P(S_{1}[n]) P(S_{1}[n] -> S_{2}[n]/t)

avec P(S_1[i]) la probabilité d’observer le nucléotide en question au i-ème site de la séquence S_1, et P(S_1[i] \rightarrow S_2[i]/t) la probabilité d’avoir muté au i-ème site du nucléotide de S_1 vers le nucléotide de S_2 pendant le temps t.

Quand on a un produit (multiplications) on aime bien le transformer en somme (additions). Pour cela on utilise la fonction logarithme:

ln L = ln P(S_1[1]) + ... + ln P(S_1[n]) + ln P(S_1[1] \rightarrow S_2[1]/t) + ... + ln P(S_1[n] \rightarrow S_2[n]/t)

Afin de trouver le maximum de cette vraisemblance on fait comme au lycée: on dérive la fonction ln L et on cherche les valeurs auxquelles la dérivée s’annule. Les n premiers termes de la formule valent une constante donc leur dérivée est nulle. Pour les n autres on peut poser que m_1 correspond aux mutations d’un nucléotide vers un autre (p_{ij}) et m_2 correspond aux mutations d’un nucléotide vers lui-même (p_{ii}).

\frac{d (ln L)}{dt} = \frac{m_{1}}{p_{ij}(t)} p_{ij}'(t) + \frac{m_{2}}{p_{ii}(t)} p_{ii}'(t) = 0

Comme on a calculé un peu plus haut f(t) et que f(t) = p_{ij}(t) on peut remplacer dans l’équation ci-dessus. Je vous épargne les calculs mais à la fin on obtient:

\hat{t} = -\frac{3}{4} ln ( 1 - \frac{4m_{1}}{3(m_{1}+m_{2})} )

Et pour résoudre enfin notre problème on peut définir p comme étant la proportion de sites différents entre nos deux séquences. Ainsi, alors qu’on a commencé par estimer la distance entre nos deux séquences par:

p = \frac{m_{1}}{m_{1}+m_{2}}

on estime maintenant cette distance par:

\hat{t} = -\frac{3}{4} ln ( 1 - \frac{4}{3} p )

Et c’est cela qu’on appelle la distance de Jukes-Cantor. Dans le cas de la figure 1, p=0.1 alors que \hat{t}=0.107. La distance de Jukes-Cantor est bien légèrement plus grande car elle prend en compte des mutations qui ont pu arriver mais qu’on ne voit pas.

Alors bien sûr, comme toujours en modélisation, on simplifie beaucoup, mais depuis l’article de Jukes Cantor en 1969 les modèles ont été perfectionnés et cela permet de bien mieux comprendre le génome des êtres vivants: vitesse d’apparition des mutations, importance fonctionnelle de certaines séquences, relations phylogénétiques entre les espèces… Mais ce serait trop pour ce billet !

ps: une bonne revue sur ces questions est disponible ici.

 

 

 

 

 

 

 

 

Laisser un Commentaire » | Science | Tagué: , , , , , | Permalien
Publié par walrus

Uncertainty

28 février 2010

Pourquoi est-ce utile de faire un modèle statistique et pourquoi peut-on y accorder quelque confiance ?

We may at once admit that any inference from the particular to the general must be attended with some degree of uncertainty, but this is not the same as to admit that such inference cannot be rigorous, for the nature and degree of uncertainty may itself be capable of rigorous expression.

Sir Ronal A. Fisher

Laisser un Commentaire » | Citations, Science | Tagué: | Permalien
Publié par walrus

Un témoin, un policier et un statisticien voit un taxi bleu

13 septembre 2009

L’homme déambulait du côté de la gare Montparnasse quand soudain il vit un taxi bleu créer un accident et s’enfuir à toute allure. D’un coup il (l’homme) passa du statut de flâneur à celui de témoin, nettement moins réjouissant… Dans ce genre de cas, des études ont montré que le témoin a raison dans 80% des cas. Le policier quant à lui sait que 85% des taxis autour de la gare Montparnasse sont bleus. Pour savoir s’il doit croire le témoin ou non, le policier sort prestement son téléphone portable de sa poche, pianote le numéro de son pote statisticien et lui pose la question suivante: “quelle est donc la probabilité pour un taxi bleu d’être impliqué dans l’accident ?” Le statisticien raccrocha, se saisit d’un crayon et d’une feuille de papier et commença à triturer ses méninges. Voici, en gros, sa démarche.

taxi à Bangkok

Pour commencer, quelques notations. Quand on écrit P(A/B) on entend: probabilité de l’évènement A sachant que l’évènement B est survenu. Le “sachant que” est noté par le “/” dans la formule. Et puis bien sûr, si A est un évènement, \bar{A} est son opposé et on a: P(A) + P(\bar{A}) = 1.

Ensuite, la première chose à faire est d’écrire les informations connues:

P( \text{le temoin dit que le taxi est bleu} / \text{le taxi est bleu} )

= P( \text{le temoin a raison} ) = 0.8

Cette probabilité renseigne le policier sur la probabilité qu’a le témoin d’avoir raison sur le fait que le taxi soit bleu.

On sait aussi que P( \text{le taxi est bleu} ) = 0.85. Cette probabilité décrit le degré de confiance qu’a le policier en le fait que le taxi soit bleu avant même d’avoir entendu le récit du témoin. En bon anglophone, on appellera cette probabilité le prior.

Maintenant ce qui intéresse vraiment le policier:

P( \text{le taxi est bleu} / \text{le temoin dit que le taxi est bleu} )

En d’autres mots, le policier veut savoir la probabilité qu’un taxi bleu soit impliqué dans l’accident sachant les données qu’il a de la part du témoin. Accrochez-vous bien: grâce aux travaux d’un révérend anglais du XVIIe siècle, Thomas Bayes, le statisticien va pouvoir fournir une réponse au policier. Sa contribution majeure fût le so-called “théorème de Bayes“:

P( A / B ) = \frac{P( B / A ) P( A )}{P( B )}

pour deux évènements A et B, avec P(B) > 0.

Et pour P(B), on peut écrire:

P(B) = P(B/A) P(A) + P(B/\bar{A})P(\bar{A})

Ce qui est utile avec ce théorème, c’est qu’il permet de “renverser” les probabilités conditionnelles: si on connait P(A) et P(B) on peut trouver P(A/B) à partir de P(B/A) (et réciproquement).

Maintenant, on peut appliquer ce théorème à notre cas. Soit A l’évènement “le taxi est bleu” et B l’évènement “le témoin dit que le taxi est bleu”. On sait déjà que P(B/A)=0.8 et que P(A)=0.85. Il est donc facile de calculer P(B/\bar{A}):

= P( \text{le temoin dit que le taxi est bleu} / \text{le taxi n'est pas bleu} )

= P( \text{le temoin a tort} )

= 1 - P( \text{le temoin a raison} )

= 1 - 0.8

= 0.2

On peut aussi calculer:

P( \bar{A} )

= P( \text{le taxi n'est pas bleu} )

= 1 - P( \text{le taxi est bleu} )

= 1 - 0.85

= 0.15

On peut alors calculer:

P(B)

= P(B/A)P(A) + P(B/\bar{A})P(\bar{A})

= (0.8 \times 0.85 ) + ( 0.2 \times 0.15 )

= 0.71

Finalement, on répond à la question initiale:

P( \text{le taxi implique dans l'accident est bleu} / \text{le temoin dit que le taxi est bleu} )

= P( A/B )

= \frac{P(B/A)P(A)}{P(B)}

= \frac{0.8 \times 0.85}{0.71}

= 0.96

Mais attention, ce n’est pas fini, et c’est même là que ça devient vraiment intéressant, alors on ne faiblit pas et on lit jusqu’au bout !

Il y a deux probabilités-clés dans la formule de Bayes ci-dessus. La première est P( \text{le temoin dit que le taxi est bleu} / \text{le taxi est bleu} ) notée P(B/A). Puisque le témoignage est la seule donnée connue de la police et que le reste est une hypothèse (85% des taxis sont bleus), on peut voir cette probabilité comme étant P( \text{donnees} / \text{hypotheses} ). Cette probabilité est la vraisemblance (likelihood) des données en fonction des hypothèses: how likely the data are given the hypotheses.

La seconde probabilité-clé est P( \text{le taxi est bleu} ) notée P(A). Comme on l’a dit précédemment, c’est le prior, le degré de confiance qu’on a en l’hypothèse, avant d’avoir vu les données.

Grâce au théorème de Bayes, on calcule la probabilité P( \text{le taxi est bleu} / \text{donnees} ) appelée le posterior. En fonction des données que l’on observe, on ajuste notre probabilité de l’hypothèse. Il est donc aussi facile de calculer la probabilité P( \text{le taxi n'est pas bleu} / \text{le temoin dit que le taxi est bleu} ):

= P( \bar{A} / B )

= 1 - P( A / B )

= 1 - 0.96

= 0.04

Comme P( \text{le taxi est bleu} / \text{le temoin dit que le taxi est bleu} ) (qui vaut 0.96) est bien plus grande que P( \text{le taxi n'est pas bleu} / \text{le temoin dit que le taxi est bleu} ) (qui vaut 0.04), il est normal pour le policier de conclure que le taxi impliqué dans l’accident est bleu.

Jusqu’à maintenant on a simplement considéré comme données le témoignage selon lequel le taxi était bleu. Mais si le témoin dit que le taxi n’était pas bleu, que se passe-t-il ? Et bien on fait le même type de calcul et on obtient comme posterior:

P( \text{le taxi est bleu} / \text{le temoin dit que le taxi n'est pas bleu} )

= 0.59

et donc aussi:

P( \text{le taxi n'est pas bleu} / \text{le temoin dit que le taxi n'est pas bleu} )

= 0.41

Et là, stupéfait, on se rend compte que, dans ce cas-là aussi, le policier conclut que le taxi doit être bleu. Si je résume, quelque soit le témoignage de notre brave homme, le policier conclut à chaque fois que le taxi est bleu ! Le prior qu’a le policier sur la couleur des taxis a en fin de compte été plus fort que les données apportées par le témoignage.

Mais pourquoi écrire tout un billet sur cette question ? Pour plusieurs raisons, tout d’abord, si vous êtes arrivés jusqu’à ces lignes, cela signifie que vous avez compris l’usage de la formule de Bayes (yes…!). Mais surtout, ce billet me permet de parler de l’utilisation des statistiques dans l’enceinte d’un tribunal. En effet, dans des affaires délicates, les experts ès “forensic statistics” sont amenés à présenter devant juges et jurés les probabilités respectives des différents scénarios possibles. Malheureusement, juges et jurés sont rarement conscients des implications d’une telle démarche. Voici le casse-tête: qui choisit les hypothèses de départ: le juge, le procureur, la défense ? Quelles données utilise-t-on ? Comment ont-elles été obtenues ces données ? Quel modèle emploie-t-on ?

Certains sont d’ailleurs allés loin dans la réflexion, voir l’article “Don’t teach statistics to lawyers !” de Robertson et Vignaux (1998). C’est un aspect du droit qui devient de plus en plus important, notamment avec les histoires de tests ADN (voir ce blog duquel est tiré l’essentiel de mon billet). Et pour la petite histoire, comment en suis-je arriver à parler de ça ? Et bien l’un de mes amis soutient sa thèse de maths dans quelques semaines et son directeur s’est beaucoup impliqué dans une affaire judiciaire au cours de laquelle des erreurs ont été commises dans l’estimation des probabilités…

Source: merci à Franz Golhen pour la photo du taxi

4 Commentaires | Science | Tagué: , , | Permalien
Publié par walrus

Significativité d’un alignement pour BLAST

12 juillet 2009

On a vu il y a quelque temps comment estimer la significativité d’un alignement en utilisant la formule de Bayes. Mais pratiquement, ce n’est pas la méthode la plus utilisée. Il est temps d’introduire BLAST…

Derrière cet acronyme percutant se cache l’algorithme le plus utilisé en alignement de séquences biologiques. Il correspond à une heuristique d’alignement local deux-à-deux présenté dans l’article “Basic local alignment search tool” dont il tire son nom. Cet article d’Altschul, Gish, Miller, Myers et Lipman, paru en 1990, a fait date: à ce jour il a été cité dans au moins 30000 publications, un record! Il est très utilisé parce qu’il est rapide et permet de trouver les alignements significatifs entre une sequence inconnue et une (très) grande banque de séquences connues.

Cet algorithme fonctionne en plusieurs étapes. Étant donné deux séquences à comparer, il commence par rechercher des segments de longueur égale, sans gap,  identiques entre eux. On appelle ça les seeds, ce sont des k-mers avec k=11 pour des séquences nucléotidiques. Puis l’algorithme essaie d’étendre ces seeds en autorisant les gaps: il obtient alors des HSPs (pour high-scoring segment pairs). Finalement, pour chaque HSP (donc pour chaque alignement local entre les deux séquences),  l’algorithme choisit ceux dont le score est considéré comme étant significatif.

Optimisation of the Smith-Waterman algorithm

Supposons que l’on aligne une séquence de fonction inconnue, la query, de taille m, avec une séquence de fonction connue, la subject, de taille n, via BLAST. On obtient N HSPs, chacun ayant un score H. Pour estimer la significativité d’un HSP, on a besoin de savoir à quel point c’est probable d’obtenir un HSP de score H par hasard, sans que les deux séquences aient quoi que ce soit à voir entre elles. Cela revient à calculer la probabilité que le score maximal des N HSPs obtenus dépasse un seuil S donné.

Maintenant considérons que les H_i sont des variables aléatoires indépendamment distribuées selon une densité de probabilité f(H).  Si on prend leur maximum X = \max_i ( H_i ), on peut écrire:

P_N ( X \le S ) = P(H_1 \le S) ... P(H_N \le S)

P_N ( X \le S ) = \left( \int_{-\infty}^{S} p(H) \, dH \right)^N

P_N ( X \le S ) = \left( 1 - \int_{S}^{+\infty} p(H) \, dH \right)^N

Si les deux séquences ont beaucoup d’HSPs entre elles (N grand) et comme S est vraisemblablement élevé, il est dans la queue de la distribution et donc la surface sous l’intégrale est petite. On peut donc appliquer l’approximation bien connue (1-x)^k \approx \exp^{-kx} quand x est petit:

P_N ( X \le S ) \approx \exp \left( - N \int_{S}^{+\infty} p(H) \, dH \right)

Maintenant, supposons pour l’instant que p(H) tombe exponentiellement sur ces extrêmes (les queues de distribution, quand S est suffisamment grand):

\int_{S}^{+\infty} p(H) dH = \int_{S}^{+\infty} a \exp^{-\lambda H} dH = \frac{a}{\lambda} \exp^{-\lambda S}

Ainsi: P_N ( X \le S ) \approx \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Finalement, la probabilité que le score maximal X des N HSPs obtenus entre nos deux séquences dépasse un seuil S donné vaut:

P_N ( X > S ) \approx 1 - \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Cette formule nous permet donc de dire quand le score X tombe dans les 1% ou 5% de la distribution. C’est bien ce qu’on cherchait, non ?

Pour la petite histoire, X suit une distribution dite “des valeurs extrêmes”. Vous pouvez trouver ici l’article de Gumbel, en français, décrivant les premiers résultats. On se sert de cette théorie dans de nombreux domaines, notamment en finance et assurance, pour estimer la variabilité d’un système et donc les risques d’atteindre des valeurs extrêmes. Les lois des valeurs extrêmes ont l’une de leur queues de distribution plus épaisse que la loi normale, donc la probabilité d’occurrence d’une valeur extrême, c’est-à-dire tombant dans cette queue de distribution plus épaisse, est plus forte que pour une loi normale. Mais certains aiment quand même bien la loi normale, et un prix Nobel ça inspire confiance, quoique après ça, certains vont réviser leur modèles…

Loi Normale et loi de GumbelLa 1e image est tirée de “The estimation of statistical parameters for local alignment score distribution”, Althschul et al., Nucleic Acids Research, 2001. Quant à la 2e image, voici le code R pour la générer:


require( distrEx )
require( Cairo )
Cairo( "normal_vs_gumbel.png", width=700, height=500, bg="white" )
x <- seq(-4,6,0.01)
par( mar=c(5,5,3,2), font=2, font.axis=2, font.lab=2, cex=1.5, lwd=2 )
plot( x, dnorm(x), type="l", col="black", main="Loi Normale et loi de Gumbel", ylab="densité de probabilité", ylim=c(0.0,0.5), xlim=c(-4,6) )
points( x, dgumbel(x), type="l", col="red" )
abline( v=0, lty=2, lwd=1 )
legend( 0.2, 0.5, c( "densité de la loi Normale", "densité de la loi de Gumbel" ), col=c("black","red"), lwd=2, bty="n" )
dev.off()

Ajout: un tutoriel sur BLAST vient de sortir dans PLoS Biology, avec une figure qui résume très bien les différentes étapes.

2 Commentaires | Science | Tagué: , , | Permalien
Publié par walrus

Vent de révolte en Iran

20 juin 2009

Il y a quelques semaines à peine tous les citoyens de l’Union Européenne ont été appelés à voter pour élire leur Parlement. La participation fût faible… C’était pourtant l’occasion d’exprimer l’un des droits les plus importants acquis de haute lutte au cours de l’histoire, l’une des occasion de réaliser pleinement le contrat social qui nous lie tous. Mais apparemment, de telles considérations n’agitent plus les foules. dans nos contrées Cependant, la situation n’est pas la même partout et dans les pays où le mot “élections libres” est interprété par le pouvoir en place, certains votants se révoltent: c’est typiquement le cas de l’Iran en ce moment même.

Heureusement dans le monde actuel, il ne faut pas aller bien loin pour suivre ce qui arrive “en ce moment même”: Twitter. C’est un site internet de “micro-blogging” (un billet ne doit pas dépasser 140 caractères) aujourd’hui utilisé par plusieurs millions de personnes à travers le monde. L’équipe derrière ce site a ajouté une fonctionnalité de recherche à la définition évocatrice: See what’s happening — right now. Alors pour suivre ce qui se passe dans les rues de Téhéran, je suis allé ici http://search.twitter.com/ puis j’ai cliqué sur le Trending topic #iranelection et là, toutes les secondes, des billets apparaissent.

Par exemple, j’ai pu suivre en temps réel l’annonce de l’attaque suicide à la bombe sur le mausolé de Khomeini, puis le fait que cette annonce ne soit diffusée que par Press TV, organe financé par le gouvernement iranien. Idem, j’ai pu suivre l’avancée des policiers dans les rues de Téhéran, et ce aidé de Google Map.


View Larger Map

J’ai également pu voir certaines vidéos postées sur YouTube, par exemple celle-ci montrant des manifestants faisant reculer les policiers en criant “Down with Khamenei”:

ou celle-là d’une réflexion déchirante avec vue des toits de Téhéran la nuit du 19 au 20 juin:

Et parmi l’avalanche de twitt (nom donné aux messages sur Twitter), j’en ai retenu un menant à cet article sur le site du Washington Post écrit par deux doctorants de Columbia University. En bons statisticiens, lorsqu’ils ont lu les nouvelles selon lesquelles les élections présidentielles en Iran étaient truquées, ils ont décidé d’en avoir le coeur net. En premier lieu ils ont récupéré sur le site du Ministère iranien de l’Intérieur le nombre de votes reçus par les quatre candidats principaux, et ce dans plusieurs régions. Puis, si par exemple Ahmadinejad a reçu 14579 votes dans une province particulière, ils se sont intéressés aux deux derniers chiffres de ce nombre: “7″ et “9″.

Et pourquoi donc? Parce que les psychologues ont montré que les humaines n’étaient pas très bons à inventer des nombres aléatoirement: on a tendance à sélectionner certains chiffres plus souvent que d’autres. S’il n’y a pas de fraude, chaque chiffre (0,1, 2, 3…9) devrait apparaître dans des proportions uniformes à la fin des nombres de votes. Dans la même veine, on a tendance à sélectionner des nombres avec des chiffres adjacents: on préfère 23 à 17 par exemple.

En regardant ce genre de choses dans les résultats des élections iraniennes, nos deux statisticiens, Bernd Beber et Alexandra Scacco, ont montré qu’effectivement on pouvait observer des irrégularités et que celles-ci étaient très improbablement dues au hasard seul. La main des Bassidjis ?

fraude des Bassidjis versus loi statistique

Laisser un Commentaire » | Politique, Science | Tagué: | Permalien
Publié par walrus

Significativité d’un alignement en bayésien

14 mai 2009

On a vu précédemment comment aligner (globalement) deux séquences l’une avec l’autre et comment calculer le score de l’alignement résultant étant donné un scoring scheme (match, mismatch, gap). Mais bien que l’on sache obtenir un alignement optimal, cela ne nous dit pas si notre alignement est significatif, c’est-à-dire s’il est biologically meaningful…

Replaçons nous dans le contexte; nous venons d’obtenir (par des moyens expérimentaux) une nouvelle séquence. Pour en savoir plus sur l’information qu’elle code, on se dit qu’on peut essayer de l’aligner avec une séquence déjà connue. C’est ce que nous faisons et on trouve l’alignement optimal A_1 de score S_1. C’est bien, mais pas suffisant. En effet, nous n’avons (pour l’instant) aucun moyen de savoir si cet alignement est dû au hasard ou bien si ces deux séquences sont homologues, c’est-à-dire si elles ont un ancêtre commun.

Et oui, on peut très bien avoir un alignement A_2 entre notre séquence inconnue et une 2e séquence que l’on connaît déjà, de score S_2 celui-là. Or on aimerait bien avoir un moyen pour juger que deux séquences sont homologues puisque dans ce cas-là on se permettra (grossièrement) d’inférer que la fonction de notre séquence inconnue est la même que la fonction de la séquence qu’on connaît déjà (à peu de choses près).  Dans notre cas, quel alignement doit-on utiliser pour inférer la fonction ? A_1 ou A_2 ?

Et comme toujours, il existe plusieurs approches pour (tenter de) répondre à cette question, l’approche bayésienne (voir ci-dessous) et la distribution du maximum de scores (prochainement dans un autre billet).

On a vu que le score S d’un alignement était un log de ratio de vraisemblances:

S = log( \frac{P(x,y/M)}{P(x,y/R)} )

Mais en fait, au lieu de P(x,y/M), la probabilité d’obtenir les séquences x et y sachant qu’elles ont un ancêtre commun (le match modèle M), ce que l’on veut vraiment, c’est calculer la probabilité que les deux séquences aient un ancêtre commun connaissant ces deux séquences: P(M/x,y). Grâce à la formule de Bayes, on peut calculer cette probabilité:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y)}

Cependant, au préalable, il faut spécifier les prior: P(M) la probabilité a priori que les deux séquences sont reliées (par un ancêtre commun) et P(R) la probabilité a priori du modèle aléatoire:  P(R) = 1 - P(M).

Ainsi:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y/M)P(M)+P(x,y/R)P(R)}

P(M/x,y) = \frac{P(x,y/M)P(M) / P(x,y/R)P(R)}{1+P(x,y/M)P(M) / P(x,y/R)P(R)}

Finalement, on pose:

S' = S + log( \frac{P(M)}{P(R)} )

On obtient alors:

P(M/x,y) = \sigma(S') avec \sigma(x)=\frac{\exp^x}{1+\exp^x}

\sigma est la fonction logistique. C’est une sigmoide tendant vers 1 quand x tend vers \infty, et vers 0 quand x tend vers -\infty. Cette fonction est bien utile pour convertir des scores faits de sommes en probabilités.

Une fois qu’on a fait tout ça, on peut comparer la valeur finale de P(M/x,y) à 0 pour savoir si nos séquences sont reliées entre elles ou non. Sans rentrer dans les détails, il faut être bien sûr de travailler avec des probabilités (leur somme sur toutes les paires possibles de séquences doit faire 1), ce qui n’est pas toujours le cas dans un scoring scheme construit de façon ad hoc. De plus, les priors prennent beaucoup d’importance quand on cherche un alignement significatif entre une séquence inconnue et une banque de séquences connues. Dans ce cas-là, il faut prendre en compte la taille de la banque puisque la probabilité d’obtenir, par chance, un bon alignement entre deux séquences non reliées entre elles augmente avec la taille de la banque.

Laisser un Commentaire » | Science | Tagué: , , , | Permalien
Publié par walrus

« Entrées Précédentes
Suivre

Get every new post delivered to your Inbox.