Le blog du Doc’

Mes tribulations mathématiques m’ont récemment poussé dans le monde assez fabuleux des tests par permutation, et je viens de consacrer une grosse demi heure à écrire pour R le test T de Student par permutations, pour m’affranchir de l’homogénéité des variances et de la normalité des données.

Vu que c’est du code qui peut servir, je le donne ici… Et éventuellement, ça permettra aux R-Guru qui traînent (j’espère) de vérifier que je ne me sois pas planté…

TTperm <- function(sample1,sample2,fac1='a',fac2='b',nperm=999){
require(car)
stu<-t.test(sample1,sample2)
stu.p<-stu$p.value
stu.t<-stu$statistic
stu.df<-stu$parameter
count <- 1
result.stu <- matrix(rep(0,3),1)
colnames(result.stu)=c("T statistic (Df)","Parametric p-value", "Permutations p-value")
for(i in 1:nperm){
NAMES<-c(rep(fac1,length(sample1)),rep(fac2,length(sample2)))
SAMPLE<-c(sample1,sample2)
PERMNAMES<-sample(NAMES)
MATR<-c(SAMPLE,PERMNAMES)
MATR<-matrix(MATR,nrow=2,byrow=T)
RS1<-RS2<-''
rownames(MATR)<-c('Value','Sample')
for(rnk in 1:length(PERMNAMES)){
if(MATR['Sample',rnk]==fac1){
RS1<-c(RS1,as.real(MATR['Value',rnk]))
}
if(MATR['Sample',rnk]==fac2){
RS2<-c(RS2,as.real(MATR['Value',rnk]))
}
}
RS1<-RS1[2:length(RS1)]
RS2<-RS2[2:length(RS2)]
stu.perm<-t.test(as.real(RS1),as.real(RS2))
stu.perm.p<-stu.perm$p.value
stu.perm.t<-stu.perm$statistic
if(abs(stu.perm.t) >= abs(stu.t)) {count <- count+1}
}
result.stu[1,1]=paste(round(stu.t,4),'(',round(stu.df,4),')',collapse='')
result.stu[1,2]=round(stu.p,6)
result.stu[1,3]=count/(nperm+1)
return(result.stu)
}

Ca commençait a me manquer, les maths. Pas vous? Alors je me suis lancé dans un petit jeu avec des lettres grecques, mon tableur préféré, et un peu d’imagination. L’objectif de ce billet, c’est surtout de me faire taper sur les doigts, puisque j’imagine que ce que j’ai produit est… faux. Ou au moins assez incomplet.

Il me fallait avant tout un problème sur lequel déchaîner la toute puissance des mathématiques, et j’en suis arrivé a repenser à la compétition apparente (une lecture rapide du billet en question peut aider a comprendre la suite). Je me suis donc armé de ma feuille et de mon stylo préféré, pour faire un petit dessin. Deux hôtes, des flèches de l’un à l’autre, et des flèches qui rentrent dans chaque hôte pour symboliser les parasites. Des flèches qui sortent de l’hôte pour la mortalité, et on avait une idée de ce qui arrivait au système.

Bien, restait à mettre le tout en équations. Par ou commencer? Les choses simples. On a des populations hôtes (H₁ et H₂), avec des parasites qui rentrent, et de la natalité. Les parasites qui rentrent sont appelés β, la natalité se voit attribuer la lettre λ. Qui dit natalité dit mortalité, ici représentée par la lettre μ.

La partie à peu près triviale étant finie, on passe à la suite. Nos populations s’échangent des parasites, c’est le principe même de la compétition apparente. On va dire que le taux d’échange est γ. H₁ transmet γ₁H₁ parasites à H₂, et réciproquement. On peut faire la première critique ici : je détermine l’échange en fonction de la taille de H₁, et pas en fonction du nombre de parasites qui infestent H₁ (différence subtile). Soit on considère que c’est à peu près équivalent, soit on complexifie le modèle plus tard, et on incorpore β à la définition de γ. J’ai même fait un essai avec une réduction du nombre de parasites (δ) au sein de l’hôte, mais ça devenait lourdingue.

Etape suivante? On va faire diminuer la taille de nos populations, à cause de la compétition (d’une part) et du parasitisme (d’autre part). Les réductions ont reçu la lettre α. On a donc α_2/1 qui est la réduction causée par la population H₂ sur la population H₁ (et que j’ai calculé comme α_2/1H₂ pour prendre en compte l’importance des tailles des populations dans la compétition), et inversement. Et puis α₁, qui correspond à la mortalité induite par le parasite sur la population 1, et ainsi de suite.

Donc, si on fait un bilan. A chaque génération, notre population 1 “gagne” (et pour la 2, vous inversez 1 et 2) : H₁+λ₁ (ce qu’on avait déjà, plus la natalité). En revanche, elle perd μ₁H₁ + α_2/1H₂ (la mortalité naturelle, plus celle liée à la compétition), et perd encore à cause du parasitisme α₁((β₁-γ₁H₁)+γ₂H₂)H₁. Soit le taux de mortalité (α₁) fois la taille de la population (H₁) fois la “quantité” (une estimation discutable) de parasites, approchée par l’apport extérieur β₁, moins ce qu’on donne à la population 2 (γ₁H₁), plus ce qu’elle nous envoie (γ₂H₂).

Le résultat est présenté sur la figure de gauche. Vous aurez remarqué que la mortalité induite par les parasites pour H1 et H2 est inversée. Problème de copier/coller des formules.

Il ne reste plus qu’à rentrer tout ça dans notre tableur préféré, et fixer des paramètres. Combien? Vous aurez compris, il y en a 12. Et comme on connaît un peu la compétition apparente, on a des indications.

Ce que j’ai choisi, c’est évidemment α_2/1<α_1/2. Et le reste? J’ai fouillé un peu pour trouver les meilleurs valeurs possibles (disons le tout de suite, celles qui donnent les plus belles courbes).

Ce qui donne :

Comme vous en mourrez d’envie, je fais tourner le modèle. Une petite indication avant. La population 1 (celle qui va utiliser le parasite a son avantage, comme vous l’aurez compris en observant les paramètres) est en rouge, et la 2 en bleu. En vert, j’ai représenté le ratio H1/H2, qui lorsqu’il est supérieur à 1 indique que H1 prédomine (il est en vert, et se lit sur l’axe de droite). Vous pouvez voir le résultat sur la figure de gauche.

Sans surprise, on voit la population 1 prendre l’avantage, pendant que 2 se fait décimer par les parasites. Avec d’autres paramètres, on arrive a stabiliser le modèle en quelques générations (ce qui m’embête un peu, parce que j’attend des variations cycliques, et que comme j’ai construit mes équations “un peu” à la va vite, elles n’y sont pas).

Oui mais voilà, ce cas un peu simple, il m’ennuie rapidement. Alors j’ai fait en sorte d’introduire le parasite après quelques générations. Regardez les paramètres. La population 1 a une mortalité plus forte, une natalité plus faible. Et comme on s’y attend, la population 2 prend l’avantage.

Le parasite est introduit à la génération 15, et la on observe ce qui était prévu : la population 1 reprend le dessus, profitant du désavantage imposé à la population 2 par le parasite (et le fort taux de transmission de 1 vers 2 y participe beaucoup). L’introduction du parasite a “renversé la tendance” dans la compétition (on peut, avec des valeurs de paramètres différentes, accentuer ou atténuer, et beaucoup d’autres choses encore).

Et on s’arrête là? Et non. Parce que penser que les parasites qui rentrent dans le système le font a une vitesse constante, ca fait des courbes lisses, et ce n’est pas tout à fait vrai (je n’en suis pas non plus à une approximation près, mais passons).

Donc, j’ai introduit une borne supérieure et inférieure pour le nombre de parasites qui peuvent rentrer dans la population 1 et 2, et pour chaque génération, j’ai demandé un tirage aléatoire. En l’occurence, β₁ entre 1 et 4, et β₂ entre 1 et 2. Ah ah, vous dites vous, il va trop loin, si il a trop de parasites sur 1, il ne prendra pas l’avantage. Détrompez vous, en montant jusqu’à 7 ou 8 la borne de β₁, j’avais encore le même type de résultats. Celui que vous pouvez voir dans l figure de gauche.

Et vous remarquerez que quand on enlève le parasite (aux environs de la génération 40), 1 perd son avantage au profit de 2.

Donc, les parasites peuvent tout à fait modifier les relations entre deux espèces/populations, et conférer un avantage important à l’une d’entre elles.

Quelques remarques quand même. Dans ce modèle, j’ai soustrait les parasites transmis de 1 vers 2 et réciproquement du total entrant. Ca implique un cycle direct, et même de la microprédation. D’autre part, je n’ai pas du tout pris en compte la population de parasites elle-même (on s’attend a des réductions dans l’hôte, ce qui obligerait à reformuler β sous la forme β-δβ, avec δ un coefficient de réduction).

Bref, il y a beaucoup a modifier sur ce modèle, c’est l’intérêt des commentaires (par ailleurs, si l’un d’entre vous a vu passer une vraie modélisation de ce phénomène, je suis prenneur).

Dans un précédent billet, je faisais part de mes remarques quant à l’utilisation de la morphologie pour identifier les parasites, dans la mesure ou plusieurs paramètres pouvaient venir perturber le phénomène. Comme j’avais un peu de temps hier après-midi au labo, et bien envie de jouer avec la version 7 de JMP, je me suis lancé dans une rapide analyse de mes données, pour voir si c’était plutôt encourageant ou pas…

Rapide analyse est à prendre au sens “Analyse en composantes principales”, cela va de soi. J’explique. Je prend des mesures de différentes parties des organes d’attachement de mes parasites. En suivant un bête raisonnement, on peut se dire que chaque valeur correspond à une coordonnée. Donc, il devient possible de représenter tout ça graphiquement. Problème, dans mon cas, on se retrouve avec un espace à 11 dimensions. Pas évident à visualiser, donc.

Le principe de l’ACP est qu’elle permet de représenter toutes les données sous forme de points, dans un espace à deux dimensions. Autrement dit, un basique repère avec un axe x et un axe y. Je passe sur les détails, il semblerait que savoir ça suffit à comprendre l’intêret de l’analyse.

Bref, cette analyse va avoir tendance à “regrouper” les points (donc, dans mon cas, les parasites) qui sont “proches”, sur l’ensemble de leurs variables. En gros, pour une tendance égale, on va avoir des “nuages de points” plus ou moins denses.

Me voila donc tout confiant, en train de lancer mon analyse. Et la, j’ai obtenu quelque chose de particulièrement intéressant: un nuage de points qui semblait s’être “coupé en deux”, avec une partie des points de chaque côté de l’axe. Allons donc. Si vous cliquez sur la miniature à côté de ce paragraphe, et que vous regardez comment se distribuent les points vert avec un marqueur +, vous verrez ce que je veux dire.

Nous avons donc la situation suivante: au sein d’une seule espèce de parasites (très clairement identifiable, avec les critères que l’on a, c’est du 0% d’erreur), on peut distinguer deux groupes sur la base de critères morphomètriques (uniquement ceux de leur organe d’attachement). C’est tout? Bien sûr que non. En jouant un peu avec les légendes, je me suis rendu compte que tous les parasites “de gauche” (ceux qui sont à part) ont été prélevés sur le même hôte. Et que ceux “de droite” (soit la majorité) sont issus de l’ensemble des autres hôtes. Quelle différence entre ces hôtes? L’hôte des parasites “de gauche” est d’une espèce (et d’un genre) différente de celle de ses petits camarades.

Bien sûr, il y a un certain nombre de bémols à apporter à ce qui semble être un résultat plutôt intéressant. D’une part, l’individu de l’espèce qui semble avoir des parasites avec une morphométrie différente était le seul de son espèce à avoir été examiné. Mais on va y remédier (pas avant janvier, je vous raconte pas le suspense). D’autre part, il est possible que les parasites que j’ai examiné soient à un stade plus précoce de leur développement. J’y crois moyennement (dans la mesure ou après avoir vu ce résultat, j’ai récupéré des parasites de l’individu en question conservés dans l’alcool, et qu’ils s’insèrent parfaitement au milieu des autres, ce qui indiquerait éventuellement une synchronicité de la population), mais ce n’est pas impossible.

En tout cas, ce résultat est intéressant. Il va maintenant falloir (en plus de la vérification sur d’autres hôtes) regarder la variabilité génétique, et qui sait, en tirer des conclusions…

Je viens de recevoir un appel sur mon portable, d’un numéro en 08. En fait d’appel, il s’agissait plutôt d’un “bip”. Etant, pour une fois, d’humeur sociale (j’ai fini mon powerpoint de soutenance, et je le fais en 20 minutes tout pile), j’ai décidé de rappeler pour voir de quoi il retournait.
Lire la suite »

Je n’imaginais pas me retrouver un jour face à un « cas de conscience » provoqué par des statistiques. Pourtant, cet après-midi, il a fallu faire un choix entre les deux choses que nous chérissons le plus au monde, notre infaillible feeling et la non moins intouchable p-value…
Lire la suite »