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.














2 comentarios:

  1. Excelente aporte muchisimas gracias por el trabajo me saco de varios apuros ya que en la red es difícil encontrar a Yule

    ResponderEliminar
  2. Muchas gracias por su comentario y me alegra mucho que le haya gustado el post ya que fue con mucha dedicacion!

    ResponderEliminar