Ejercicio modelos Log-Lineales
Sandra Katteryne Rodríguez Hurtado y Esmeralda Lagos Moreno
14/12/2020
Se quiere hacer un estudio interacción de presencia de hormigas, áfidos y una enfermedad presente en plantas en el municipio Nuevo Progreso - Meta.
Tabla de datos
Tabla de contingencia: tabcontin
Funciones a utilizar: -Función array: Tabla con más de dos dimensiones. -Función list: Para identificar los niveles de cada variable categórica. -Función data: Proporciona el conteo de celdas. -dinames: Especificar los nombres. -Argumento dim: Crear una tabla con 2 filas, 2 columnas y 2 capas (tabla 2x2). Tres variables con dos niveles cada una.
1. Subir los datos a R
Las dimensiones son: Sí, No.
Presencia de hormigas vs la presencia de áfidos para cada nivel de presencia de la enfermedad.
tabcontin <- array(data = c(180,10, 108,90,50, 2, 12, 20),
dim = c(2,2,2),
dimnames = list("Hormigas" = c("S","N"),
"Afido" = c("S","N"),
"Enfermedad" = c("S","N")))
tabcontin## , , Enfermedad = S
##
## Afido
## Hormigas S N
## S 180 108
## N 10 90
##
## , , Enfermedad = N
##
## Afido
## Hormigas S N
## S 50 12
## N 2 20Usar la función ftable para aplanar la tabla con el fin de que se parezca a la fuente de datos. Argumento row.vars especifica qué variables se quieren en las filas y el orden en que se quieren, para este caso será Enfermedad y Hormiga.
ftable(tabcontin, row.vars = c("Enfermedad","Hormigas"))## Afido S N
## Enfermedad Hormigas
## S S 180 108
## N 10 90
## N S 50 12
## N 2 201. Explorar los datos
Función addmargins: Proporciona totales marginales. En este caso, 472 plantas fueron analizadas, lo anterior a la sumatoria de los valores dispuestos en la tabla.
addmargins(tabcontin)## , , Enfermedad = S
##
## Afido
## Hormigas S N Sum
## S 180 108 288
## N 10 90 100
## Sum 190 198 388
##
## , , Enfermedad = N
##
## Afido
## Hormigas S N Sum
## S 50 12 62
## N 2 20 22
## Sum 52 32 84
##
## , , Enfermedad = Sum
##
## Afido
## Hormigas S N Sum
## S 230 120 350
## N 12 110 122
## Sum 242 230 472La relación de plantas enfermas con respecto a las sanas es de (5:1).
Función prop.table: Permite calcular las proporciones de los valores.
options(digits=3)
prop.table(tabcontin, margin = c(1,3))## , , Enfermedad = S
##
## Afido
## Hormigas S N
## S 0.625 0.375
## N 0.100 0.900
##
## , , Enfermedad = N
##
## Afido
## Hormigas S N
## S 0.8065 0.194
## N 0.0909 0.909Según los resultados, se puede apreciar que en las plantas enfermas que tienen hormigas el 62.5% presenta áfidos y de las plantas enfermas sin hormigas el 10% tiene áfidos.
El 80% de las plantas sanas que presentan hormigas, presentan afidos. Por su parte, el 90% de las plantas sanas que no tienen hormigas tampoco tienen afidos.
Modelado Loglineal
- Tabla de frecuencia
Uso de la función glm (modelos lineales generalizados), se debe transformar los datos en un formato diferente. Para esto, lo primero que se debe hacer es transformar los datos en un data.frame por medio de la función as.data.frame, y se convierte en una tabla con la función as.table. Se tiene una tabla de contingencia 2^3.
tabcontin.df <- as.data.frame(as.table(tabcontin))
tabcontin.df[,-4] <- lapply(tabcontin.df[,-4], relevel, ref = "N")
tabcontin.df## Hormigas Afido Enfermedad Freq
## 1 S S S 180
## 2 N S S 10
## 3 S N S 108
## 4 N N S 90
## 5 S S N 50
## 6 N S N 2
## 7 S N N 12
## 8 N N N 20De 472 plantas analizadas, 20 de estas no estaban enfermas y no presentaron hormigas ni áfidos; 180 estan enfermas y presentan hormigas y áfidos; 2 presentan áfidos, pero no estan ni presentan hormigas; 108 estan enfermas y presentan hormigas, pero no áfidos.
2. Modelo de independencia Loglineal.
Modelo Loglineal: Permiten determinar si dos variables están asociadas y de qué manera. La tabla de contingencia que se trabajará en esta ocasión será de tres vías. Como se está modelando conteos, en el argumento family se pone poisson.
Función glm: Modela la frecuencia como una función de las tres variables (enfermedad, hormigas y áfidos).
mod0 <- glm(Freq ~ Hormigas + Afido + Enfermedad,
data=tabcontin.df, family = poisson)
summary(mod0)##
## Call:
## glm(formula = Freq ~ Hormigas + Afido + Enfermedad, family = poisson,
## data = tabcontin.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 2.58 -7.08 -2.83 5.26 2.95 -3.38 -3.80 2.58
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3590 0.1422 16.59 <2e-16 ***
## HormigasS 1.0539 0.1051 10.02 <2e-16 ***
## AfidoS 0.0509 0.0921 0.55 0.58
## EnfermedadS 1.5302 0.1203 12.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 461.07 on 7 degrees of freedom
## Residual deviance: 133.63 on 4 degrees of freedom
## AIC: 183.2
##
## Number of Fisher Scoring iterations: 5El valor arrojado de criterio de información (AIC) es 183.2. Con el valor de la desviación residual (133.63) hay un problema y es que se presenta sobredispersión, es decir que tiene más variabilidad que la esperada en el modelo de Poisson. Según lo anterior, en los modelos loglineales se busca calcular la desviación residual y que el resultado sea cercano a los grados de libertad.
Pruebas para el modelo.
- Calcular p_valor.
La estadística de deviancia residual tiene una distribución aproximada de Chi-cuadrado, por lo que se usa la función pchiq y la función df.residual: Extrae los grados de libertad.
pchisq(deviance(mod0), df = df.residual(mod0), lower.tail = F)## [1] 6.53e-28Cuando el p_valor es un valor cercano a cero indica que el modelo es malo debido a que se rechaza la Ho de que las frecuencias esperadas satisfacen el modelo loglineal dado.
- Comparar valores ajustados con observados.
cbind(mod0$data, fitted(mod0))## Hormigas Afido Enfermedad Freq fitted(mod0)
## 1 S S S 180 147.5
## 2 N S S 10 51.4
## 3 S N S 108 140.2
## 4 N N S 90 48.9
## 5 S S N 50 31.9
## 6 N S N 2 11.1
## 7 S N N 12 30.4
## 8 N N N 20 10.6El modelo ajustado predice que 51.4 plantas están enfermas y presentan áfidos, pero no tienen hormigas; sin embargo el valor observado corresponde a 10 plantas. Del mismo modo, el modelo predice que 10 plantas no estan enfermas y no presentan hormigas ni áfidos, pero el valor observado es de 20 plantas.
library(ggplot2)
ggplot(mod0$data,aes(fitted(mod0),Freq)) +
geom_point() +
labs(x = "Ajustados", y = "Observados") +
geom_abline(intercept = 0, slope = 1, color = "blue") +
theme_bw()
Los valores estimados para el modelo ajustado difiren demasiado de los valores originales observados. A partir de lo anterior se confirma que el modelo no es adecuado.
- Análisis odds
Es posible exponenciar los coeficientes para que puedan ser interpretados como probabilidades. Se debe aclarar que el análisis odds no debe confundirse con probabilidades, son posibilidades o también son conocidas como coeficiente de dispariedad. Su fórmula es la siguiente, que puede pensarse también como (Éxito / Fracaso) o (Sí/No):
\[odds=\frac{p}{1-p}\]
Se puede exponenciar el coeficiente de presencia de áfidos, el cual arrojó un resultado de 1.05, de la sigueinte manera:
exp(coef(mod0)[3])## AfidoS
## 1.05Este valor de 1.05 significa que existe una relación 1:1 de presentar áfidos, en presencia de hormigas y cuando las plantas estan enfermas, pero cuando se observa la tabla de frecuencia original se observa que este valor no se ajusta a la realidad, puesto que la relación es de 180/108, aproximadamente 2:1. Este es otro indicador que informa que el modelo no es bueno.
3. Modelo de asociación homogénea
Se le llama asociación homogénea al ajustar un modelo más complejo que permita la asociación de las variables entre sí. La relación que existe entre dos variables es idependiente de la existencia de una tercera variable, por tanto se pueden hacer combinaciones de variables por pares.
mod1 <- glm(Freq ~ (Hormigas + Afido + Enfermedad)^2,
data = tabcontin.df, family = poisson)
summary(mod1)##
## Call:
## glm(formula = Freq ~ (Hormigas + Afido + Enfermedad)^2, family = poisson,
## data = tabcontin.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.116 0.522 0.151 -0.164 0.223 -0.902 -0.432 0.358
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.915 0.222 13.15 < 2e-16 ***
## HormigasS -0.307 0.315 -0.97 0.330
## AfidoS -1.645 0.373 -4.41 1.1e-05 ***
## EnfermedadS 1.602 0.240 6.67 2.6e-11 ***
## HormigasS:AfidoS 2.918 0.329 8.88 < 2e-16 ***
## HormigasS:EnfermedadS 0.458 0.332 1.38 0.167
## AfidoS:EnfermedadS -0.739 0.297 -2.49 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 461.0688 on 7 degrees of freedom
## Residual deviance: 1.5132 on 1 degrees of freedom
## AIC: 57.11
##
## Number of Fisher Scoring iterations: 4Según la tabla anterior, se puede apreciar que la deviancia residual tiene un valor de 1,51 el cual no difiere de los grados de libertad por lo que este modelo ya no tiene problemas de sobredispersión.
Pruebas para el modelo.
- Calcular el p-valor de la siguiente manera:
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 0.219El p_valor obtenido no rechaza la Ho de que las frecuencias esperadas satisfacen el modelo loglineal dado, por tanto se puede decir que el modelo es bueno.
2.Comparar valores ajustados con observados
Ahora se volverá a comparar los valores del modelo ajustado y los valores observados, dando como resultado la siguiente tabla:
cbind(mod1$data, fitted(mod1))## Hormigas Afido Enfermedad Freq fitted(mod1)
## 1 S S S 180 181.56
## 2 N S S 10 8.44
## 3 S N S 108 106.44
## 4 N N S 90 91.56
## 5 S S N 50 48.44
## 6 N S N 2 3.56
## 7 S N N 12 13.56
## 8 N N N 20 18.44Se puede apreciar que ahora los valores ajustados y observados se parecen mucho más. Por ejemplo, las plantas sanas que no presentan hormigas ni áfidos tienen un valor 18.44 en el modelo ajustado y en el observado es de 20. Así mismo, la presencia de áfidos y hormigas, en plantas sanas, el valor observado es de 50 y el ajustado de 48 plantas aproximadamente.
library(ggplot2)
ggplot(mod1$data,aes(fitted(mod1),Freq)) +
geom_point() +
labs(x = "Ajustados", y = "Observados") +
geom_abline(intercept = 0, slope = 1, color = "orange") +
theme_bw()
Este gráfico indica que este es un modelo ajustado con buenas predicciones, debido a que los valores ajustados en el modelo concuerdan con los valores observados.
Coeficiente de interacciones para los modelos loglineal de asociación homogenea.
Sirve para describir la asociación entre las variables y el efecto que tienen entre sí. Lo anteior se hace exponenciando los coeficientes.
Odds<- exp(mod1$coefficients);Odds## (Intercept) HormigasS AfidoS
## 18.440 0.735 0.193
## EnfermedadS HormigasS:AfidoS HormigasS:EnfermedadS
## 4.965 18.504 1.581
## AfidoS:EnfermedadS
## 0.477Odds_df <- data.frame(Odds)La posibilidad de que halla presencia de hormigas frente a la presencia de áfidos es de 18.5, lo cual indica que por cada 19 plantas que presentan hormigas, 1 presenta áfidos. Por cada 2 plantas que tienen hormigas, una esta esferma; y por cada planta que presenta áfidos, 2 estan enfermas.
Intervalos de confianza para los modelos loglineal de asociación homogenea.
Cálculo de intervalos de confianza para las estimaciones de odds ratios (posibilidades) con la función confint.
inter_confianza<- exp(confint(mod1, parm = c("HormigasS:AfidoS", "HormigasS:EnfermedadS", "AfidoS:EnfermedadS")))## Waiting for profiling to be done...data.frame(inter_confianza)## X2.5.. X97.5..
## HormigasS:AfidoS 10.092 36.965
## HormigasS:EnfermedadS 0.822 3.037
## AfidoS:EnfermedadS 0.262 0.844A partir de la información anterior se puede concluir que la posibilidad de que una planta presente áfidos cuando presenta hormigas es 10 a 37 veces mayor, lo cual indica la presencia de áfidos tiene un fuerte relación con la presencia de hormigas. Casualmente estos resultados coninciden con una relación mutualista que existe entre estos dos insectos, en la cual las hormigas protegen a los áfidos contra depredadores a cambio de secreciones azucaradas que utilizan para el cultivo de hongos.
La posibilidad de que las plantas esten enfermas es de 0.82 a 3.037 veces más alta cuando las plantas presentan hormigas, y de 0.26 a 0.84 veces más alta cuando presentan áfidos. Por lo tanto, se puede concluir que la enfermedad de las plantas tiene una relación mas fuerte con la presencia de hormigas que de áfidos.
Ecuación final
\[log (\mu_{ijk})=log (n\pi_{i}\pi_{j}\pi_{k})= log n + log \pi_{i} + log \pi_j + log \pi_k + \epsilon_{ijk}\]