Ejercicio planteado.

La tabla 1 resume las resultados de un estudio realizado en Pereira, dónde se observó la presencia o ausencia de malezas, patógenos y enfermedades en los culivos de café.

knitr::include_graphics("/cloud/project/tabla2.PNG")

Tabla 1 Tabla con los resultados del estudio.

Interesa comprender la relación entre las tres variables sin que una sea necesariamente la “respuesta”, de este modo se usa un modelo loglineal.

data <- array (data = c (911, 44, 538, 456, 3, 2, 43, 279), 
                         dim = c (2,2,2), 
                         dimnames = list ("Malezas" = c ("Sí", "No"),
                                         "Patogenos" = c ("Sí", "No"),
                                         "Enfermedades" = c ("Sí", "No")))

# array: Función para crear una tabla con más de dos dimensiones. En esta ocasión el argumento dim dice que queremos crear una tabla con 2 filas, 2 columnas y 2 capas.

# dimnames, es un argumento que proporciona nombres para las dimensiones. Requiere un objeto de lista, por lo que  la entrada del argumento en la función list.

data
## , , Enfermedades = Sí
## 
##        Patogenos
## Malezas  Sí  No
##      Sí 911 538
##      No  44 456
## 
## , , Enfermedades = No
## 
##        Patogenos
## Malezas Sí  No
##      Sí  3  43
##      No  2 279

Al imprimir se obtienen dos tablas de 2x2.

ftable (data, row.vars = c ("Malezas", "Patogenos")) 
##                   Enfermedades  Sí  No
## Malezas Patogenos                     
## Sí      Sí                     911   3
##         No                     538  43
## No      Sí                      44   2
##         No                     456 279
# Para aplanar la tabla se usa ftable, el argumento row.vars especifica qué variables queremos en las filas y su respectivo orden.
addmargins (data) 
## , , Enfermedades = Sí
## 
##        Patogenos
## Malezas  Sí  No  Sum
##     Sí  911 538 1449
##     No   44 456  500
##     Sum 955 994 1949
## 
## , , Enfermedades = No
## 
##        Patogenos
## Malezas Sí  No Sum
##     Sí   3  43  46
##     No   2 279 281
##     Sum  5 322 327
## 
## , , Enfermedades = Sum
## 
##        Patogenos
## Malezas  Sí   No  Sum
##     Sí  914  581 1495
##     No   46  735  781
##     Sum 960 1316 2276
# Antes de continuar con el modelado, se pueden explorar los datos. La función addmargins proporciona totales marginales.

Se observa que se realizarón 2276 observaciones en el estudio de la región y la mayoría de los resultados muestran presencia de malezas en los cultivos de café. La relación de presencia de enfermedades respecto a la ausencia es de 1949:327.

 prop.table(data, margin = c(1,3))
## , , Enfermedades = Sí
## 
##        Patogenos
## Malezas        Sí        No
##      Sí 0.6287095 0.3712905
##      No 0.0880000 0.9120000
## 
## , , Enfermedades = No
## 
##        Patogenos
## Malezas          Sí        No
##      Sí 0.065217391 0.9347826
##      No 0.007117438 0.9928826
# La función prop.table permite calcular proporciones de celda

En el estudio aquellos cultivos con presencia de malezas y patogenos, el 62% también tenia presencia de enfermedades. Asimismo, de aquellos que no tenian malezas ni patogenos, el 99% tampoco tenia enfermedades. Definitivamente parece haber algún tipo de relación aquí.

 prop.table(data, margin = c(1,3))
## , , Enfermedades = Sí
## 
##        Patogenos
## Malezas        Sí        No
##      Sí 0.6287095 0.3712905
##      No 0.0880000 0.9120000
## 
## , , Enfermedades = No
## 
##        Patogenos
## Malezas          Sí        No
##      Sí 0.065217391 0.9347826
##      No 0.007117438 0.9928826
# La función prop.table permite calcular proporciones de celda
Marco de datos
data.df <- as.data.frame(as.table(data))
data.df[,-4] <- lapply(data.df[,-4], relevel, ref = "No")
data.df
##   Malezas Patogenos Enfermedades Freq
## 1      Sí        Sí           Sí  911
## 2      No        Sí           Sí   44
## 3      Sí        No           Sí  538
## 4      No        No           Sí  456
## 5      Sí        Sí           No    3
## 6      No        Sí           No    2
## 7      Sí        No           No   43
## 8      No        No           No  279
# as.data.frame vuelve un objeto de tabla en R a un marco de datos con una columna para frecuencias de celda, donde cada fila representa una combinación única de variables. 

# as.table permite convertir la  matriz en una tabla.

# Los coeficientes del modelo se expresarán en términos de respuestas “sí”.
Modelo de independencia - Se asume que las tres variables son independienters entre sí

Supone que las tres variables son independientes entre sí, lo que parece muy poco probable.

mod0 <- glm(Freq ~ Patogenos + Enfermedades + Malezas, 
            data = data.df, family = poisson)
summary(mod0)
## 
## Call:
## glm(formula = Freq ~ Patogenos + Enfermedades + Malezas, family = poisson, 
##     data = data.df)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  14.522  -17.683   -7.817    3.426  -12.440   -8.832   -8.436   19.639  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     4.17254    0.06496  64.234  < 2e-16 ***
## PatogenosSí    -0.31542    0.04244  -7.431 1.08e-13 ***
## EnfermedadesSí  1.78511    0.05976  29.872  < 2e-16 ***
## MalezasSí       0.64931    0.04415  14.707  < 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: 2851.5  on 7  degrees of freedom
## Residual deviance: 1286.0  on 4  degrees of freedom
## AIC: 1343.1
## 
## Number of Fisher Scoring iterations: 6
# La función glm se usa para trabajar con  modelos loglineales. Se modela con  Freq como una función de las tres variables usando glm. 

# En family se establece el argumento en poisson, ya que estamos modelando recuentos.

Mirando el resumen, parece que este es un gran modelo. Vemos coeficientes altamente significativos y valores p cercanos a 0. Pero para los modelos loglineales se quiere verificar la desviación residual. Como regla general, nos gustaría que tuviera un valor cercano a los grados de libertad. Aquí tenemos 1286 en 4 grados de libertad.Esto indica un mal ajuste.

Calculo del p.valor
pchisq(deviance(mod0), df = df.residual (mod0), lower.tail = F)
## [1] 3.574437e-277

Este p.valor nos indica que el modelo no es bueno, por ende se procede a realizarlo como un modelo de valores estimados.

Valores estimados del modelo

Una forma más intuitiva de investigar el ajuste es comparar los valores ajustados con los valores observados.

cbind(mod0$data, fitted(mod0))
##   Malezas Patogenos Enfermedades Freq fitted(mod0)
## 1      Sí        Sí           Sí  911    539.98258
## 2      No        Sí           Sí   44    282.09123
## 3      Sí        No           Sí  538    740.22612
## 4      No        No           Sí  456    386.70007
## 5      Sí        Sí           No    3     90.59739
## 6      No        Sí           No    2     47.32880
## 7      Sí        No           No   43    124.19392
## 8      No        No           No  279     64.87990

El modelo ajustado esta muy lejos de los valores observados. Por ejemplo, observamos que 3 resultados del estudio que tenian presencia de malezas y patogenos pero no de enfermedadades, y el modelo predice alrededor de 90 cultivos bajo estas condiciones.

margin.table(data, margin = 2)/sum(margin.table(data, margin = 2))
## Patogenos
##        Sí        No 
## 0.4217926 0.5782074
0.4217926 / 0.5782074
## [1] 0.7294832

Según este modelo, las probabilidades de tener enfermedades son de aproximadamente 0,73 a 1, independientemente de si hay presencia de malezas o patogenos Sin embargo, sabemos que probablemente no sea una buena estimación.

Modelo de asociación homogeneo - Modelo que trabajan con dos variables, fijando una tercera

Ajustemos un modelo más complejo que permite asociar variables entre sí, pero mantiene la misma asociación independientemente del nivel de la tercera variable. Esto quiere decir que por ejemplo, la presencia de patogenos y enfermedades tienen algún tipo de relación, pero esa relación es la misma independientemente de si hay presencia o no de malezas.

mod1 <- glm(Freq ~ (Patogenos + Enfermedades + Malezas)^2, 
            data = data.df, family = poisson)
summary(mod1)
## 
## Call:
## glm(formula = Freq ~ (Patogenos + Enfermedades + Malezas)^2, 
##     family = poisson, data = data.df)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
##  0.02044  -0.09256  -0.02658   0.02890  -0.33428   0.49134   0.09452  -0.03690  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 5.63342    0.05970  94.361  < 2e-16 ***
## PatogenosSí                -5.30904    0.47520 -11.172  < 2e-16 ***
## EnfermedadesSí              0.48772    0.07577   6.437 1.22e-10 ***
## MalezasSí                  -1.88667    0.16270 -11.596  < 2e-16 ***
## PatogenosSí:EnfermedadesSí  2.98601    0.46468   6.426 1.31e-10 ***
## PatogenosSí:MalezasSí       2.84789    0.16384  17.382  < 2e-16 ***
## EnfermedadesSí:MalezasSí    2.05453    0.17406  11.803  < 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: 2851.46098  on 7  degrees of freedom
## Residual deviance:    0.37399  on 1  degrees of freedom
## AIC: 63.417
## 
## Number of Fisher Scoring iterations: 4

Este modelo encaja mucho mejor. Observe la desviación residual (0.37399) en comparación con los grados de libertad (1).

pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)
## [1] 0.5408396

El alto valor p dice que no tenemos pruebas suficientes para rechazar la hipótesis nula de que las frecuencias esperadas satisfacen nuestro modelo.

Una vez más, podemos comparar los valores ajustados y observados y ver qué tan bien coinciden:

cbind(mod1$data, fitted(mod1))
##   Malezas Patogenos Enfermedades Freq fitted(mod1)
## 1      Sí        Sí           Sí  911    910.38317
## 2      No        Sí           Sí   44     44.61683
## 3      Sí        No           Sí  538    538.61683
## 4      No        No           Sí  456    455.38317
## 5      Sí        Sí           No    3      3.61683
## 6      No        Sí           No    2      1.38317
## 7      Sí        No           No   43     42.38317
## 8      No        No           No  279    279.61683
Coeficientes de las interacciones
exp(coef(mod1)["EnfermedadesSí:MalezasSí"])
## EnfermedadesSí:MalezasSí 
##                 7.803201

Los cultivos con presencia de enfermedades tienen probabilidades estimadas de tener presencia de malezas que son 7 veces mayores que las probabilidades estimadas de los cultivos que no tienen enfermedades. Debido a que ajustamos un modelo de asociación homogéneo, la razón de posibilidades es la misma independientemente de si hay presnecia de patógenos o no.

Es una buena idea calcular un intervalo de confianza para estas estimaciones de razón de posibilidades. Podemos hacer eso con la función confint:

exp(confint(mod1, parm = c("PatogenosSí:EnfermedadesSí",
                                "PatogenosSí:MalezasSí",
                                "EnfermedadesSí:MalezasSí")))
## Waiting for profiling to be done...
##                                2.5 %   97.5 %
## PatogenosSí:EnfermedadesSí  8.814046 56.64360
## PatogenosSí:MalezasSí      12.645760 24.06925
## EnfermedadesSí:MalezasSí    5.601452 11.09715

Preguntas

1. Cree una tabla que rompa por presencia de patógenos y que se vean las tablas 2*2 para las otras dos variables. ¿precibe alguna relación? ¿Aquellos cultivos con presencia de enfermedades son más propensos a tener patógenos? ¿Aquellos cultivos con presencia de malezas son más propensos a tener patogenos?

data1 <- array (data = c (911, 44, 538, 456, 3, 2, 43, 279), 
                         dim = c (2,2,2), 
                         dimnames = list ("Malezas" = c ("Sí", "No"),
                                         "Enfermedades" = c ("Sí", "No"),
                                         "Patógenos" = c ("Sí", "No")))
data1
## , , Patógenos = Sí
## 
##        Enfermedades
## Malezas  Sí  No
##      Sí 911 538
##      No  44 456
## 
## , , Patógenos = No
## 
##        Enfermedades
## Malezas Sí  No
##      Sí  3  43
##      No  2 279
 prop.table(data1, margin = c(1,3))
## , , Patógenos = Sí
## 
##        Enfermedades
## Malezas        Sí        No
##      Sí 0.6287095 0.3712905
##      No 0.0880000 0.9120000
## 
## , , Patógenos = No
## 
##        Enfermedades
## Malezas          Sí        No
##      Sí 0.065217391 0.9347826
##      No 0.007117438 0.9928826

De la tabla se puede evidenciar que aquellos cultivos de café que tengan presencia de una enfermedad son mas propensos a tener malezas, donde el 62% de estos tendrían a su vez patógenos. Adémas se puede observar que un 6% de los cultivos que no tienen patógenos, cuentan con presencia de enfermedad y malezas, frente a un 99% que no tienen enfermedad ni malezas, y los patógenos estan ausentes. De estos análisis, podría decirse (por ahora) que se evidencia algún tipo de relación entre estas tres variables.

2. Elimine interacciones dobles innecesarias y repita el proceso. ¿Mejora el modelo? Estime todos los valores y use el AIC para seleccionar el mejor modelo. Interprete las relaciones encontradas.

En el modelo todas las interacciones dobles son significativas, debido a que el p.valor es inferior al 5%, por ende no se realizó este punto.

3.* ¿Qué sucede si ajustamos un modelo con una interacción de tres vías?

El modelo con interacción de tres vías permite que las asociaciones por pares cambien según el nivel de la tercera variable, en un modelo saturado.

Modelo saturado.
mod2 <- glm(Freq ~ (Patogenos * Enfermedades * Malezas)^2, 
            data = data.df, family = poisson)
summary(mod2)
## 
## Call:
## glm(formula = Freq ~ (Patogenos * Enfermedades * Malezas)^2, 
##     family = poisson, data = data.df)
## 
## Deviance Residuals: 
## [1]  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                           5.63121    0.05987  94.060  < 2e-16 ***
## PatogenosSí                          -4.93806    0.70964  -6.959 3.44e-12 ***
## EnfermedadesSí                        0.49128    0.07601   6.464 1.02e-10 ***
## MalezasSí                            -1.87001    0.16383 -11.414  < 2e-16 ***
## PatogenosSí:EnfermedadesSí            2.59976    0.72698   3.576 0.000349 ***
## PatogenosSí:MalezasSí                 2.27548    0.92746   2.453 0.014149 *  
## EnfermedadesSí:MalezasSí              2.03538    0.17576  11.580  < 2e-16 ***
## PatogenosSí:EnfermedadesSí:MalezasSí  0.58951    0.94236   0.626 0.531600    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.8515e+03  on 7  degrees of freedom
## Residual deviance: 1.5388e-13  on 0  degrees of freedom
## AIC: 65.043
## 
## Number of Fisher Scoring iterations: 3

La desviación de este modelo es básicamente 0 en 0 grados de libertad. Los recuentos ajustados coinciden con los recuentos observados, como se observa a continuación:

deviance(mod2)
## [1] 1.538769e-13
cbind(mod2$data, fitted(mod2))
##   Malezas Patogenos Enfermedades Freq fitted(mod2)
## 1      Sí        Sí           Sí  911          911
## 2      No        Sí           Sí   44           44
## 3      Sí        No           Sí  538          538
## 4      No        No           Sí  456          456
## 5      Sí        Sí           No    3            3
## 6      No        Sí           No    2            2
## 7      Sí        No           No   43           43
## 8      No        No           No  279          279

En igualdad de condiciones, preferimos un modelo más simple. Por lo general, no queremos terminar con un modelo saturado que se ajuste perfectamente a nuestros datos. Sin embargo, podemos verificar que el modelo de asociación homogéneo se ajusta tan bien como el modelo saturado realizando una prueba de razón de verosimilitud:

anova (mod1, mod2)
## Analysis of Deviance Table
## 
## Model 1: Freq ~ (Patogenos + Enfermedades + Malezas)^2
## Model 2: Freq ~ (Patogenos * Enfermedades * Malezas)^2
##   Resid. Df Resid. Dev Df Deviance
## 1         1    0.37399            
## 2         0    0.00000  1  0.37399
pchisq(0.373, df = 1, lower.tail = F)
## [1] 0.5413735

Esto dice que no rechazamos la hipótesis nula de que mod1 se ajusta tan bien como mod2.

REFERENCIAS

-Universidad de Virginia. 2016. Introducción a los modelos loglineales. Disponible en: https://data.library.virginia.edu/an-introduction-to-loglinear-models/