Jesús David Ortiz

Marly Daniela Quiroga

Erika Milena Aparicio

Cristian Danilo Romero

Ana María Montaño Hernández

En la siguiente tabla se registran los resultados obtenidos en un cultivo de maíz con respecto a la presencia de tres enfermedades: Manchas foliares o tizón, pudrición de la raíz y pudrición del tallo.

Se tienen tres variables, cada una con 2 niveles (presencia o ausencia).

Tabla de contingencia

Uno de los objetivos es ver si hay propensión a la pudrición del tallo como consecuencia de la presencia de manchas foliare o la pudrición de la raíz. Observar la probabilidad de una variable (la presencia de una enfermedad) con respecto a las otras.

Se van a probar independencias, mirando la relación entre las tres variables, sin que necesariamente una sea variable respuesta.

tabcon <- array(data = c(180, 10, 108, 90, 50, 2, 12, 20),
                dim = c(2,2,2),
                dimnames = list("raiz" = c("S", "N"), 
                                "tallo" = c("S", "N"),
                                "manchas"= c("S", "N")))   # Se cargan los datos 
tabcon
## , , manchas = S
## 
##     tallo
## raiz   S   N
##    S 180 108
##    N  10  90
## 
## , , manchas = N
## 
##     tallo
## raiz  S  N
##    S 50 12
##    N  2 20
ftable(tabcon, row.vars = c("manchas", "raiz"))
##              tallo   S   N
## manchas raiz              
## S       S          180 108
##         N           10  90
## N       S           50  12
##         N            2  20
addmargins(tabcon)
## , , manchas = S
## 
##      tallo
## raiz    S   N Sum
##   S   180 108 288
##   N    10  90 100
##   Sum 190 198 388
## 
## , , manchas = N
## 
##      tallo
## raiz   S  N Sum
##   S   50 12  62
##   N    2 20  22
##   Sum 52 32  84
## 
## , , manchas = Sum
## 
##      tallo
## raiz    S   N Sum
##   S   230 120 350
##   N    12 110 122
##   Sum 242 230 472

Se observa que hubo un total de 472 plantas analizadas.

Además, se observa que la relación entre las plantas que presentaron manchas con respecto a las que no es de 388:84, aproximadamente 5:1. Por cada 5 plantas que presentan manchas 1 no tiene.

options(digits = 3)
prop.table(tabcon, margin = c(1,3))   # se convierten los datos anteriores en proporciones
## , , manchas = S
## 
##     tallo
## raiz     S     N
##    S 0.625 0.375
##    N 0.100 0.900
## 
## , , manchas = N
## 
##     tallo
## raiz      S     N
##    S 0.8065 0.194
##    N 0.0909 0.909

De las que presentan manchas, el porcentaje más alto está en las que no presentan pudrición ni de raíz ni de tallo. Ocurre lo mismo en las plantas sin manchas foliares.

Pregunta: Cree una tabla que rompa por presencia de pudrición en la raíz y que se vean las tablas 2*2 para las otras dos variables. ¿ precibe alguna relación? ¿Aquellas plantas que presentan pudrición en la raíz son más propensos a la presencia de pudrición en el tallo? ¿Aquellas plantas que tienen manchas foliares son más propensos a la pudrición de la raíz?

tabcon2 <- array(data = c(180, 50, 108, 12, 10, 2, 90, 20),
                dim = c(2,2,2),
                dimnames = list("manchas" = c("S", "N"), 
                                "tallo" = c("S", "N"),
                                "raiz"= c("S", "N")))
addmargins(tabcon2)
## , , raiz = S
## 
##        tallo
## manchas   S   N Sum
##     S   180 108 288
##     N    50  12  62
##     Sum 230 120 350
## 
## , , raiz = N
## 
##        tallo
## manchas  S   N Sum
##     S   10  90 100
##     N    2  20  22
##     Sum 12 110 122
## 
## , , raiz = Sum
## 
##        tallo
## manchas   S   N Sum
##     S   190 198 388
##     N    52  32  84
##     Sum 242 230 472

En la tabla anterior se observa que la relación entre las plantas que presentaron pudrición de la raíz con respecto a las que no es de 350:122, aproximadamente 3:1. Por tanto, por cada 3 plantas que presentan pudrición en la raíz, 1 no se ve afectada por esta enfermedad.

options(digits = 3)
prop.table(tabcon2, margin = c(1,3))
## , , raiz = S
## 
##        tallo
## manchas     S     N
##       S 0.625 0.375
##       N 0.806 0.194
## 
## , , raiz = N
## 
##        tallo
## manchas      S     N
##       S 0.1000 0.900
##       N 0.0909 0.909

Asimismo, con las proporciones se puede observar, que de las plantas que presentan pudrición en la raíz, el porcentaje más alto se encuentra en las que a su vez presentaron pudrición en el tallo y manchas foliares con un 80.6%. Por tanto, las plantas con pudrición en la raíz son más propensas a presentar las otras dos enfermedades. Además, el 19% (el menor porcentaje) de las plantas con pudrición de la raíz no se vieron afectadas por las otras dos enfermedades.

Si se analiza sólo la propensión a la pudrición del tallo en el caso que presente pudrición de la raíz, sin importar la presencia o ausencia de manchas foliares, se observa que presentan los dos más altos porcentajes; 62.5% y 80.6%, demostrando que la pudrición de la raíz hace más propensa a la planta a la presencia, que a la ausencia de pudrición del tallo.

Por otra parte, si se desea analizar la propensión de las plantas a la pudrición de la raíz, cuando hay manchas foliares, se observa, que de las plantas con manchas, el 62.5% y el 37.5% (con respecto a la pudrición del tallo) presentaron pudrición de la raíz, obteniendo valores más altos que aquellas sin pudrición de la raíz, siendo propensas a esta enfermedad. Sin embargo, es importante analizar las tras variables al tiempo ya que los porcentajes también van a depender de la presencia o ausencia de pudrición en el tallo.

# Es una tabla de contingencia 2^3

tabcon.df <- as.data.frame(as.table(tabcon))
tabcon.df[,-4] <- lapply(tabcon.df[,-4], relevel, ref = "N")  # Se agrega una columna con los conteos, frecuencias. 

tabcon.df
##   raiz tallo manchas 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   20

Modelo loglineal de INDEPENDENCIA para tablas de contingencia de dos vías

Modelo sin interacción. Tienen que ser independientes las tres variables.

\[ log(\mu_{ij})=log(n \pi_{i}\pi_{j})= log(n) + log\pi_{i}+log\pi_{j}\]

Para este caso, tenemos una tabla de contingencia de tres vías. Modelo de independencia.

mod0 <- glm(Freq ~ raiz + tallo + manchas, data = tabcon.df, family = poisson)
summary(mod0) # Se llama análisis de devianza 
## 
## Call:
## glm(formula = Freq ~ raiz + tallo + manchas, family = poisson, 
##     data = tabcon.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 ***
## raizS         1.0539     0.1051   10.02   <2e-16 ***
## talloS        0.0509     0.0921    0.55     0.58    
## manchasS      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: 5

Según el análisis de devianza, parece que el efecto de la pudrición del tallo no es significativo.

Cuando la devianza residual (133.63) comparada con los grados de libertad (4), presenta valores muy alejados, es señal de que el modelo está mal y no modela adecuadamente estos datos. Se tiene un problema de sobredispersión, es decir, los datos tienen más variabilidad que la que se espera en un modelo Poisson.

pchisq(deviance(mod0), df= df.residual(mod0), lower.tail = F) # se obtiene un p_valor cercano a 0, señal de que el modelo es MALO. No sirve. 
## [1] 6.53e-28
cbind(mod0$data, fitted(mod0))
##   raiz tallo manchas 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.6

En la tabla anterior se combinan los datos originales con los valores ajustados por el modelo. Se obtiene que, el modelo ajustado está muy lejos de los datos observados, corroborando un problema de sobredispersión muy grande. Por ejemplo, se observaron 180 plantas con la presencia de las tres enfermedades, y el modelo solamente predice 147.

exp(coef(mod0)[3])  # Coeficiente de la pudrición de tallos. Aplica anti-logaritmo
## talloS 
##   1.05

El 1.05 se conoce como odd o posibilidad. Es el odd de la presencia de pudrición de tallos, lo que significa que la posibilidad de la presencia o ausencia de esta enfermedad es casi igual, independientemente de las otras dos enfermedades. Sin embargo, puede que no sea una buena estimación si se compara con los datos de la tabla original de contingencia. No ajusta correctamente.

prueba_odd = margin.table(tabcon, margin = 2)/sum(margin.table(tabcon, margin = 2))
prueba_odd[1]/prueba_odd[2]
##    S 
## 1.05

Ahora, se va a ajustar un modelo más complejo que permita que las variables se asocien entre sí pero que mantengan la misma asociación independientemente del nivel de la tercera variable. Se deja una variable fija y se estudia la independencia de las otras dos. Esto se conoce como asociación homogénea:

mod1aju <- glm(Freq ~ (raiz + tallo + manchas)^2, data = tabcon.df, family = poisson)
summary(mod1aju)
## 
## Call:
## glm(formula = Freq ~ (raiz + tallo + manchas)^2, family = poisson, 
##     data = tabcon.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 ***
## raizS             -0.307      0.315   -0.97    0.330    
## talloS            -1.645      0.373   -4.41  1.1e-05 ***
## manchasS           1.602      0.240    6.67  2.6e-11 ***
## raizS:talloS       2.918      0.329    8.88  < 2e-16 ***
## raizS:manchasS     0.458      0.332    1.38    0.167    
## talloS:manchasS   -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: 4

Se observa que en este caso, la devianza residual (1.5132) si tiene un valor cercano a los grados de libertad (1). El modelo se ajusta mejor a los datos. No hay problema de sobredispersión y se puede observar que el p_valor es más cercano a 1.

Además, se observa un AIC menor que en el modelo anterior, corroborando que este modelo es más preciso.

pchisq(deviance(mod1aju), df= df.residual(mod1aju), lower.tail = F)
## [1] 0.219
cbind(mod1aju$data, fitted(mod1aju))
##   raiz tallo manchas Freq fitted(mod1aju)
## 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.44

En este caso, el modelo se ajusta más con los datos. Se observan valores cercanos entre la frecuencia y el valor que predice el valor.

Coeficientes:

exp(coef(mod1aju)["talloS:manchasS"])
## talloS:manchasS 
##           0.477

Este odd de 0.477 indica aproximadamente una relación de 1:2, es decir, por cada planta que presenta pudrición de tallos dos presentan manchas foliares.

mod1aju$coefficients
##     (Intercept)           raizS          talloS        manchasS    raizS:talloS 
##           2.915          -0.307          -1.645           1.602           2.918 
##  raizS:manchasS talloS:manchasS 
##           0.458          -0.739

Ahora, con la función confint se va a calcular un intervalo de confianza para estas estimaciones de odds ratios (posibilidades)

exp(confint(mod1aju, parm = c("raizS:talloS", "raizS:manchasS", "talloS:manchasS")))
## Waiting for profiling to be done...
##                  2.5 % 97.5 %
## raizS:talloS    10.092 36.965
## raizS:manchasS   0.822  3.037
## talloS:manchasS  0.262  0.844

Un odd de 10.092 indica que por cada 10 plantas que tienen pudrición en la raíz, 1 presenta pudrición en el tallo. Sin embargo, este odd puede alcanzar un valor de 36.965. Se evidencia algún tipo de relación entre la pudrición de raíces y la de tallos.

Asimismo, los odds en los otros dos casos (interacciones), son más bajos que el de raíz:tallo. Un odd de 0.822 representa una relación de 4:5, es decir por cada 4 plantas con pudrición en la raíz, 5 presentan manchas foliares.

Otra pregunta:

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. Haga algunos gráficos definidos por ustedes para visualizar el ajuste, y para comparar las posibilidades de presentar pudrición en el tallo si se tiene pudrición de la raíz independientemente de la presencia de manchas. Escriba la ecuación final del modelo ajustado.

mod2aju <- glm(Freq ~ (raiz + tallo)^2 + (manchas + tallo)^2 , data = tabcon.df, family = poisson)
summary(mod2aju)
## 
## Call:
## glm(formula = Freq ~ (raiz + tallo)^2 + (manchas + tallo)^2, 
##     family = poisson, data = tabcon.df)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.0431   0.1866   0.4586  -0.4866   0.0821  -0.3752  -1.2106   1.1456  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.728      0.190   14.38  < 2e-16 ***
## raizS              0.087      0.132    0.66    0.510    
## talloS            -1.781      0.367   -4.86  1.2e-06 ***
## manchasS           1.823      0.190    9.57  < 2e-16 ***
## raizS:talloS       2.866      0.324    8.84  < 2e-16 ***
## talloS:manchasS   -0.527      0.247   -2.14    0.033 *  
## ---
## 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:   3.4093  on 2  degrees of freedom
## AIC: 57
## 
## Number of Fisher Scoring iterations: 4

Al eliminar la interacción doble raízS:manchasS que se observó en el modelo pasado que no era significativa, se oberva que aunque el AIC disminuyó un poco, es casi el mismo al modelo anterior, con un valor de 57. Asimismo, la devianza residual (3.4093) también presentó un valor cercano a los grados de libertad (2). Al presentar ambos modelos un AIC igual, podría decirse que el modelo anterior se ajusta más a los datos originales ya que su devianza residual (1.5132) tiene un valor más cercano a los grados de libertad (1), que cuando se elimina la interacción. Además, el modelo anterior presenta un p_valor mayor, menos cercano a 0, indicando que va a existir menos problema de sobredispersión.

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

Ya no se obtiene un pvalor cercano a 0. Es bueno el modelo.

cbind(mod2aju$data, fitted(mod2aju))
##   raiz tallo manchas Freq fitted(mod2aju)
## 1    S     S       S  180          180.58
## 2    N     S       S   10            9.42
## 3    S     N       S  108          103.30
## 4    N     N       S   90           94.70
## 5    S     S       N   50           49.42
## 6    N     S       N    2            2.58
## 7    S     N       N   12           16.70
## 8    N     N       N   20           15.30

Ahora, si comparamos los valores originales con los valores ajustados, de ambos modelos, se observa que ambos modelos son bastante buenos y se acercan a los datos originales.

mod2aju$coefficients
##     (Intercept)           raizS          talloS        manchasS    raizS:talloS 
##           2.728           0.087          -1.781           1.823           2.866 
## talloS:manchasS 
##          -0.527
exp(confint(mod2aju, parm = c("raizS:talloS", "talloS:manchasS")))
## Waiting for profiling to be done...
##                 2.5 % 97.5 %
## raizS:talloS    9.662 34.790
## talloS:manchasS 0.361  0.953

Ahora, al aplicarle intervalos de confianza al modelo ajustado eliminando la interacción se obtiene:

Primero, un odd de 9.662 indica que por cada 10 plantas aproximadamente que presentan pudrición de la raíz, 1 presenta pudrición en el tallo. Sin embargo, este odd puede alcanzar un valor de 34.790. Se evidencia algún tipo de relación entre la pudrición de raíces y la de tallos.

El otro odd de 0.361 indica que por cada planta con pudrición en los tallos, 3 presentan manchas foliares (una relación 1:3) y este odd puede alcanzar un valor de 0.953.

library(ggplot2)
ggplot(tabcon.df, aes(x = raiz, y = Freq, fill = tallo)) + 
  stat_identity(geom = "bar", position = "dodge") +
  facet_wrap(~ tallo) +
  labs(xlab ="Pudrición de la raíz", title = "Relación entre la pudrición del tallo y la raíz", fill = "Pudrición del tallo") +
  theme_bw()

En la gráfica se observa una mayor propensión a presentar pudrición del tallo cuando se presenta pudrición de la raíz, ya que es la mayor frecuencia de plantas las que presentaron estas dos enfermedades. Lo anterior, independientemente de si presentan manchas foliares o no.

plot(mod2aju)

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

## Warning in sqrt(crit * p * (1 - hh)/hh): Se han producido NaNs

Modelo de asociación homogénea: El siguiente modelo permite todas las asociaciones condicionales de las tres variables (las tres enfermedades) por parejas, lo que corresponde a una asociación homogénea (XY, XZ, YZ)

\[ log(\mu_{ijk})= \lambda + \lambda_{i}^X + \lambda_{j}^Y + \lambda_{k}^Z +\lambda_{ij}^{XY} +\lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \] donde,

X: pudrición de las raíces

Y: pudrición del tallo

Z: presencia de manchas

Al no existir la interacción raízS:manchasS, se elimina el término \(\lambda_{ik}^{XZ}\), y el modelo ajustado queda:

\[ log(\mu_{ijk})= \lambda + \lambda_{i}^X + \lambda_{j}^Y + \lambda_{k}^Z +\lambda_{ij}^{XY} + \lambda_{jk}^{YZ} \]