jueves, 16 de junio de 2011

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.














miércoles, 8 de junio de 2011

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