Coeficiente de asociación de Yule (calculo Bayesiano en R)

Recientemente Helena Brentani y colaboradores han introducido una rutina en R para calcular el coeficiente de asociación de G. Udny Yule de 1912, el cual mide la asociación entre dos atributos en una tabla de contingencia 2X2; este coeficiente de Yule es una modificación del método de "razón de productos cruzados" (CPR), el cual fue criticado por Yule, debido a la falta de limitación de los valores de asociación y al grado de desbalance que producía.

De este modo tenemos dos atributos en la tabla 1, "g" y "G" es un sistema de clasificación binario (por ejemplo: g= cuaderno grapado y G= cuaderno cosido), mientras que "t" y "T" son otro sistema de clasificación binario (por ejemplo: t= hojas cuadriculadas y T= hojas rayadas), y con el (CPR) se pretende medir el grado de asociación entre estos dos sistemas de clasificación (en el ejemplo nuestro: "que tan asociado están la forma en que se pegan las hojas con los trazos de las hojas del cuaderno".

CPR: en donde pi1 es la suma de n1 + p1, y asi sucesivamente para todos los valores.



El problema de esta medida de asociación es que no esta limitada ya que puede llegar hasta infinito y que esta ademas desbalanceada en sus resultados; Es por esta razon que Yule reformulo esta ecuacion para obtener es coeficiente de asociacion lambda (el que calcularemos en este post).

Coeficiente de asociación (lambda) de Yule.

Este coeficiente es limitado y balanceado, va de -1 a 1, (-1 a 0= asociación negativa; 0 a 1= asociación positiva; 0=no asociación = independencia).

Entonces, tenemos los anteriores datos, y queremos saber si se encuentran asociadas las dos variables anteriores.

Empezamos:

> set.seed(456)

> m<-10000000

#Introducimos el rotulo de los numeros generados al azar (456) y el tamaño de la muestra de Monte Carlo.

> a1<-0.5

> a12<-0.5

> a21<-0.5

> a3<-0.5

#Damos los valores de los parametros a priori, segun la regla de Jeffrey de 1931, en donde todos los valores son iguales (1/2).

> n1<-42

> n12<-12

> n21<-10

> n3<-35

#Introducimos los valores de las frecuencias observadas en nuestros datos (por ejemplo= 42 cuadernos cuadriculados y grapados).

> q1<-rgamma(m,a1+n1,2)

> q2<-rgamma(m,a12+n12,2)

> q3<-rgamma(m,a21+n21,2)

> q4<-rgamma(m,a3+n3,2)

#Genero mis valores a posteriori, dados mis datos (n) y mis parámetros (a), que haya incluido, utilizando el numero de permutaciones (m), utilizando la distribución gamma, que es una función de distribución de probabilidades con una función de densidad determinada.  

> lambda<- ((q1*q4)^.5 - (q2*q3)^.5)/((q1*q4)^.5 + (q2*q3)^.5)

> lambda<-c(-1,lambda,1)

#calculo el Coeficiente de asociación de lambda para todas las permutaciones utilizadas y realizo un vector entre -1, los valores de lambda y 1.

> freq<-hist(lambda,breaks=200,main="Posterior of lambda",probability=T,ylab="Density")

#Realizo un histograma de frecuencias para ver mis resultados graficamente.

> count_pos<-sum(freq$count[1:100])

> position_pos<-200

>while (position_pos>=101){

   if (freq$count[position_pos]<=freq$count[100])count_pos<-count_pos+freq$count[position_pos]
   position_pos<-position_pos-1
}

#Defino cuando mi función tiene un comportamiento de asociación positivo, utilizando mis datos de frecuencias del histograma.

> count_neg<-sum(freq$count[101:200])

> position_neg<-1

>while (position_neg<=100){

   if (freq$count[position_neg]<=freq$count[101])count_neg<-count_neg+freq$count[position_neg]
   position_neg<-position_neg+1
}

#Defino cuando mi función tiene un comportamiento de asociación negativo, utilizando mis datos de frecuencias del histograma.

> evalue_pos<-(count_pos-2)/m
> evalue_neg<-(count_neg-2)/m
## Defino mis valores de probabilidad "e" (que es el área bajo la curva del histograma de densidad posterior que se encuentra (o se pasa) mas allá del valor de Lambda=0), se utiliza "evalue_pos" sí esta el histograma hacia los valores positivos de lambda y "evalue_neg" sí esta el histograma en los valores negativos de lambda.

> alpha<-0.05
# (1-alpha) es la credibilidad, en este caso del 95%.

>n.HPD<-m*alpha 
>LI<-1
LS<-200 
>count.HPD<-freq$count[LI]+freq$count[LS] 
>while (count.HPD < n.HPD){
   if (freq$count[LI]<=freq$count[LS])
      {LI<-LI+1
       count.HPD<-count.HPD+freq$count[LI]
      }
   if (freq$count[LI]> freq$count[LS])
      {LS<-LS-1
       count.HPD<-count.HPD+freq$count[LS]
      }
}
# Para definir el intervalo fidedigno o de "confianza", dados los valores de frecuencia.

> mean(lambda)
[1] 0.5498831
# Via promedio, obtengo el valor del coeficiente de asociación, estimado por el metodod Bayesiano.


> (LI-100)/100
 [1] 0.39 
> (LS-100)/100
 [1] 0.71
# Obtenemos los limites fidedignos mas bajo (LI-100)/100  y  mas alto (LS-100)/100.

> evalue_pos
[1] -1e-07
> evalue_neg
[1] 0.9999999
# Evalúa la probabilidad "evalue_pos" sí lambda es positivo, "evalue_neg" sì lambda es negativo.

Eso es todo, nos vemos al siguiente.

Mas allá de Hardy-Weinberg: Coeficiente de desequilibrio, una perspectiva Bayesiana en R!

Recientemente Helena Brentani y Colaboradores de la Universidad de Sao Paulo y de Brasilia propusieron un nuevo coeficiente de medida de desviación del equilibrio de Hardy-Weinberg, basado en métodos Bayesianos y modificado del coeficiente propuesto por Pereira & Rogatko en 1984, en el articulo de Brentani et al. 2011, ellos introducen dos rutinas para calcular el coeficiente de desequilibrio planteado en el articulo y el coeficiente de asociación de Yule en 1912.
En el presente post, explicare con tanto detalle como pueda la forma de seguir la rutina al hacer el calculo del coeficiente de desequilibrio, la rutina utilizada puede ser revisada en la pagina de uno de los autores (Carlos Alberto Pereira).

>set.seed(1234) 

>m<-10000000

# Para empezar, tenemos que generar un numero de permutaciones o el tamaño (m) de la muestra de Monte Carlo, que sea reproducible, y utilizamos "set.seed", para rotular (1234), la operación de generación de números al azar, de este modo, al volver a generar los números al azar y anteponer "set.seed (1234), se generaran los mismos, números (mas información).

>a1<-0.5    
>a2<-0.5    
>a3<-0.5    

# Escojo los parámetros a priori (los cuales en el articulo se sugiere que sean escogidos según la regla de Jeffreys de 1931; (a1, a2, a3=1/2).

>n1<-4    
>n2<-18    
>n3<-94

# Incluyo los valores de las frecuencias absolutas de los genotipos del gen observado; n1=AA (homocigoto dominante), n2=Aa (Heterocigoto) y 

n3=aa (Homocigoto recesivo).

>p1<-rgamma(m,a1+n1,2)
>p2<-rgamma(m,a2+n2,2)
>p3<-rgamma(m,a3+n3,2)

# Genero mis valores a posteriori, dados mis datos (n) y los parámetros
 (a) que haya incluido, utilizando el numero de permutaciones (m) que 
había incluido, utilizando la distribución gamma, que es una función de distribución de probabilidades, con una función de densidad determinada (mas información).
>theta<-4*p1*p3/(p2^2)
# Calculo el valor de theta de acuerdo a la siguiente formula, pero aplicada para valores de genotipos en donde pi12 y pi21 corresponden a los individuos heterocigotos, por eso son cambiados por pi2="p2": 


>lambda<- (theta^0.5 - 1)/(theta^0.5 + 1)

# Finalmente calculamos el coeficiente de desequilibrio (lambda) a partir del valor de theta que obtuvimos (para mas infomracion, ver el articulo), segun la formula:


# La unión de estas dos formulas y su adaptación a la genética de poblaciones resultaría en la formula general, planteada por los autores, como formula general de calculo de coeficiente de desequilibrio:

>lambda<-c(-1,lambda,1)

# Concateno, todos los valores de lambda que obtuve en todas las permutaciones que genere al azar y las pongo entre -1 y 1 (los limites de correlación del coeficiente).

>freq<-hist(lambda,breaks=200,main="Posterior of lambda",probability=T,ylab="Density")

# Genero un histograma de frecuencias, para realizar una interpretación de mis resultados mas gráfica, en el que se divide el histograma en 200 columnas y se rotulan, los ejes y el gráfico.
# En el diagrama podemos observar que este se encuentra por fuera del equilibrio HW, ya que la densidad máxima de la distribución posterior se haya muy alejada del valor de Lambda=0=Equilibrio HW.
>count_pos<-sum(freq$count[1:100])
>position_pos<-200
>while (position_pos>=101){
   if (freq$count[position_pos]<=freq$count[100])count_pos<-count_pos+freq$count[position_pos]
   position_pos<-position_pos-1
}

# Defino cuando mi funcion tiene un comportamiento de asociación positivo (Sensu Brentani et al, 2011), utilizando mis datos de frecuencias de histogramas!!!

>count_neg<-sum(freq$count[101:200])
>position_neg<-1
>while (position_neg<=100){
   if (freq$count[position_neg]<=freq$count[101])count_neg<-count_neg+freq$count[position_neg]
   position_neg<-position_neg+1
}

# Y realizo lo mismo para definir sì es negativa la asociación o la desviación de HW, mas exactamente!!

>evalue_pos<-(count_pos-2)/m
>evalue_neg<-(count_neg-2)/m
# Defino mis valores de probabilidad "e" (que es el área bajo la curva del histograma de densidad posterior que se encuentra (o se pasa) mas allá del valor de Lambda=0), se utiliza "evalue_pos" sí esta el histograma hacia los valores positivos de lambda y "evalue_neg" sí esta el histograma en los valores negativos de lambda.
>alpha<-0.05
# (1-alpha) es la credibilidad, en este caso del 95%.
>n.HPD<-m*alpha 
>LI<-1
>LS<-200
count.HPD<-freq$count[LI]+freq$count[LS]
while (count.HPD < n.HPD){
   if (freq$count[LI]<=freq$count[LS])
      {LI<-LI+1
       count.HPD<-count.HPD+freq$count[LI]
      }
   if (freq$count[LI]> freq$count[LS])
      {LS<-LS-1
       count.HPD<-count.HPD+freq$count[LS]
      }
}
# Para definir el intervalo fidedigno o de "confianza", dados los valores de frecuencia.
>mean(lambda) 
# Via promedio, obtengo el valor del coeficiente de desequilibrio, estimado por el metodod Bayesiano.
>(LI-100)/100  
>(LS-100)/100  

# Obtenemos los limites fidedignos mas bajo (LI-100)/100  y  mas alto (LS-100)/100.
>evalue_pos  
>evalue_neg  
# Evalúa la probabilidad "evalue_pos" sí lambda es positivo, "evalue_neg" sì lambda es negativo.
Eso es todo, y espero les haya gustado este post, ya que es bastante actual y tiene varias cosas interesantes!!!

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!!!!

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!!!!

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!!!