silico.biotoul.fr
 

M1 Traitement de Donnees Biologiques - TP 3 R

From silico.biotoul.fr

Revision as of 14:50, 19 September 2016 by Bonhomme (Talk | contribs)
Jump to: navigation, search

Contents

Tests statistiques

La plupart des tests statistiques courants (et bien d’autres) sont programmés dans R. Dans ce TP, nous allons aborder quelques types de tests statistiques très courants:

  • Les tests de Chi2 (test de conformité, test d'indépendance), basés sur les effectifs de différentes classes
  • Les tests d'adéquation (ici nous verrons l'adéquation à la loi Normale)
  • Les tests d'homogénéité (comparaison de moyennes, comparaison de variances)

ATTENTION: pour chaque test statistique, il est crucial de bien connaître l'hypothèse nulle du test (H0) pour pouvoir correctement interpréter !

Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP3) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : M1.TDB.TP_tests_R.Rmd (click droit de la souris -- enregistrer la cible sous...). Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.

Tests de Chi2

Tests de Chi2 de conformité

On cherche à tester la conformité d'une (ou plusieurs) valeur(s) par rapport à une (des) valeur(s) standard(s) (ou théoriques). Le test de Chi2 de conformité s'effectue sur les effectifs des différentes classes d'une variable qualitative (ou quantitative).

Notre exemple porte sur les proportion mendéliennes dans une descendance issue d'autofécondation d'une F1 hétérozygote. Image:mendel.jpeg

Pour un gène avec 2 allèles, les proportions génotypiques attendues (relatives) dans la descendance sont:

freq_exp =c(?,?,?) # mettre sous forme de fréquences relatives

Suite à un croisement effectué en laboratoire, les proportions génotypiques observées (absolues) parmis 100 individus descendants sont:

freq_obs =c(22,53,25)

Question: les proportions observées sont-elles conformes aux proportions mendéliennes? Pour cela il faut répondre aux questions suivantes:

  • Quelle est l'hypothèse nulle (H0)du test ? L'hypothèse alternative (H1) ?
  • Qu'indique la p-valeur ?
  • La p-valeur est-elle inférieure ou supérieure au seuil alpha = 5% ?

Conclusion:

Ajoutez ces parties à votre compte rendu.




Lecture et exploration du tableau de données

Lecture:

genomes=read.table("bacterial_genomes.txt", sep="\t", header=TRUE)

Accéder directement aux variables simplement en donnant leurs noms:

attach(genomes)
names(genomes)
Genome_size
ORF_number

Faire un graphique:

plot(Genome_size,ORF_number,pch=16)

Quantifier la relation entre ces 2 variables:

  • covariance
cov(Genome_size,ORF_number)
  • coefficient de corrélation r de Pearson
cor(Genome_size,ORF_number)

Régression linéaire

Calculer une fonction linéaire qui relie les 2 variables, avec la commande lm():

lm(ORF_number ~ Genome_size) # remarque: c'est une régression de "y" sur "x", d'où lm(y~x)

Gardons en mémoire le résultat de la régression:

reglin=lm(ORF_number ~ Genome_size)

On peut vérifier la significativité des coefficients de la droite avec:

summary(reglin)

Quelle est l'équation de la droite de régression ?

Calculez le coefficient de détermination R2 (% de variance expliquée par le modèle linéaire ==> bien si > 70%):

cor(Genome_size,ORF_number)^2 (stocké aussi dans summary(reglin))

Représenter le nuage de points avec la droite de régression:

plot(Genome_size,ORF_number,pch=16)			
abline(reglin,col="red",lwd=2)

Ajoutez ces parties à votre compte rendu.


Calculs de probabilité

Calculs basés sur des données

L'exemple porte sur la taille des arbres dans une forêt de sequoia. Téléchargez et sauvegardez le fichier sequoia.txt (click droit de la souris -- enregistrer la cible sous...).

Lecture du tableau de données:

sequoia=read.table("sequoia.txt", sep="\t", header=TRUE)
attach(sequoia);names(sequoia)

Représentation graphique des données:

hist(taille_arbre)

Quelle est la probabilité qu'un arbre mesure 80m ?

length(taille_arbre[taille_arbre==80])/length(taille_arbre)

Quelle est la probabilité qu'un arbre mesure plus de 100 ?

length(taille_arbre[taille_arbre>100])/length(taille_arbre)

Calculs basés sur des lois (ou distributions) de probabilité

Différentes lois sont programmées (binomiale,Poisson,normale,chi2...) Il y a plusieurs fonctions pour chaque loi. Par exemple, pour la loi normale:

  • dnorm() : fonction densité (density)
  • pnorm() : fonction de répartition (probability)
  • qnorm() : fonction quantile (quantile)
  • rnorm() : générateur aléatoire (random)

Calculs basés sur un loi normale centrée réduite N(0,1) (0 et 1 = moyenne et écart-type de la variable X)

Image:loi_normale.jpeg

Prob[X=-1]:

dnorm(-1,0,1)

Prob[X<-1]:

pnorm(-1,0,1)

Prob[X>-1]:

1-pnorm(-1,0,1) # ou pnorm(-1,0,1, lower.tail=F)

Echantillonage + histogramme de 100 valeurs dans une N(0,1):

hist(rnorm(100,0,1))

Calculs basés sur un loi Binomiale B(n,p)

Exemple: pour un lot de n=100 graines d'Arabidopsis ayant chacune une probabilité de germination p=0.8 (B(100,0.8)), je veux:

La probabilité que k=80 graines germent (k=nombre de succès)

dbinom(80,100,prob=0.8)

La probabilité que au maximum 80 graines germent

pbinom(80,100,prob=0.8)

La probabilité que plus de 80 graines germente

pbinom(80,100,prob=0.8, lower.tail=F)

Ajoutez ces parties à votre compte rendu.


Estimer un intervalle de confiance (IC) d'un paramètre populationnel (µ, p,...) à partir d'un échantillon

IC de la moyenne µ de la population, à partir de la moyenne m calculée sur un échantillon

Calculer l'IC du poids moyen d'un récolte de tomate cerises, à partir d'un échantillon. Téléchargez et sauvegardez le fichier le fichier tomates_cerises.txt (click droit de la souris -- enregistrer la cible sous...).

Lecture du tableau de données:

tomates=read.table("tomates_cerises.txt", sep="\t", header=TRUE)
attach(tomates);names(tomates)
poids_tomate

Script pour calculer l'IC:

m = mean(poids_tomate)  # moyenne calculée sur l'échantillon
s = sd(poids_tomate)    # écart-type calculé sur l'échantillon
n = length(poids_tomate)# taille de l'échantillon

# si n<100, on utilise la distribution de Student pour calculer l'erreur.

error <- qt(0.975,df=n-1)*s/sqrt(n) # qt(0.975) = quantile 97,5% de la distribution de student; s/sqrt(n) =erreur standard de la moyenne
left  <- m-error
right <- m+error

L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est:

left;right

# si n>100 on utilise la distribution normale pour calculer l'erreur.

error <- qnorm(0.975)*s/sqrt(n) # qnorm(0.975) = quantile 97,5% de la distribution normale
left  <- m-error
right <- m+error

L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est alors

left;right

REMARQUE IMPORTANTE: si l'on veut un IC à 99% ==> remplacer 0.975 par 0.995.


IC de la fréquence f d'un caractère dans la population, à partir de la proportion p calculée sur un échantillon

Calculer l'IC de la fréquence de pommes rouges dans une récolte comportant des pommes rouges ET vertes, à partir d'un échantillon de n=125 pommes contenant p=40% de pommes rouges On utilise une approximation normale de la distribution d'une proportion

p <- 0.4	# proportion de pommes rouges dans l'échantillon
n <- 125	# taille de l'échantillon
error <- qnorm(0.975)*sqrt((p*(1-p)/n))
left  <- p-error
right <- p+error

L'intervalle ayant 95% de chance de contenir la vraie proportion de la population est:

left;right

Faites varier la taille de l'échantillon ainsi que la confiance ...


Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail (bonhomme@lrsv.ups-tlse.fr ou barriot@biotoul.fr). Le compte rendu est à envoyer avant de commencer le TP4. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP3 TDB de -et votre Nom et Prénom-".

Liens