Universidad Nacional de Colombia
Facultad de Ciencias Agrarias
Diseño de experimentos
Taller 4. Modelos Log-lineales
Introducción a modelos Log-lineales
Los modelos Log-lineales, también denominados “modelos lineales logarítmicos” y “modelos log-linear”, se presentan como una técnica que analiza la relación que se produce entre un conjunto de variables cualitativas (siempre más de dos, pues de ser éste el caso el análisis de tablas de contingencia). En consecuencia, los modelos log-lineales se presentan como la técnica más apropiada en aquellos casos en los que nos interese valorar la relación que se produce entre las variables de una tabla de contingencia de múltiples entradas.
Tomado de: [Análisis de tablas de Contingencia]Ejercicio
A continuación, se plantea un escenario donde se analizan los resultados a través de un modelo Log-lineal. Partimos de dos variedades de Maracuyá donde se evalúa la presencia de enfermedades y plagas.
I. Datos
Variable I: Variedad de Maracuyá. Passiflora edulis flavicarpa ‘PEF’ (Si) y Passiflora edulis f. edulis ‘PEFE’ (No)
Variable II: Presencia de la enfermedad Alternaria passiflorae
Variable III: Presencia de plagas, nematodos. Meloidogyne sp
datos=array(data=c(190,20,118,100,60,12,22,30),dim=c(2,2,2),dimnames=list('Alternaria passiflorae'=c('Si','No'),'Meloidogyne sp'=c('Si','No'),'Variedad PEF'=c('Si','No')));datos## , , Variedad PEF = Si
##
## Meloidogyne sp
## Alternaria passiflorae Si No
## Si 190 118
## No 20 100
##
## , , Variedad PEF = No
##
## Meloidogyne sp
## Alternaria passiflorae Si No
## Si 60 22
## No 12 30A continuación, se hace un arreglo para organizar los datos en una distribución cómoda para el lector.
ftable(datos,row.vars = c('Alternaria passiflorae','Meloidogyne sp'))## Variedad PEF Si No
## Alternaria passiflorae Meloidogyne sp
## Si Si 190 60
## No 118 22
## No Si 20 12
## No 100 30II. Análisis descriptivo
addmargins(datos)## , , Variedad PEF = Si
##
## Meloidogyne sp
## Alternaria passiflorae Si No Sum
## Si 190 118 308
## No 20 100 120
## Sum 210 218 428
##
## , , Variedad PEF = No
##
## Meloidogyne sp
## Alternaria passiflorae Si No Sum
## Si 60 22 82
## No 12 30 42
## Sum 72 52 124
##
## , , Variedad PEF = Sum
##
## Meloidogyne sp
## Alternaria passiflorae Si No Sum
## Si 250 140 390
## No 32 130 162
## Sum 282 270 552Como resultado del análisis de la suma de las plantas evaluadas por cada una de las variables, encontramos que la variedad Passiflora edulis flavicarpa ‘PEF’ presentan mayor daño ante presencia de la enfermedad Alternaria passiflorae y la plaga Meloidogyne sp, donde encontramos que de las 428 plantas analizadas 190 presentan esta condición, mientras que con otra variedad Passiflora edulis f. edulis el número de plantas afectadas por las dos condiciones se reduce a 60 plantas, donde el total analizado son 124 plantas. A nivel general observamos que de las 552 plantas objeto de estudio 250 son atacadas por las dos condiciones, mientras que 130 no son afectadas por alguno de los patógenos o enfermedades, mientras que 140 de las plantas son afectadas solo por la enfermedad Alternaria passiflorae y 32 de las plantas son afectadas por Meloidogyne sp.
options(digits = 2)
prop.table(datos,margin=c(1,3))## , , Variedad PEF = Si
##
## Meloidogyne sp
## Alternaria passiflorae Si No
## Si 0.62 0.38
## No 0.17 0.83
##
## , , Variedad PEF = No
##
## Meloidogyne sp
## Alternaria passiflorae Si No
## Si 0.73 0.27
## No 0.29 0.71Del análisis de proporciones observamos que a nivel de la variedad PEF el 62% es afectada por las dos condiciones, mientras que el 38% solo es afectada por la enfermedad Alternaria passiflorae. Por otro lado, cuando no hay presencia de la enfermedad, pero sí de la plaga, se reporta un 17% y sin presencia de alguna enfermedad o patógeno es del 83%.
Para el caso de la variedad PEFE el 73% es afectada por las dos condiciones, mientras que el 27 solo es afectada por la enfermedad Alternaria passiflorae. Por otro lado, cuando no hay presencia de la enfermedad, pero sí de la plaga, se reporta un 26% y sin presencia de alguna enfermedad o patógeno es del 71%.
Cree una tabla que rompa por enfermedad y que se vean las tablas 2x2 para las otras variables.
datos1=array(data=c(170,20,50,15,40,12,30,2),dim=c(2,2,2),dimnames=list('Variedad'=c('Si','No'),'Meloidogyne sp'=c('Si','No'),'Alternaria passiflorae'=c('Si','No')));datos1## , , Alternaria passiflorae = Si
##
## Meloidogyne sp
## Variedad Si No
## Si 170 50
## No 20 15
##
## , , Alternaria passiflorae = No
##
## Meloidogyne sp
## Variedad Si No
## Si 40 30
## No 12 21. ¿Percibe alguna relación?
Con la realización de la nueva tabla observamos una estructura similar a la de los datos obtenidos cuando se realizó la partición por variedad, donde la variedad PEF presenta mayor daño respecto a la variedad PEFE.
2. ¿Aquellos que tienen presencia de plagas son más propensas a la enfermedad?
Con la nueva tabla observamos que 190 son propensas a la enfermedad cuando existe la presencia de estos nematodos, mientras que 65 no tienen presencia de la plaga, pero sí de la enfermedad.
3. ¿Aquellas plantas que son de la variedad PEF son más propensas a la enfermedad?
Así es, hay 220 plantas de esta variedad que tienen dicha enfermedad respecto al total de 255 de plantas que tienen la enfermedad.
III. Modelado Log lineal
Primero se transforman los datos a un data frame, para este análisis como nivel de referencia se fija el No, ya que estamos interesados en saber qué pasa con las plantas que tienen mayor daño en la variedad Passiflora edulis flavicarpa ‘PEF’
datos.df=as.data.frame(as.table(datos))
datos.df[,-4]=lapply(datos.df[,-4],relevel,ref='No')
datos.df## Alternaria.passiflorae Meloidogyne.sp Variedad.PEF Freq
## 1 Si Si Si 190
## 2 No Si Si 20
## 3 Si No Si 118
## 4 No No Si 100
## 5 Si Si No 60
## 6 No Si No 12
## 7 Si No No 22
## 8 No No No 301. Modelo de independencia
Modelo
\[log(\mu_{ijk})=\lambda+\lambda_i+\lambda_j+\lambda_k\\i=1,2\\j=1,2\\k=1,2\]
\(\lambda\) : Conteo
\(\lambda_i\) : Efecto de la enfermedad
\(\lambda_j\) : Efecto de la plaga
\(\lambda_k\) : Efecto de la variedad
Para este modelo se utiliza la función glm (Modelos lineales generalizados)
Análisis de los datos
mod1=glm(Freq~Alternaria.passiflorae+Meloidogyne.sp+Variedad.PEF,data=datos.df,family = poisson)
summary(mod1)##
## Call:
## glm(formula = Freq ~ Alternaria.passiflorae + Meloidogyne.sp +
## Variedad.PEF, family = poisson, data = datos.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 2.76 -6.46 -2.55 4.51 2.16 -1.64 -3.52 2.63
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.8792 0.1197 24.06 <2e-16 ***
## Alternaria.passifloraeSi 0.8786 0.0935 9.40 <2e-16 ***
## Meloidogyne.spSi 0.0435 0.0851 0.51 0.61
## Variedad.PEFSi 1.2388 0.1020 12.15 <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: 377.20 on 7 degrees of freedom
## Residual deviance: 102.77 on 4 degrees of freedom
## AIC: 156.2
##
## Number of Fisher Scoring iterations: 5Con este modelo podemos observar que la enfermedad y la variedad tienen efectos significativos (p-valor<0.05) pero la plaga no tiene efectos significativos (p-valor>0.05).
Cuando la deviancia residual comparada con el número de grados de libertad es muy alto, entonces esto es señal de que el modelo no es tan bueno para este tipo de datos. En este caso hay sobre dispersión porque los datos tienen mayor variabilidad que el modelo Poisson (102.77/4)=25.7 y se espera que la relación sea 1.
Comprobar el modelo
pchisq(deviance(mod1),df=df.residual(mod1),lower.tail = F)## [1] 2.5e-21La hipótesis nula (Ho) es que las frecuencias esperadas satisfacen el modelo log-lineal dado.
Para este caso como el p-valor de 2.5e-21<0.05 entonces se rechaza la hipótesis nula (Ho), así que el modelo no se ajusta.
cbind(mod1$data,fitted(mod1))## Alternaria.passiflorae Meloidogyne.sp Variedad.PEF Freq fitted(mod1)
## 1 Si Si Si 190 154
## 2 No Si Si 20 64
## 3 Si No Si 118 148
## 4 No No Si 100 61
## 5 Si Si No 60 45
## 6 No Si No 12 19
## 7 Si No No 22 43
## 8 No No No 30 18Los valores del modelo están dispersos de los valores observados, por lo tanto, se descarta este modelo.
Ahora se extrae el coeficiente de la plaga para interpretar como probabilidad si lo exponenciamos.
exp(coef(mod1)[3])## Meloidogyne.spSi
## 1prueba_odd=margin.table(datos,margin=2)/sum(margin.table(datos,margin=2))
prueba_odd[1]/prueba_odd[2]## Si
## 1El 1 se conoce como odd o posibilidad, en este caso este odd es el asociado con la presencia de la plaga Meloidogyne.sp. Como observamos con el resultado esta posibilidad es 1, pero a través del cuerpo de la tabla de datos originales observamos que si hay diferencias cuando se presenta esta plaga, respecto a los demás variables.
Grafico
ggplot(data = datos.df, aes(x = Alternaria.passiflorae, y = Freq, color = Meloidogyne.sp))+geom_boxplot()+theme_bw()+ylab('Frecuencia')+xlab('Enfermedad')+facet_wrap(~Meloidogyne.sp)
Con este grafico observamos que la presencia de la enfermedad y de la plaga si influye, dado el alto número de plantas que se encuentran en estas condiciones, mientras que si no hay presencia de la enfermedad se disminuye drásticamente la presencia de la plaga.
2. Modelo de asociación homogénea
Modelo ajustado
\[log(\mu_{ijk})=\lambda+\lambda_i+\lambda_j+\lambda_k+\lambda_{ij}+\lambda_{ik}+\lambda_{jk}\\i=1,2\\j=1,2\\k=1,2\]
\(\lambda\) : Conteo
\(\lambda_i\) : Efecto de la enfermedad
\(\lambda_j\) : Efecto de la plaga
\(\lambda_k\) : Efecto de la variedad
\(\lambda_{ij}\) : interacción entre enfermedad y plaga
\(\lambda_{ik}\) : interacción entre enfermedad y variedad
\(\lambda_{jk}\) : interacción entre plaga y variedad
mod2=glm(Freq~(Alternaria.passiflorae+Meloidogyne.sp+Variedad.PEF)^2,data=datos.df,family=poisson)
summary(mod2)##
## Call:
## glm(formula = Freq ~ (Alternaria.passiflorae + Meloidogyne.sp +
## Variedad.PEF)^2, family = poisson, data = datos.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.0475 -0.1445 -0.0600 0.0655 -0.0840 0.1922 0.1407 -0.1184
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.423 0.169 20.29 < 2e-16
## Alternaria.passifloraeSi -0.362 0.234 -1.55 0.122
## Meloidogyne.spSi -0.994 0.253 -3.93 8.4e-05
## Variedad.PEFSi 1.176 0.189 6.23 4.5e-10
## Alternaria.passifloraeSi:Meloidogyne.spSi 2.038 0.228 8.95 < 2e-16
## Alternaria.passifloraeSi:Variedad.PEFSi 0.539 0.246 2.19 0.028
## Meloidogyne.spSi:Variedad.PEFSi -0.577 0.232 -2.49 0.013
##
## (Intercept) ***
## Alternaria.passifloraeSi
## Meloidogyne.spSi ***
## Variedad.PEFSi ***
## Alternaria.passifloraeSi:Meloidogyne.spSi ***
## Alternaria.passifloraeSi:Variedad.PEFSi *
## Meloidogyne.spSi:Variedad.PEFSi *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 377.19728 on 7 degrees of freedom
## Residual deviance: 0.10885 on 1 degrees of freedom
## AIC: 59.54
##
## Number of Fisher Scoring iterations: 3Con este modelo observamos como primera medida las interacciones presentes, donde se evidencia que todas las tres interacciones son significativas (Alternaria.passiflorae:Meloidogyne.sp, Alternaria.passiflorae:Variedad.PEF y Meloidogyne.sp:Variedad.PEF). Como segunda medida observamos el valor de deviancia de los residuales, donde se tiene que es un valor inferior a uno, lo que indica que es un modelo que se ajusta a los datos obtenidos.
pchisq(deviance(mod2),df=df.residual(mod2),lower.tail = F)## [1] 0.74Para este caso como el p-valor de 0.74>0.05 entonces no se rechaza la hipótesis nula (Ho), así que el modelo se ajusta.
cbind(mod2$data,fitted(mod2))## Alternaria.passiflorae Meloidogyne.sp Variedad.PEF Freq fitted(mod2)
## 1 Si Si Si 190 189
## 2 No Si Si 20 21
## 3 Si No Si 118 119
## 4 No No Si 100 99
## 5 Si Si No 60 61
## 6 No Si No 12 11
## 7 Si No No 22 21
## 8 No No No 30 31Los valores del modelo no difieren demasiado respecto a los valores observados, por lo tanto, no se descarta este modelo.
3. Modelo saturado
mod3=glm(Freq~Alternaria.passiflorae*Variedad.PEF*Meloidogyne.sp,data=datos.df,family=poisson)
summary(mod3)##
## Call:
## glm(formula = Freq ~ Alternaria.passiflorae * Variedad.PEF *
## Meloidogyne.sp, family = poisson, data = datos.df)
##
## Deviance Residuals:
## [1] 0 0 0 0 0 0 0 0
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 3.401 0.183
## Alternaria.passifloraeSi -0.310 0.281
## Variedad.PEFSi 1.204 0.208
## Meloidogyne.spSi -0.916 0.342
## Alternaria.passifloraeSi:Variedad.PEFSi 0.476 0.312
## Alternaria.passifloraeSi:Meloidogyne.spSi 1.920 0.423
## Variedad.PEFSi:Meloidogyne.spSi -0.693 0.420
## Alternaria.passifloraeSi:Variedad.PEFSi:Meloidogyne.spSi 0.166 0.503
## z value Pr(>|z|)
## (Intercept) 18.63 < 2e-16 ***
## Alternaria.passifloraeSi -1.10 0.2692
## Variedad.PEFSi 5.78 7.3e-09 ***
## Meloidogyne.spSi -2.68 0.0073 **
## Alternaria.passifloraeSi:Variedad.PEFSi 1.53 0.1272
## Alternaria.passifloraeSi:Meloidogyne.spSi 4.54 5.6e-06 ***
## Variedad.PEFSi:Meloidogyne.spSi -1.65 0.0991 .
## Alternaria.passifloraeSi:Variedad.PEFSi:Meloidogyne.spSi 0.33 0.7409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3.7720e+02 on 7 degrees of freedom
## Residual deviance: 8.8818e-15 on 0 degrees of freedom
## AIC: 61.44
##
## Number of Fisher Scoring iterations: 3pchisq(deviance(mod3),df=df.residual(mod3),lower.tail = F)## [1] 0cbind(mod3$data,fitted(mod3))## Alternaria.passiflorae Meloidogyne.sp Variedad.PEF Freq fitted(mod3)
## 1 Si Si Si 190 190
## 2 No Si Si 20 20
## 3 Si No Si 118 118
## 4 No No Si 100 100
## 5 Si Si No 60 60
## 6 No Si No 12 12
## 7 Si No No 22 22
## 8 No No No 30 30Este modelo cumple con los requisitos del modelo 2, sin embargo, los valores son exactos y realmente no es bueno, este tipo de modelo se utiliza más como comparación para evaluar otros modelos respecto a este.
IV. Gráficos
plot(mod1,1)
plot(mod2,1)
plot(mod3,1)
A través de las gráficas se observa el comportamiento de los residuales en los 3 modelos, donde en el modelo 1 se aleja de los datos y no coincide con los diferentes puntos, por otro lado, en el modelo 3 se tiene la tendencia logarítmica bien definida, pero al dar tan preciso mejor se toma como referencia, por otro lado, el modelo 2 es el que más se acerca a esa curva del modelo 3.
V. Prueba de verosimilitud
anova(mod2,mod3)## Analysis of Deviance Table
##
## Model 1: Freq ~ (Alternaria.passiflorae + Meloidogyne.sp + Variedad.PEF)^2
## Model 2: Freq ~ Alternaria.passiflorae * Variedad.PEF * Meloidogyne.sp
## Resid. Df Resid. Dev Df Deviance
## 1 1 0.109
## 2 0 0.000 1 0.109pchisq(0.109,df=1,lower.tail = F)## [1] 0.74Con este valor de 0.74, no se rechaza la hipótesis nula de que el modelo 2 (Asociación homogénea) ajusta tan bien como el modelo 3 (saturado) p-valor>0.05. Como en el modelo 2 tenemos que todas las interacciones son significativas, entonces no se hace necesario realizar un ajuste al modelo. Otra forma de comparar cual es el mejor modelo es a través del ACI, como se hace a continuación.
mod2$aic;mod3$aic## [1] 60## [1] 61Con esta comparación nos damos cuenta de que el modelo 2 se ajusta mejor que el modelo 3.
Modelo seleccionado: Modelo 2