Normalisation de mes donnees et aire sous la courbe

[Myriam Croze]

Bonjour,

J'ai quelques notion sur R mais helas insuffisante pour realiser un graphe a partir de mes donnees.
Voila, j'ai un tableau de ce trype:
Start End Midpoint NumSites Tajima_D FuLi_Dstar
1 500 251 500 -1.5347037 -1.6652469
251 750 501 500 -1.3106336 -1.409583
501 1000 751 488 -1.3258685 -1.403907
751 1250 1001 464 -1.3387919 -1.3992572
1001 1500 1251 476 -0.787162 -1.1232101
1251 1750 1501 474 -0.7785037 -1.1259229
1501 2000 1751 474 -1.#IND000 -1.#IND000
1751 2250 2001 500 -1.#IND000 -1.#IND000

Je veux faire un graphe de type normal pour les valeurs de Tajima_D.
Normalement si je prend la frequence de chaque valeur ou leur densite, je devrait obtenir une courbe de type normal.

Pouvez vous me dire comment faire ca sur R?

De plus, j'aimerai, une fois que j'ai ma courbe, delimiter l'aire sous la courbe.
Concretement, je voudrait savoir a partir de qu'elle seuil je n'ai que 5 % de valeurs maximales.

Merci par avance pour l'aide que vous pouvez m'apporter.

[Pierre COLIN]

Bonjour,

pour comparer des densités (f1 la densité d'une loi normale, f2 la densité empirique)

Code : Tout sélectionner

tab=read.table("data.txt",header=T)

m=mean(tab$Tajima_D)
sd=sd(tab$Tajima_D)
x1=seq(qnorm(0.005,mean=m,sd=sd),qnorm(0.995,mean=m,sd=sd),length.out=1000)
f1=dnorm(x1,mean=m,sd=sd)

dens=density(tab$Tajima_D,n=1000,from=qnorm(0.005,mean=m,sd=sd),to=qnorm(0.995,mean=m,sd=sd))

x2=dens$x
f2=dens$y

plot(x1,f1,ylim=c(0,max(f1,f2)),type="l",col="blue",xlab="Tajima_D",ylab="Densité")
lines(x2,f2,col="red")

pour ce qui est de limiter la courbe, il faut utiliser la fonction quantile

Code : Tout sélectionner

quantile(tab$Tajima_D,probs=0.95)

[Myriam Croze]

Bonjour et merci pour la reponse.

Mais je n'arrive pas a faire marcher ce code. Quand je le lance, voila ce que j'obtiens:

Code : Tout sélectionner

> tab=read.table("D:/PhD/Myriam_TajimasD.csv",header=T)
 
> m=mean(tab$Tajima_D)
Warning message:
In mean.default(tab$Tajima_D) :
  argument is not numeric or logical: returning NA

> sd=sd(tab$Tajima_D)
Error in var(as.vector(x), na.rm = na.rm) : 'x' is NULL

>x1=seq(qnorm(0.005,mean=m,sd=sd),qnorm(0.995,mean=m,sd=sd),length.out=1000)
Error in qnorm(p, mean, sd, lower.tail, log.p) :
  Non-numeric argument to mathematical function

> f1=dnorm(x1,mean=m,sd=sd)
Error in dnorm(x1, mean = m, sd = sd) : object 'x1' not found
>dens=density(tab$Tajima_D,n=1000,from=qnorm(0.005,mean=m,sd=sd),to=qnorm(0.995,mean=m,sd=sd))
Error in density.default(tab$Tajima_D, n = 1000, from = qnorm(0.005, mean = m,  :
 argument 'x' must be numeric
 
> x2=dens$x
Error: object 'dens' not found
> f2=dens$y
Error: object 'dens' not found

>plot(x1,f1,ylim=c(0,max(f1,f2)),type="l",col="blue",xlab="Tajima_D",ylab="Densité")
Error in plot(x1, f1, ylim = c(0, max(f1, f2)), type = "l", col = "blue",  :
  object 'x1' not found

> lines(x2,f2,col="red")
Error in lines(x2, f2, col = "red") : object 'x2' not found

Peut etre que cela vient du fait que j'ai parfois des valeurs = NA que je ne veux pas prendre em compte justement.

Sinon, j'ai reussi a faire un histogramme "Barplot" :

Code : Tout sélectionner

TajimasD_table<-read.table("D:/PhD/Myriam_TajimasD.csv",
                                          sep=",",
                                          header=TRUE)

TajimasD<-TajimasD_table$Tajima_D

hist(TajimasD, breaks=100)
rug(TajimasD)

Serait-il possible a partir de cet histogramme de tracer une courbe ou si ce n'est pas possible, juste obtenir les 5% de valeurs maximum toujours sans prendre en compte les valeurs = NA.

Merci.

[Pierre COLIN]

l'erreur vient surement des valeurs "-1.#IND000" présente dans tes données.

Pour gérer les valeurs manquantes tu peux utiliser la fonction "is.na"

Code : Tout sélectionner

Tajima_D=tab$Tajima_D[!is.na(tab$Tajima_D)]

ou bien l'argument "na.rm=TRUE" qui est présent dans pas mal de fonction

[Myriam Croze]

Merci mais helas, je n'y arrive toujours pas.
Maintenant, j'obtiens ca:

Code : Tout sélectionner

> tab<-read.table("D:/PhD/Myriam_TajimasD.csv",sep=",",header=T)
 
> Tajima_D<-tab$Tajima_D[!is.na(tab$Tajima_D)]
 
> m=mean(tab$Tajima_D)
> sd=sd(tab$Tajima_D)
 >x1=seq(qnorm(0.005,mean=m,sd=sd),qnorm(0.995,mean=m,sd=sd),length.out=1000)
Error in if (from == to) rep.int(from, length.out) else as.vector(c(from,  :
  missing value where TRUE/FALSE needed

> f1=dnorm(x1,mean=m,sd=sd)
Error in dnorm(x1, mean = m, sd = sd) : object 'x1' not found

 >dens=density(tab$Tajima_D,n=1000,from=qnorm(0.005,mean=m,sd=sd),to=qnorm(0.995,mean=m,sd=sd))
Error in density.default(tab$Tajima_D, n = 1000, from = qnorm(0.005, mean = m,  :
  'x' contains missing values
 
> x2=dens$x
Error: object 'dens' not found
> f2=dens$y
Error: object 'dens' not found
 >plot(x1,f1,ylim=c(0,max(f1,f2)),type="l",col="blue",xlab="Tajima_D",ylab="Densité")
Error in plot(x1, f1, ylim = c(0, max(f1, f2)), type = "l", col = "blue",  :
  object 'x1' not found
> lines(x2,f2,col="red")


Je ne vois pas ce qui pourrait gener.

Sinon, j'ai essaye d'appliquer le code des quantiles sur mon histogramme,
mais la encore ca ne marche pas et j'obtiens ceci:

Code : Tout sélectionner

> TajimasD_table<-read.table("D:/PhD/Myriam_TajimasD.csv",
+                                           sep=",",
+                                           header=TRUE, na.rm=TRUE)
Error in read.table("D:/PhD/Myriam_TajimasD.csv", sep = ",", header = TRUE,  :
  unused argument(s) (na.rm = TRUE)
>
> TajimasD<-TajimasD_table$Tajima_D
>
> hist(TajimasD, breaks=100)
> rug(TajimasD)
>
> quantile(TajimasD_table$Tajima_D,probs=0.95)
Error in quantile.default(TajimasD_table$Tajima_D, probs = 0.95) :
  missing values and NaN's not allowed if 'na.rm' is FALSE


La encore je ne vois pas ce qui gene.

Par ailleurs, j'ai supprimer toute les valeurs=NA dans mon tableau et meme en faisant ca, le code ne marche toujours pas.

Avez-vous une idee du probleme?
Merci.

[Serge Rapenne]

Bonjour,

il n'y a pas de paramètre na.rm dans read.table mais un paramètre na.strings qui permet d'indiquer quelle chaine de caractère doit être considérée comme NA donc ton read.table doit s’écrire :

Code : Tout sélectionner

tab<-read.table("D:/PhD/Myriam_TajimasD.csv",sep=",",header=T,na.strings="-1.#IND000")

par contre la fonction quantile a elle un paramètre na.rm et donc il doit s'écrire

Code : Tout sélectionner

quantile(TajimasD_table$Tajima_D,probs=0.95,na.rm=TRUE)

pour obtenir l'aide et donc la syntaxe et les paramètres d'une fonction utilise ?

Code : Tout sélectionner

?read.table #aide de read.table
?quantile  #aide de quantile

Serge

[Myriam Croze]

Merci pour ces informations.
Je vais regarder l'aide en esperant finir par y arriver.

[Myriam Croze]

Bonjour,
j'ai finalement reussi a faire mon histogramme et a avoir la valeur des 5%.
De meme, j'arrive a tracer la courbe des densites.

Voila mon code:

Code : Tout sélectionner

# read in the data
TajimasD_table<-read.table("D:/PhD/Myriam_TajimasD.csv", sep=",", header=TRUE)
# store Tajima"s D data in a vector
TajimasD<-TajimasD_table$Tajima_D
# create a layout for 2 plots, with 75% of vertical space going to the upper, 25% to the lower
layout(matrix(seq(2)),heights=c(0.75,0.25))
# reduce the size of the lower margin
par(mar=c(0,4.1,4.1,2.1))
# plot the histogram with 100 bins, topographic colours, no axes or axis names, and custom title
hist(TajimasD, breaks=100, col=topo.colors(53), axes=FALSE, xlab="", ylab="", main="Tajima's D, Run Mode 12")
# add the y-axis
axis(2)
# add the label to the y-axis, left side (2), with 2.5 margin
mtext("Frequency", side=2, line=2.5, font=2)
# add the rug plot of a desired size, colour it, bottom side (1), position it under histogram
rug(TajimasD, ticksize=0.03, col="blue", side=1, pos=-50)
# adjust margins and axis location
par(mar=c(5.1,4.1,0,2.1), mgp=c(3,0.5,0.0))
# add a boxplot underneath, a horizontal one, without its own axes
boxplot(TajimasD, horizontal=TRUE, axes=FALSE)
# add the x-axis
axis(1)
# add the x-axis label
mtext("Tajima's D", side=1,line=2.5, font=2)

dens<-density(TajimasD, na.rm=TRUE)
plot (dens, na.rm=TRUE)

quantile(TajimasD_table$Tajima_D,probs=0.95,na.rm=TRUE)

Le probleme, c'est que je voudrais tracer la courbe de densite par dessus l'histogramme (pour voir les 2 ensembles).
malgre de nombreuses tentatives, je n'arrive qu'a avoir mes deux graphes separes.

De plus, existe t'il une fonction qui me permettrai de colorier les derniers 5% de la courbe ou de l'histogramme?

Merci.

[Dominique Soudant]

Code : Tout sélectionner

x <- rnorm(100)
hist(x)
par(new=TRUE)
density(x)

[Vincent Guillemot]

Code : Tout sélectionner

x <- rnorm(100)
hist(x)
par(new=TRUE)
density(x)

Je ne comprends pas cet exemple, pour afficher la densité et l'histogramme en même temps, est-ce qu'il ne vaudrait mieux pas utiliser le code suivant ?

Code : Tout sélectionner

x <- rnorm(100)
hist(x,freq=F)
lines(density(x))

V.

[Myriam Croze]

Merci pour l'aide.
J'ai reussi a tracer ma courbe avec le code de Vincent.