jueves, 26 de mayo de 2011

Calculo de nùmero de Clases por la regla de Sturges

Muchas veces es necesario encontrar el numero de clases que es estadisticamente optimo para realizar un histograma. Herbert Sturger en 1926, desarrollo un método  basado en el N de las muestras para encontrar este numero de clases y su amplitud de rango.
Este metodo es útil también para definir estadìos de desarrollo, basados por ejemplo en medidas de distancia.
Obviamente esta función también ha sido escrita para R y es realmente muy sencilla, así que manos a la obra:


>estadios<-c(18.30,18.70,19.60,20.20,21.00,22.70,22.75,22.50,22.89,22.90,23.00,23.20,23.28,23.30,23.40,24.07,24.40,24.78,24.90,25.41,25.60,26.10,26.20,26.60,25.60,26.13,26.30,26.31,26.56,26.60,27.04,26.80,27.50,28.00,28.32,29.50,28.90,29.00,29.10,29.30,30.20,30.35,30.50,30.50,30.50,29.60,30.50,30.60,30.90,31.20,32.00,31.40,31.73,32.60,32.04,32.70,34.00,34.10,42.40,42.00,45.10,45.50,47.90,47.00,47.80,48.60,53.70,55.80,52.90,57.30,56.20,58.40,60.40,59.90,62.30,63.00,63.90,65.30,71.20,68.70,70.40,71.30,73.30,78.60,80.20,92.10,72.02,89.96,67.59,73.23,71.96,73.13,74.20,78.92,85.36,88.55,88.71,88.25,97.38,98.09,99.01,97.94,46.09,55.26,56.91,64.63,62.84,67.01,65.65,68.63,71.19,71.28,71.62,72.45,75.37,80.33,86.07,89.90,99.57,89.52,84.22,78.47,74.48,75.72,71.26,69.75,67.83,67.00,63.68,65.50,58.66,57.70,53.68,38.24,33.78,31.12,32.69,30.94,29.83,30.36,29.80,28.54,28.82,28.04,26.45,26.33,25.59,26.64,24.48,22.64,24.13,22.99,21.67,100.63,98.32,88.10,87.39,81.44,83.92,74.31,69.22,60.43,110.10,99.65,62.39,92.00,88.64,83.74,79.71,77.74,75.69,71.64,73.57,72.70,70.79,67.99,65.28,61.86,61.30,60.69,56.09,51.70,49.66,49.25,45.36,40.22,29.99,31.68,27.96,27.99,27.82,26.41,27.98,25.36,27.41,26.56,26.92,24.14,25.28,23.20,23.80,20.91,21.65,14.23,21.03,16.62,16.63)


#Escribimos un vector con la lista de medidas que tenemos, en este caso se utilizan las medidas de longitudes estándar de individuos de una especie de Pez de la familia Loricariidae, colectados en un Rio.


> nclass.Sturges(estadios)


#Obtenemos el Nùmero de clases en los que se debe dividir la muestra, para que sea consistente, en este caso serà 9.
[1] 9


> summary (estadios)


#Obtenemos los estadísticos descriptivos básicos que serán utilizados para hallar la amplitud de rango de cada estadìo o clase.
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  14.23   26.72   45.36   49.56   71.23  110.10



> RangoTotal<-(110.10-14.23)


#Obtenemos el rango total que tiene muestra muestra.


> RangoTotal
[1] 95.87



> amplitudtotal<-(RangoTotal+1)


#Sumamos 1 al rango total, para obtener la amplitud de rango total, necesaria para calcular la amplitud de rango para nuestras clases.


> amplitudtotal
[1] 96.87



> amplitudderango<-(amplitudtotal/9)


#Hallamos la amplitud de rango de nuestras clases o estadìos.


> amplitudderango
[1] 10.76333



Y listo!!!! de este modo, tendremos en este ejemplo 9 clases o estadios de amplitud de rango de 10.76333


así:


Estadio 1: 14.23 a 24.99
Estadio 2: 25.00 a 35.75
Estadio 3: 35.76 a 46.51
Estadio 4: 46.52 a 57.27
Estadio 5: 57.28 a 68.03
Estadio 6: 68.04 a 78.79
Estadio 7: 78.80 a 89.55
Estadio 8: 89.56 a 100.31
Estadio 9: 100.31 a 111.07



Que les sirva y que les rinda!!!!







jueves, 5 de mayo de 2011

Càlculo de distancias euclidianas a partir de landmarks en R.

NOTA: este post puede parecer muy largo, pero no lo es!!; lo que pasa es que se muestran los resultados para que sean comparados con los que obtendrás en casa al hacer la practica.

Muchas veces se necesitan realizar análisis de morfometrìa en los cuales aun se utilizan distancias euclidianas para llevar a cabo análisis de variaciones en tamaño como parte de estudios de tipo "tradicional".
Este tipo de análisis, a pesar de "boom" de el estudio de forma sin relación con el tamaño, esta aun vigente y es importante conocer aspectos básicos acerca del uso de técnicas de distancias "tradicionales".

Básicamente para el calculo de distancias euclidianas, se utiliza (valga la redundancia) la formula de distancias en el espacio euclidiano en un sistema de coordenadas .

Asi a partir de la toma de landmarks y los archivos que podamos tener de landmarks de varios individuos, podemos medir varias distancias en varios individuos con un solo Script, siguiendo una serie de comandos y ponerlos todos en un solo archivo. (los archivos todos deben haber sido tomados utilizando la misma escala (objeto de referencia de distancia conocida), puesto que estamos trabajando con tamaño)

De tal modo que empezare explicando el calculo de distancias en un solo individuo,(1) primero tomando solo una medida y después  (2)tomando varias medidas entre varios landmarks y poniéndolas todas en un un solo archivo.
(3)Finalmente les mostrare como se calculas varias distancias en varios individuos y ponerlas en un archivo como una matriz de distancias.

Primero necesitamos los archivos de los ladmarks en formato .txt, /aunque podria ser tambien .csv), llamaremos 7 archivos .txt de 7 individuos de Diachlorus curvipes a los cuales se les ha tomado los mismos landmarks (12), los archivos seran los siguientes y estaran sobre el directorio sobre el cual estemos trabajando en R.


dcp1.txt




"x" "y"
"lan1" 0.310384615384615 3.01932692307692
"lan2" 1.87769230769231 2.72125
"lan3" 1.96423076923077 3.01932692307692
"lan4" 5.06038461538461 2.58663461538461
"lan5" 5.13730769230769 2.41355769230769
"lan6" 4.40653846153846 1.75009615384615
"lan7" 4.07961538461538 1.57701923076923
"lan8" 3.06038461538462 2.05778846153846
"lan9" 3.16615384615385 1.33663461538461
"lan10" 2.07961538461538 1.32701923076923
"lan11" 2.06038461538462 1.37509615384615
"lan12" 1.74307692307692 2.25009615384615


dcp2.txt




"x" "y"
"lan1" 0.719286944871943 3.44351836802884
"lan2" 2.58291646916267 3.09794468140540
"lan3" 2.74336139509498 3.48054412016707
"lan4" 6.23612401346766 2.99920934237012
"lan5" 6.2854916829853 2.78939674692018
"lan6" 5.42155746642668 2.02419786939683
"lan7" 5.08832569718265 1.85141102608511
"lan8" 3.87881779400059 2.45616497767614
"lan9" 4.06394655469172 1.57988884373812
"lan10" 2.65696797343912 1.53052117422048
"lan11" 2.64462605605971 1.60457267849693
"lan12" 2.37310387371272 2.56724223409082


dcp3.txt




"x" "y"
"lan1" 1.10064827940554 3.7794857454421
"lan2" 2.84403223991721 3.48688983598560
"lan3" 2.92937271350869 3.86482621903358
"lan4" 6.42833213075936 3.30401739257528
"lan5" 6.44052362698672 3.18210243030174
"lan6" 5.57492739484456 2.43842116043313
"lan7" 5.29452298161541 2.29212320570488
"lan8" 4.18509682492618 2.86512352839053
"lan9" 4.31920328342707 1.97514430379367
"lan10" 3.07567066823694 1.91418682265689
"lan11" 3.03909617955487 1.97514430379367
"lan12" 2.6855427889616 2.99922998689142


dcp4.txt




"x" "y"
"lan1" 0.941808934714676 3.94869970156446
"lan2" 2.78204119586783 3.50996883135577
"lan3" 2.90391088203692 3.92432576433065
"lan4" 6.37719693785579 3.24185552178378
"lan5" 6.41375784370651 3.05905099253016
"lan6" 5.45098732297076 2.27908500104803
"lan7" 5.23162188786641 2.18158925211276
"lan8" 4.0494859320263 2.76656374572436
"lan9" 4.18354258681229 1.91347594254078
"lan10" 2.86734997618619 1.98659775424223
"lan11" 2.84297603895237 2.04753259732677
"lan12" 2.50174091767894 3.03467705529634

dcp5.txt

"x" "y"
"lan1" 0.585132420157185 3.30044582245595
"lan2" 2.00234100248879 2.90246259043132
"lan3" 2.09941008347041 3.18396292527801
"lan4" 4.87558579954464 2.49477245030853
"lan5" 4.90470652383913 2.34916882883611
"lan6" 4.19610223267333 1.79587506724089
"lan7" 3.93401571402296 1.68909907816111
"lan8" 3.02156635279576 2.27151356405081
"lan9" 3.10892852567922 1.56290927288501
"lan10" 2.03146172678327 1.65997835386663
"lan11" 2.02175481868511 1.68909907816111
"lan12" 1.78878902432923 2.53360008270118

dcp6.txt

"x" "y"
"lan1" 0.508475803224269 4.05301396510432
"lan2" 2.30040803313004 3.47574720647696
"lan3" 2.45675111359162 3.80045975820485
"lan4" 5.94440444696526 2.65795263175486
"lan5" 5.93237805616053 2.36931925244118
"lan6" 5.07850430902422 1.8521844478375
"lan7" 4.72973897568686 1.75597332139961
"lan8" 3.64736380326055 2.56174150531697
"lan9" 3.67141658487002 1.68381497657119
"lan10" 2.30040803313004 1.90029001105645
"lan11" 2.27635525152057 1.99650113749434
"lan12" 2.07190660784004 2.98266518348275

dcp7.txt

"x" "y"
"lan1" 0.137282309611791 3.15335927001144
"lan2" 1.72352748369761 2.92263269923532
"lan3" 1.87734519754835 3.2687225553995
"lan4" 5.06906275995132 2.84572384230995
"lan5" 5.12674440264535 2.68229252134353
"lan6" 4.32881501204461 1.99011280901518
"lan7" 4.07886122703715 1.86513591651145
"lan8" 3.0117508371976 2.40349791498906
"lan9" 3.18479576527969 1.62479573861966
"lan10" 1.95425405447372 1.51904606034727
"lan11" 1.92541323312671 1.59595491727264
"lan12" 1.61777780542522 2.45156595056741


Entonces...

>dcp1<-read.table ("dcp1.txt")

#llamo el archivo "dcp1.txt" de mi directorio y lo leo en R, poniéndole como nombre "dcp1" (Diachlorus curvipes, individuo 1).

> dcp1

#verifico que contiene mi nuevo objeto que acabo de llamar.

              x        y
lan1  0.3103846 3.019327
lan2  1.8776923 2.721250
lan3  1.9642308 3.019327
lan4  5.0603846 2.586635
lan5  5.1373077 2.413558
lan6  4.4065385 1.750096
lan7  4.0796154 1.577019
lan8  3.0603846 2.057788
lan9  3.1661538 1.336635
lan10 2.0796154 1.327019
lan11 2.0603846 1.375096
lan12 1.7430769 2.250096

>dist1_5<-list((sqrt(sum((dcp1[1,]-dcp1[5,])^2))))

#tomo la distancia entre el landmark 1  ("dcp1[1,]") de coordenadas  (x)0.3103846 (y)3.019327 y 5  ("dcp1[5,]") de coordenadas (X)5.1373077 (y)2.413558.

>dist1_5

#miro el valor de distancia que obtuve entre los dos puntos:

[[1]]
[1] 4.864786

> disdcp1<-list((sqrt(sum((dcp1[1,]-dcp1[5,])^2))),(sqrt(sum((dcp1[2,]-dcp1[8,])^2))),(sqrt(sum((dcp1[1,]-dcp1[11,])^2))),(sqrt(sum((dcp1[6,]-dcp1[7,])^2))),(sqrt(sum((dcp1[3,]-dcp1[10,])^2))),(sqrt(sum((dcp1[2,]-dcp1[12,])^2))),(sqrt(sum((dcp1[3,]-dcp1[4,])^2))),(sqrt(sum((dcp1[5,]-dcp1[10,])^2))),(sqrt(sum((dcp1[4,]-dcp1[6,])^2))))

#ahora tomo 9 distancias a partir de 12 landmarks, tomados de una foto. las distancias serian entre los landmarks 1y5, 2y8, 1y11, 6y7, 3y10, 2y12, 3y4, 5y10,   4y6. Y nombor este objeto "disdcp1".

> disdcp1

#reviso las distancias que obtuve:

[[1]]
[1] 4.864786

[[2]]
[1] 1.356076

[[3]]
[1] 2.401249

[[4]]
[1] 0.3699112

[[5]]
[1] 1.696237

[[6]]
[1] 0.4900074

[[7]]
[1] 3.126242

[[8]]
[1] 3.245004

[[9]]
[1] 1.061749

> distdcp1<-matrix(disdcp1, 1, 9)

#convierto mis resultados de la forma de lista a la forma de matriz!!!
> distdcp1

#reviso como se organizaron mis resultados en la nueva forma que los organice (matriz):

     [,1]     [,2]     [,3]     [,4]      [,5]     [,6]      [,7]     [,8]    
[1,] 4.864786 1.356076 2.401249 0.3699112 1.696237 0.4900074 3.126242 3.245004
     [,9]    
[1,] 1.061749

>write.table(distdcp1, file = "distdcp1.txt")

#escribo un archivo en formato .txt que contenga mi matriz de distancias halladas, que quedara sobre el directorio el cual estemos trabajando, al abrir este archivo en un edito de texto (geany, wordpad, blog de notas, etc) tendrá esta apariencia:

"V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9"
"1" 4.86478599236177 1.35607607002421 2.40124859656080 0.36991123195522 1.69623669778461 0.490007396393886 3.1262423565983 3.24500355576432 1.06174921263570

> dcp1.txt<-read.table ("dcp1.txt")
> dcp2.txt<-read.table ("dcp2.txt")
> dcp3.txt<-read.table ("dcp3.txt")
> dcp4.txt<-read.table ("dcp4.txt")
> dcp5.txt<-read.table ("dcp5.txt")
> dcp6.txt<-read.table ("dcp6.txt")
> dcp7.txt<-read.table ("dcp7.txt")

#leo de mi directorio 7 archivos que contienen los landmarks de 7 individuos de la misma especie (no necesariamente tienen que ser de la misma especie, valga la aclaración) 

>disdcp1al7<-list((sqrt(sum((dcp1.txt[1,]-dcp1.txt[5,])^2))),(sqrt(sum((dcp1.txt[2,]-dcp1.txt[8,])^2))),(sqrt(sum((dcp1.txt[1,]-dcp1.txt[11,])^2))),(sqrt(sum((dcp1.txt[6,]-dcp1.txt[7,])^2))),(sqrt(sum((dcp1.txt[3,]-dcp1.txt[10,])^2))),(sqrt(sum((dcp1.txt[2,]-dcp1.txt[12,])^2))),(sqrt(sum((dcp1.txt[3,]-dcp1.txt[4,])^2))),(sqrt(sum((dcp1.txt[5,]-dcp1.txt[10,])^2))),(sqrt(sum((dcp1.txt[4,]-dcp1.txt[6,])^2))),(sqrt(sum((dcp2.txt[1,]-dcp2.txt[5,])^2))),(sqrt(sum((dcp2.txt[2,]-dcp2.txt[8,])^2))),(sqrt(sum((dcp2.txt[1,]-dcp2.txt[11,])^2))),(sqrt(sum((dcp2.txt[6,]-dcp2.txt[7,])^2))),(sqrt(sum((dcp2.txt[3,]-dcp2.txt[10,])^2))),(sqrt(sum((dcp2.txt[2,]-dcp2.txt[12,])^2))),(sqrt(sum((dcp2.txt[3,]-dcp2.txt[4,])^2))),(sqrt(sum((dcp2.txt[5,]-dcp2.txt[10,])^2))),(sqrt(sum((dcp2.txt[4,]-dcp2.txt[6,])^2))),(sqrt(sum((dcp3.txt[1,]-dcp3.txt[5,])^2))),(sqrt(sum((dcp3.txt[2,]-dcp3.txt[8,])^2))),(sqrt(sum((dcp3.txt[1,]-dcp3.txt[11,])^2))),(sqrt(sum((dcp3.txt[6,]-dcp3.txt[7,])^2))),(sqrt(sum((dcp3.txt[3,]-dcp3.txt[10,])^2))),(sqrt(sum((dcp3.txt[2,]-dcp3.txt[12,])^2))),(sqrt(sum((dcp3.txt[3,]-dcp3.txt[4,])^2))),(sqrt(sum((dcp3.txt[5,]-dcp3.txt[10,])^2))),(sqrt(sum((dcp3.txt[4,]-dcp3.txt[6,])^2))),(sqrt(sum((dcp4.txt[1,]-dcp4.txt[5,])^2))),(sqrt(sum((dcp4.txt[2,]-dcp4.txt[8,])^2))),(sqrt(sum((dcp4.txt[1,]-dcp4.txt[11,])^2))),(sqrt(sum((dcp4.txt[6,]-dcp4.txt[7,])^2))),(sqrt(sum((dcp4.txt[3,]-dcp4.txt[10,])^2))),(sqrt(sum((dcp4.txt[2,]-dcp4.txt[12,])^2))),(sqrt(sum((dcp4.txt[3,]-dcp4.txt[4,])^2))),(sqrt(sum((dcp4.txt[5,]-dcp4.txt[10,])^2))),(sqrt(sum((dcp4.txt[4,]-dcp4.txt[6,])^2))),(sqrt(sum((dcp5.txt[1,]-dcp5.txt[5,])^2))),(sqrt(sum((dcp5.txt[2,]-dcp5.txt[8,])^2))),(sqrt(sum((dcp5.txt[1,]-dcp5.txt[11,])^2))),(sqrt(sum((dcp5.txt[6,]-dcp5.txt[7,])^2))),(sqrt(sum((dcp5.txt[3,]-dcp5.txt[10,])^2))),(sqrt(sum((dcp5.txt[2,]-dcp5.txt[12,])^2))),(sqrt(sum((dcp5.txt[3,]-dcp5.txt[4,])^2))),(sqrt(sum((dcp5.txt[5,]-dcp5.txt[10,])^2))),(sqrt(sum((dcp5.txt[4,]-dcp5.txt[6,])^2))),(sqrt(sum((dcp6.txt[1,]-dcp6.txt[5,])^2))),(sqrt(sum((dcp6.txt[2,]-dcp6.txt[8,])^2))),(sqrt(sum((dcp6.txt[1,]-dcp6.txt[11,])^2))),(sqrt(sum((dcp6.txt[6,]-dcp6.txt[7,])^2))),(sqrt(sum((dcp6.txt[3,]-dcp6.txt[10,])^2))),(sqrt(sum((dcp6.txt[2,]-dcp6.txt[12,])^2))),(sqrt(sum((dcp6.txt[3,]-dcp6.txt[4,])^2))),(sqrt(sum((dcp6.txt[5,]-dcp6.txt[10,])^2))),(sqrt(sum((dcp6.txt[4,]-dcp6.txt[6,])^2))),(sqrt(sum((dcp7.txt[1,]-dcp7.txt[5,])^2))),(sqrt(sum((dcp7.txt[2,]-dcp7.txt[8,])^2))),(sqrt(sum((dcp7.txt[1,]-dcp7.txt[11,])^2))),(sqrt(sum((dcp7.txt[6,]-dcp7.txt[7,])^2))),(sqrt(sum((dcp7.txt[3,]-dcp7.txt[10,])^2))),(sqrt(sum((dcp7.txt[2,]-dcp7.txt[12,])^2))),(sqrt(sum((dcp7.txt[3,]-dcp7.txt[4,])^2))),(sqrt(sum((dcp7.txt[5,]-dcp7.txt[10,])^2))),(sqrt(sum((dcp7.txt[4,]-dcp7.txt[6,])^2))))

#Tomo las mismas 9 medidas de cada uno de los siete individuos que tengo de la forma como lo hacíamos, con un solo individuo las nueve medidas, pero esta vez se hacen con los 7 individuos poniendo los comandos uno tras de otro. 


>disdcp1al7

#reviso las medidas que obtuve:

[[1]]
[1] 4.864786

[[2]]
[1] 1.356076

[[3]]
[1] 2.401249

[[4]]
[1] 0.3699112

[[5]]
[1] 1.696237

[[6]]
[1] 0.4900074

[[7]]
[1] 3.126242

[[8]]
[1] 3.245004

[[9]]
[1] 1.061749

[[10]]
[1] 5.604508

[[11]]
[1] 1.446113

[[12]]
[1] 2.662452

[[13]]
[1] 0.3753648

[[14]]
[1] 1.951936

[[15]]
[1] 0.5706719

[[16]]
[1] 3.525773

[[17]]
[1] 3.840697

[[18]]
[1] 1.270498

[[19]]
[1] 5.373187

[[20]]
[1] 1.478191

[[21]]
[1] 2.64825

[[22]]
[1] 0.3162748

[[23]]
[1] 1.956118

[[24]]
[1] 0.512768

[[25]]
[1] 3.543617

[[26]]
[1] 3.595809

[[27]]
[1] 1.215548

[[28]]
[1] 5.543798

[[29]]
[1] 1.469376

[[30]]
[1] 2.688656

[[31]]
[1] 0.2400554

[[32]]
[1] 1.938073

[[33]]
[1] 0.5517885

[[34]]
[1] 3.539701

[[35]]
[1] 3.705019

[[36]]
[1] 1.335961

[[37]]
[1] 4.423081

[[38]]
[1] 1.198715

[[39]]
[1] 2.158778

[[40]]
[1] 0.2830026

[[41]]
[1] 1.525499

[[42]]
[1] 0.4262206

[[43]]
[1] 2.860443

[[44]]
[1] 2.954745

[[45]]
[1] 0.9747592

[[46]]
[1] 5.67922

[[47]]
[1] 1.627789

[[48]]
[1] 2.711944

[[49]]
[1] 0.3617925

[[50]]
[1] 1.906591

[[51]]
[1] 0.5434545

[[52]]
[1] 3.67002

[[53]]
[1] 3.66213

[[54]]
[1] 1.182813

[[55]]
[1] 5.01165

[[56]]
[1] 1.388892

[[57]]
[1] 2.37127

[[58]]
[1] 0.2794568

[[59]]
[1] 1.751366

[[60]]
[1] 0.4827907

[[61]]
[1] 3.219626

[[62]]
[1] 3.379029

[[63]]
[1] 1.131387 


> distdcp1al7<-matrix(disdcp1al7, 9, 7)

#organizo mis medidas en una matriz de 9x7 (9 medidas correspondientes a 7 individuos).

> distdcp1al7

#ahora mis resultados se verán asi:

[,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]     
 [1,] 4.864786  5.604508  5.373187  5.543798  4.423081  5.67922   5.01165  
 [2,] 1.356076  1.446113  1.478191  1.469376  1.198715  1.627789  1.388892 
 [3,] 2.401249  2.662452  2.64825   2.688656  2.158778  2.711944  2.37127  
 [4,] 0.3699112 0.3753648 0.3162748 0.2400554 0.2830026 0.3617925 0.2794568
 [5,] 1.696237  1.951936  1.956118  1.938073  1.525499  1.906591  1.751366 
 [6,] 0.4900074 0.5706719 0.512768  0.5517885 0.4262206 0.5434545 0.4827907
 [7,] 3.126242  3.525773  3.543617  3.539701  2.860443  3.67002   3.219626 
 [8,] 3.245004  3.840697  3.595809  3.705019  2.954745  3.66213   3.379029 
 [9,] 1.061749  1.270498  1.215548  1.335961  0.9747592 1.182813  1.131387 


>write.table(distdcp1al7, file = "distdcp1al7.txt")

#por ultimo escribo mis resultados en una tabla en formato .txt para depues realizar los correspondientes analisis que desee. Este archivo se vera asi:

"V1" "V2" "V3" "V4" "V5" "V6" "V7"
"1" 4.86478599236177 5.60450803209322 5.37318672231976 5.54379830875036 4.42308131914733 5.67922032807822 5.01165001367009
"2" 1.35607607002421 1.44611252391904 1.47819063805238 1.46937649380901 1.19871472360978 1.62778876652987 1.38889176421965
"3" 2.40124859656080 2.66245224222281 2.64825004477358 2.68865630315053 2.15877795206031 2.71194445994619 2.37126981121088
"4" 0.36991123195522 0.375364762937514 0.316274764273176 0.240055441885765 0.283002570840330 0.361792535285579 0.279456827255025
"5" 1.69623669778461 1.95193578609124 1.95611787637005 1.93807289360278 1.52549859163370 1.90659073395059 1.75136598391422
"6" 0.490007396393886 0.570671895921929 0.512768012343564 0.551788472462498 0.426220596622018 0.54345449005379 0.48279071672441
"7" 3.1262423565983 3.52577280558723 3.54361729640782 3.53970078643772 2.86044316799144 3.67002020536417 3.21962558528581
"8" 3.24500355576432 3.84069681415740 3.59580942516654 3.70501885435476 2.95474519622204 3.66212980073973 3.37902905259225
"9" 1.06174921263570 1.27049834001773 1.21554781081638 1.33596082513795 0.97475918549504 1.18281250208974 1.13138718773756



De este modo, se pueden llegar a tomar muchísimas distancias de muchísimos individuos, haciendo los scripts y corriéndolos en R, y se obtienen mucho mas fácil,  confiables y rápido. 

que les rinda!!!!



miércoles, 16 de marzo de 2011

Análisis de Componentes Principales en R.

En esta ocasión me interesa un poco volver a lo básico y general, y alejarme un poco de cosas especificas para la Biología.
So, presento pues la forma básica de realizar el muy conocido "análisis de componentes principales" en R.
Básicamente este análisis reduce la dimensionalidad o cantidad de un conjunto de datos y como resultado se obtiene una serie de datos "virtuales" llamados componentes principales en los cuales esta contenida la mayor cantidad de variación de los datos primarios de entrada ordenados por importancia. 
Como un elemento adicional este análisis puede servir para explorar de una manera rapida (no excesivamente eficiente) como se agrupan las variables utilizadas y/o los individuos evaluados (mas información).


Para explicar rápidamente como realizar este análisis en R, utilizaremos un set de datos de medidas tomadas a varios organismos de distintas poblaciones (misma especie) de tortugas con tres variables (longitud (lon), ancho (anc), y altura (alt) del caparazón). 
Como siempre estos datos se encontraran en formato .CSV (separado por comas)
y debe tener la siguiente apariencia y al que llamaremos "tortugas_ninja.csv":




lon,anc,alt
98,81,38
103,84,38
103,86,42
105,86,42
109,88,44
123,95,46
134,100,48
136,102,49
123,92,50
133,99,51
133,102,51
133,102,51
138,98,51
138,99,51
141,105,53
149,107,55
153,107,56
147,108,57
155,117,60
158,155,62
155,115,63
159,118,63
162,124,61
177,132,67


Manos a la Obra:






>tortugas_ninja <- read.csv("tortugas_ninja.csv",header=T)

#(llamo a R el archivo que cree en formato .CSV (separado por comas) !!se puede crear en excel y guardarlo en formato .csv pero se deben cambiar los ";" por "," y los decimales deben estar separados por "." no por ","en home)


>tortugas_ninja


#(asi puedo ver lo que tiene mi archivo tortugas_ninja y garantizar que no hayan errores en el archivo)


   lon anc alt
1   98  81  38
2  103  84  38
3  103  86  42
4  105  86  42
5  109  88  44
6  123  95  46
7  134 100  48
8  136 102  49
9  123  92  50
10 133  99  51
11 133 102  51
12 133 102  51
13 138  98  51
14 138  99  51
15 141 105  53
16 149 107  55
17 153 107  56
18 147 108  57
19 155 117  60
20 158 155  62
21 155 115  63
22 159 118  63
23 162 124  61
24 177 132  67



>attach (tortugas_ninja)

#(para tomar mis columnas como objetos)


>tortugas.cp <- princomp(tortugas_ninja[,-1], scale = TRUE)

#(realizo mi análisis de componentes principales, el argumento "[,-1]" indica que no se debe tomar la primer columna del set de datos que corresponde a una variable  altamente correlacionada con las demás(longitud) (esto se sabe por previos análisis de correlación entre variables p.ej. un correlograma))


>summary (tortugas.cp)

#(inspecciono mis resultados del análisis de componentes realizado y que tanto explican mis componentes)

Importance of components:
                          Comp.1     Comp.2
Standard deviation     17.898105 3.39881269
Proportion of Variance  0.965194 0.03480607
Cumulative Proportion   0.965194 1.00000000


>plot (tortugas.cp)

#(gráfico las varianzas para los componentes explicativos encontrados)
>biplot(tortugas.cp)

# (realizo mi gráfico de componentes principales)

>loadings(tortugas.cp)

# (cargo los datos de las variables que mis componentes me están explicando)
Loadings:
    Comp.1 Comp.2
anc -0.915  0.404
alt -0.404 -0.915

               Comp.1 Comp.2
SS loadings       1.0    1.0
Proportion Var    0.5    0.5
Cumulative Var    0.5    1.0
>tortugas.cp$scores
# (de este modo obtengo mis scores de los componentes para cada individuo de la muestra inicial, para realizar posteriores análisis con estas nuevas variables)
          Comp.1      Comp.2
 [1,]  26.942316  3.44129819
 [2,]  24.198509  4.65436178
 [3,]  20.751886  1.80466153
 [4,]  20.751886  1.80466153
 [5,]  18.113972  0.78416594
 [6,]  10.903047  1.78544299
 [7,]   5.521326  1.97801099
 [8,]   3.287767  1.87211772
 [9,]  12.029436 -3.08602990
[10,]   5.222865 -1.17015052
[11,]   2.479058  0.04291307
[12,]   2.479058  0.04291307
[13,]   6.137467 -1.57450505
[14,]   5.222865 -1.17015052
[15,]  -1.073458 -0.57322800
[16,]  -3.711372 -1.59372359
[17,]  -4.115726 -2.50832592
[18,]  -5.434683 -3.01857372
[19,] -14.879168 -2.12318993
[20,] -50.442765 11.41307753
[21,] -14.263027 -5.67570597
[22,] -17.006834 -4.46264238
[23,] -21.685739 -0.20731055
[24,] -31.428684 -2.46008828




Y finalmente eso es todo amigos....
Buenos deseos!!!

viernes, 4 de marzo de 2011

Como manejo mis coordenadas? : Array's en R.

Retomando la Morfometrìa......

En el post de toma de landmarkas en R,  pudimos ver como obtener a partir de imágenes y/o fotografías la configuración espacial de un objeto o región X determinada...Pero, después de tener todos nuestros archivos en formato .txt en el que se encuentran almacenadas las configuraciones de los landmarks individuales de cada imagen.. como analizarlos en conjunto para poder hacer comparaciones??

El primer paso es "concatenar" mis archivos en un formato u objeto en R que sea manejable para realizar comparaciones y/o análisis, un tipo de objeto de este tipo es el objeto "array" que es la generalización de una matriz de dos dimensiones al caso multidimensional. (mas información).

Una vez tenemos nuestros datos almacenados en un "array" o varios "array's"; (dependiendo de sí es un "array" por cada especie o por cada sitio o etc... para ser comparados) podemos empezar con los distintos análisis morfometricos.

Entonces, supongamos que hemos tomado 12 landmarks en el ala de 4 individuos de Diachlorus ferrugatus  (Díptera: Tabanidae) y la configuracion geometrica del ala de esta especie sera comparada con la configuración de otras especies del mismo Genero, para ello necesitamos construir un "array" para cada especie. (Como ejemplo, utilizaremos el de la especie mencionada)

Partiríamos de 4 archivos .txt que contienen las coordenadas de los landmarks de los 4 individuos que utilizamos de la misma especie, estos tendrían una apariencia como esta:

Archivo 1 (df1.txt):


"x" "y"
"lan1" 0.987216506531992 4.65296820402765
"lan2" 3.44297896059167 4.30214499630484
"lan3" 3.72363752676992 4.81083864750291
"lan4" 9.33680885033491 4.40739195862168
"lan5" 9.58238509574088 3.98640410935431
"lan6" 8.02122182137436 2.75852288232447
"lan7" 7.56515165133471 2.53048779730464
"lan8" 5.72332981078995 3.44262813738395
"lan9" 6.0916941788989 2.14458226880955
"lan10" 3.98675493256203 2.00425298572042
"lan11" 3.75871984754220 2.33753503305709
"lan12" 3.21494387557185 3.56541626008693

Archivo 2 (df2.txt):


"x" "y"
"lan1" 0.451428571428571 5.23678571428571
"lan2" 3.09428571428571 4.80821428571429
"lan3" 3.29071428571429 5.3975
"lan4" 9.41571428571429 4.93321428571429
"lan5" 9.68357142857143 4.48678571428571
"lan6" 7.84428571428572 3.09392857142857
"lan7" 7.32642857142857 2.84392857142857
"lan8" 5.36214285714286 3.91535714285714
"lan9" 5.88 2.50464285714286
"lan10" 3.79071428571429 2.21892857142857
"lan11" 3.54071428571429 2.55821428571429
"lan12" 2.79071428571429 4.04035714285714

Archivo 3 (df3.txt):

"x" "y"
"lan1" 0.149340776528574 5.28933286531034
"lan2" 2.99908023202552 4.91057002628859
"lan3" 3.17944348870254 5.50576877332276
"lan4" 9.6003754264045 5.07289695729791
"lan5" 9.90699296275545 4.54984351293455
"lan6" 8.04925141898212 3.17908276218918
"lan7" 7.58030695162187 2.98068317984446
"lan8" 5.56023847683922 3.95464476590038
"lan9" 6.02918294419948 2.54781136381961
"lan10" 3.77464223573671 2.27726647880408
"lan11" 3.48606102505347 2.65602931782582
"lan12" 2.72853534700998 4.11697169690970

Archivo 4 (df4.txt):

"x" "y"
"lan1" 0.74741186771418 4.70685628896825
"lan2" 3.18693026280777 4.35330289837498
"lan3" 3.38138462763407 4.90131065379455
"lan4" 9.21501557242309 4.56543493273094
"lan5" 9.46250294583838 4.19420387260801
"lan6" 7.87151268816865 2.95676700553155
"lan7" 7.41189328039739 2.72695730164592
"lan8" 5.50270497119371 3.55780776954011
"lan9" 5.96232437896497 2.33804857199332
"lan10" 3.92939238305364 2.07288352904836
"lan11" 3.68190500963835 2.33804857199332
"lan12" 3.01015356751113 3.62851844765877

En estos archivos se presentan las coordenadas de 12 landmarks tomados en cada uno de los 4 individuos de Diachlorus ferrugatus que utilizamos, es por eso que los archivos se llaman "df1.txt", "df2.txt", "df3.txt" y "df4.txt" (por las iniciales de la especie y el numero del individuo (pero se puede utilizar el nombre que se desee)) respectivamente.
Los diminutivos "lan1" a "lan12" indican el numero del landmark y se presentan las coordenadas en los ejes "x" y "y".
Estos archivos deben estar contenidos en el directorio sobre el cual estamos trabajando en R!!!

Para hacer el array necesitamos un script que puede estar en formato ".txt" o en formato ".R" así (para este caso) o escribir todos los comandos uno por uno!!:

 a<-read.table("df1.txt")
 b<-read.table("df2.txt")
 c<-read.table("df3.txt")
 d<-read.table("df4.txt")
 a<-as.matrix(a)
 b<-as.matrix(b)
 c<-as.matrix(c)
 d<-as.matrix(d)
 Dferrugatus<-array(cbind(a, b, c, d), dim=c(12, 2, 4 ))

#En este Script leemos primero cada uno de de nuestros archivos y les damos nombres: a, b, c y d. :

> a <- read.table ("df1.txt")
> b <- read.table ("df2.txt")
> c <- read.table ("df3.txt")
> d <- read.table ("df4.txt")

#Después los tomamos como matrices:
 
>a<-as.matrix(a)
>b<-as.matrix(b)
>c<-as.matrix(c)
>d<-as.matrix(d)

#Y finalmente construimos el array utilizando la función "cbind" y diciéndole cuales objetos utilizar (a,b,c,d) y dándole las dimensiones del array (12= que es el numero de filas sin contar los encabezados y en este caso es el numero de landmarks), (2= que es el numero de columnas que en este caso es el numero de variables "x" y "y"), y (4= que es el numero de archivos que estamos utilizando que este caso es 4):

>Dferrugatus<-array(cbind(a, b, c, d), dim=c(12, 2, 4 ))

#Podemos revisar nuestro array dandole el nombre de este que en este caso es 
"Dferrugatus":

>Dferrugatus

Y lo que obtendríamos (array) seria algo así: 

, , 1

           [,1]     [,2]
 [1,] 0.9872165 4.652968
 [2,] 3.4429790 4.302145
 [3,] 3.7236375 4.810839
 [4,] 9.3368089 4.407392
 [5,] 9.5823851 3.986404
 [6,] 8.0212218 2.758523
 [7,] 7.5651517 2.530488
 [8,] 5.7233298 3.442628
 [9,] 6.0916942 2.144582
[10,] 3.9867549 2.004253
[11,] 3.7587198 2.337535
[12,] 3.2149439 3.565416

, , 2

           [,1]     [,2]
 [1,] 0.4514286 5.236786
 [2,] 3.0942857 4.808214
 [3,] 3.2907143 5.397500
 [4,] 9.4157143 4.933214
 [5,] 9.6835714 4.486786
 [6,] 7.8442857 3.093929
 [7,] 7.3264286 2.843929
 [8,] 5.3621429 3.915357
 [9,] 5.8800000 2.504643
[10,] 3.7907143 2.218929
[11,] 3.5407143 2.558214
[12,] 2.7907143 4.040357

, , 3

           [,1]     [,2]
 [1,] 0.1493408 5.289333
 [2,] 2.9990802 4.910570
 [3,] 3.1794435 5.505769
 [4,] 9.6003754 5.072897
 [5,] 9.9069930 4.549844
 [6,] 8.0492514 3.179083
 [7,] 7.5803070 2.980683
 [8,] 5.5602385 3.954645
 [9,] 6.0291829 2.547811
[10,] 3.7746422 2.277266
[11,] 3.4860610 2.656029
[12,] 2.7285353 4.116972

, , 4

           [,1]     [,2]
 [1,] 0.7474119 4.706856
 [2,] 3.1869303 4.353303
 [3,] 3.3813846 4.901311
 [4,] 9.2150156 4.565435
 [5,] 9.4625029 4.194204
 [6,] 7.8715127 2.956767
 [7,] 7.4118933 2.726957
 [8,] 5.5027050 3.557808
 [9,] 5.9623244 2.338049
[10,] 3.9293924 2.072884
[11,] 3.6819050 2.338049
[12,] 3.0101536 3.628518


Cabe recordar que estos array no son utilizados solamente para morfometrìa sino para múltiples análisis en R y que su construcción es importante para todo usuario de R!!!!
(mas adelante utilizaremos estos array)

Sean felices y hasta la próxima !!!!