Se quiere determinar la relación entre las siguientes variables en las plantas:
set.seed(2020)
datos_mod_log <- c(190,30,116, 95,40, 4, 23, 12)
tabcon_plant <- array( data = datos_mod_log ,
dim = c(2,2,2),
dimnames = list("enferma" = c("SI","NO"),
"clorosis" = c("SI","NO"),
"estres" = c("SI","NO")))
tabcon_plant## , , estres = SI
##
## clorosis
## enferma SI NO
## SI 190 116
## NO 30 95
##
## , , estres = NO
##
## clorosis
## enferma SI NO
## SI 40 23
## NO 4 12Para obtener una mejor tabla donde se relacionen las tres variables, aplanamos la tabla anterior.
ftable(tabcon_plant, row.vars = c("estres", "enferma"))## clorosis SI NO
## estres enferma
## SI SI 190 116
## NO 30 95
## NO SI 40 23
## NO 4 12Calculamos los totales marginales
addmargins(tabcon_plant)## , , estres = SI
##
## clorosis
## enferma SI NO Sum
## SI 190 116 306
## NO 30 95 125
## Sum 220 211 431
##
## , , estres = NO
##
## clorosis
## enferma SI NO Sum
## SI 40 23 63
## NO 4 12 16
## Sum 44 35 79
##
## , , estres = Sum
##
## clorosis
## enferma SI NO Sum
## SI 230 139 369
## NO 34 107 141
## Sum 264 246 510Analizamos las proporciones entre las variables
options(digits = 3)
prop.table(tabcon_plant, margin = c(1,3))## , , estres = SI
##
## clorosis
## enferma SI NO
## SI 0.621 0.379
## NO 0.240 0.760
##
## , , estres = NO
##
## clorosis
## enferma SI NO
## SI 0.635 0.365
## NO 0.250 0.750Se debe convertir la tabla de contingencia en un dataframe para poder ejecutar el modelo
tabcon_plant.df <- as.data.frame(as.table(tabcon_plant))
tabcon_plant.df[,-4] <- lapply(tabcon_plant.df[,-4], relevel, ref = "NO" )
tabcon_plant.df## enferma clorosis estres Freq
## 1 SI SI SI 190
## 2 NO SI SI 30
## 3 SI NO SI 116
## 4 NO NO SI 95
## 5 SI SI NO 40
## 6 NO SI NO 4
## 7 SI NO NO 23
## 8 NO NO NO 12Modelamos la frecuencia como una función de tres variables
mod1 <- glm(Freq ~ enferma + clorosis + estres,
data = tabcon_plant.df, family = poisson)
summary(mod1)##
## Call:
## glm(formula = Freq ~ enferma + clorosis + estres, family = poisson,
## data = tabcon_plant.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 2.187 -4.485 -2.925 4.520 1.816 -2.510 -0.896 0.441
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3547 0.1410 16.69 <2e-16 ***
## enfermaSI 0.9620 0.0990 9.72 <2e-16 ***
## clorosisSI 0.0706 0.0886 0.80 0.43
## estresSI 1.6967 0.1224 13.86 <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: 438.013 on 7 degrees of freedom
## Residual deviance: 64.479 on 4 degrees of freedom
## AIC: 115.9
##
## Number of Fisher Scoring iterations: 5Al ver el residual de la deviancia vemos que el modelo no podría funcionar del todo bien, pues lo idel es que el valor de la deviancia residual se acerque al número de grados de libertad y claramente 64.5 está muy alejado de 4.
Sin embargo, podemos calcular un p-valor para ver si el modelo se ajusta adecuadamente, se debe tener en cuenta que la estadística de la deviancia residual tiene una distribución aproximada de chi-cuadrado.
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 3.31e-13teniendo en cuenta que \(H_0 : las~frecuencias~esperadas~satisfacen~el~modelo~loglineal\) el p-valor rechaza \(H_0\) lo que quiere decir que claramente el modelo no sirve para las frecuencias que se obtuvieron.
Una forma más intuitiva de verificar el ajuste es comparar los valores ajustados con los valores observados.
Hallamos los valores estimados del modelo.
cbind(mod1$data, fitted(mod1))## enferma clorosis estres Freq fitted(mod1)
## 1 SI SI SI 190 161.4
## 2 NO SI SI 30 61.7
## 3 SI NO SI 116 150.4
## 4 NO NO SI 95 57.5
## 5 SI SI NO 40 29.6
## 6 NO SI NO 4 11.3
## 7 SI NO NO 23 27.6
## 8 NO NO NO 12 10.5El modelo ajustado está muy alejado de los datos observados, por ejemplo para las plantas que solo presentaron estrés observamos 95 pero el modelo predice 57. Lo que se confirma que el modelo no se ajusta para nada, por lo cual es ineficiente para tratar los datos que se obtuvieron.
Ajustemos un modelo mas complejo que permita que las variables se asocien entre si, pero que mantenga la misma asociación independiente del nivel de la tercera variable. Ajustar este modelo significa que añadimos interacciones para cada combinación de variables por pares, se podría decir “todas las variables y todas las interacciones por pares”.
mod2 <- glm(Freq ~ (enferma + clorosis + estres)^2,
data = tabcon_plant.df, family = poisson)
summary(mod2)##
## Call:
## glm(formula = Freq ~ (enferma + clorosis + estres)^2, family = poisson,
## data = tabcon_plant.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.00093 0.00234 0.00119 -0.00131 0.00203 -0.00640 -0.00267 0.00370
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4838 0.2613 9.51 < 2e-16 ***
## enfermaSI 0.6522 0.3030 2.15 0.0313 *
## clorosisSI -1.0943 0.3025 -3.62 0.0003 ***
## estresSI 2.0702 0.2733 7.57 3.6e-14 ***
## enfermaSI:clorosisSI 1.6469 0.2247 7.33 2.3e-13 ***
## enfermaSI:estresSI -0.4527 0.3155 -1.43 0.1513
## clorosisSI:estresSI -0.0589 0.2599 -0.23 0.8208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.3801e+02 on 7 degrees of freedom
## Residual deviance: 7.5379e-05 on 1 degrees of freedom
## AIC: 57.43
##
## Number of Fisher Scoring iterations: 3Este modelo se ajusta muchísimo mejor, fijémonos en el valor de la deviancia residual se asemeja al número de grados de libertad, lo que quiere decir que la razón entre ellos se acerca a 1.
Calculamos el p-valor para confirmar el ajuste del modelo
pchisq(deviance(mod2), df = df.residual(mod2), lower.tail = F)## [1] 0.993Obtenemos un p-valor muy alto con lo cual no tenemos suficientes pruebas para rechazar \(H_0\), por lo cual las frecuencias esperadas satisfacen a nuestro modelo.
Nuevamente comparamos los valores ajustados con los observados para ver que tanto coinciden.
cbind(mod2$data, fitted(mod2))## enferma clorosis estres Freq fitted(mod2)
## 1 SI SI SI 190 190.01
## 2 NO SI SI 30 29.99
## 3 SI NO SI 116 115.99
## 4 NO NO SI 95 95.01
## 5 SI SI NO 40 39.99
## 6 NO SI NO 4 4.01
## 7 SI NO NO 23 23.01
## 8 NO NO NO 12 11.99Este modelo se ajusta muy bien, vemos que coinciden casi perfecto los valores observados con los valores estimados por el modelo.
Ya tenemos un modelo que se ajusta perfectamente a los datos que obtuvimos, ahora si vale la pena analizar que efecto tienen las variables entre si, para ello observamos los coeficientes de las interacciones.
Si exponenciamos los coeficientes obtenemos los odds ratios, por ejemplo veamos el coeficiente para “enferma_SI:clorosis_SI”.
exp(coef(mod2)["enfermaSI:clorosisSI"])## enfermaSI:clorosisSI
## 5.19calculamos todos los odd ratios
exp(mod2$coefficients)## (Intercept) enfermaSI clorosisSI
## 11.987 1.920 0.335
## estresSI enfermaSI:clorosisSI enfermaSI:estresSI
## 7.926 5.191 0.636
## clorosisSI:estresSI
## 0.943También es buena idea calcular un intervalo de confianza para la estimación de los odds ratios.
exp(confint(mod2, parm = c("enfermaSI:clorosisSI","enfermaSI:estresSI","clorosisSI:estresSI")))## Waiting for profiling to be done...## 2.5 % 97.5 %
## enfermaSI:clorosisSI 3.373 8.15
## enfermaSI:estresSI 0.334 1.16
## clorosisSI:estresSI 0.564 1.57Las asociaciones son casi todas muy fuertes. Vemos que las probabilidades de que presente una planta clorosis si está enferma son 3 y hasta 8 veces más altas que las probabilidades de que presente clorosis si no se encuentra enferma, es decir, que la clorosis es una señal de que la planta se encuentra enferma.
Por su parte el estrés hídrico tiene influencia aunque no tan fuerte en la presencia de la enfermedad ya que por cada 3 plantas con estrés, se presenta un con enfermedad, entonces se puede decir que el estrés hídrico aumenta la probabilidad de que sea infectada, siendo el estrés una puerta para los patógenos para causar alguna enfermedad.
En el modelo anterior se vio que la única interacción necesaria es la de “enferma_SI:clorosis_SI” por lo que ahora no tendremos en cuenta las otras interacciones.
mod3 <- glm(Freq ~ (enferma + clorosis)^2,
data = tabcon_plant.df, family = poisson)
summary(mod3)##
## Call:
## glm(formula = Freq ~ (enferma + clorosis)^2, family = poisson,
## data = tabcon_plant.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 6.39 2.84 5.08 5.11 -8.09 -3.80 -6.49 -6.87
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.9797 0.0967 41.17 < 2e-16 ***
## enfermaSI 0.2616 0.1286 2.03 0.042 *
## clorosisSI -1.1465 0.1969 -5.82 5.8e-09 ***
## enfermaSI:clorosisSI 1.6501 0.2243 7.36 1.9e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 438.01 on 7 degrees of freedom
## Residual deviance: 270.01 on 4 degrees of freedom
## AIC: 321.4
##
## Number of Fisher Scoring iterations: 5Curiosamente este modelo se desajusta bastante a los que se había logrado en el modelo anterior donde se tenía en cuenta todas las interacciones, el valor de la deviancia residual está muy alejado del número de grados de libertad. Además el valor del AIC es mucho más elevado tanto que es peor que el primer modelo ejecutado, por lo que no vale la pena seguir valorando su ajuste, este modelo no sirve para las frecuencias observadas.
Es necesario evaluar cómo se comparta un modelo que tenga en cuenta TODAS las interacciones entre las variables, la manera de expresa que se ajusten todas las interacciones es mediante la fórmula “estrés * enferma * clorosis”
mod4 <- glm(Freq ~ enferma*clorosis*estres,
data = tabcon_plant.df, family = poisson)
summary(mod4)##
## Call:
## glm(formula = Freq ~ enferma * clorosis * estres, family = poisson,
## data = tabcon_plant.df)
##
## Deviance Residuals:
## [1] 0 0 0 0 0 0 0 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.48491 0.28868 8.61 < 2e-16 ***
## enfermaSI 0.65059 0.35611 1.83 0.0677 .
## clorosisSI -1.09861 0.57735 -1.90 0.0571 .
## estresSI 2.06897 0.30637 6.75 1.4e-11 ***
## enfermaSI:clorosisSI 1.65200 0.63389 2.61 0.0092 **
## enfermaSI:estresSI -0.45087 0.38205 -1.18 0.2379
## clorosisSI:estresSI -0.05407 0.61416 -0.09 0.9298
## enfermaSI:clorosisSI:estresSI -0.00588 0.67790 -0.01 0.9931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.3801e+02 on 7 degrees of freedom
## Residual deviance: 8.4377e-15 on 0 degrees of freedom
## AIC: 59.43
##
## Number of Fisher Scoring iterations: 3No se evidencia que la interacción triple sea necesaria, sin embargo vamos a comprar los valores estimados con los observados.
cbind(mod4$data, fitted(mod4))## enferma clorosis estres Freq fitted(mod4)
## 1 SI SI SI 190 190
## 2 NO SI SI 30 30
## 3 SI NO SI 116 116
## 4 NO NO SI 95 95
## 5 SI SI NO 40 40
## 6 NO SI NO 4 4
## 7 SI NO NO 23 23
## 8 NO NO NO 12 12Sorprende como los valores estimados coinciden perfectamente con los valores observados, esto es característico del modelo saturado.
En igualdad de condiciones, se prefiere un modelo más simple. Normalmente no se quiere terminar con un modelo saturado que se ajuste perfectamente a nuestros datos, sin embargo es útil ajustar un modelo saturado para verificar que la interacción de orden superior no es estadísticamente significativa.
Para verificar que el modelo de asociación homogénea se ajusta tan bien como el modelo saturado se realiza una prueba de la razón de verosimilidad.
anova(mod2,mod4)## Analysis of Deviance Table
##
## Model 1: Freq ~ (enferma + clorosis + estres)^2
## Model 2: Freq ~ enferma * clorosis * estres
## Resid. Df Resid. Dev Df Deviance
## 1 1 7.54e-05
## 2 0 0.00e+00 1 7.54e-05pchisq(7.54e-05, df = 1, lower.tail = F)## [1] 0.993Teniendo en cuenta que $ H_0 : elmodelomássencilloeselmejor $ con ese p-valor no se rechaza \(H_0\) lo que quiere decir que aunque el modelo saturado coincida perfectamente con los datos observados es mejor el modelo de asociación homogénea.
library(plotly)## Warning: package 'plotly' was built under R version 4.0.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.0.3##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutfig_mods <- plot_ly(data=tabcon_plant.df,y=~tabcon_plant.df$Freq,type='bar',name = "Observado")
fig_mods <- fig_mods %>% add_trace(y=~mod4$fitted.values,name = "Estimado Modelo 4")
fig_mods <- fig_mods %>% add_trace(y=~mod2$fitted.values,name = "Estimado Modelo 2")
fig_mods <- fig_mods %>% add_trace(y=~mod1$fitted.values,name = "Estimado Modelo 1")
fig_mods <- fig_mods %>% add_trace(y=~mod3$fitted.values,name = "Estimado Modelo 3")
fig_mods <- fig_mods %>% layout(xaxis=list(title="Categoria"),
yaxis=list(title="Frecuencias"))
fig_mods## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Se debe tener en cuenta que no se trata de seleccionar modelos que solo hagan predicciones, sino seleccionar modelos que aunque la predicción no sea perfecta pero si muy cercana, las variables que queden dentro del modelo evidencien relaciones estadísticas con la variable respuesta.
Por lo que se concluye que el mejor modelo es el número 2, que corresponde al modelo de asociación homogénea ya que este modelo se ajusta muy bien y aunque no es perfecto como el modelo saturado (modelo 4) se acerca a ello, además se determinó que la interacción tripla es innecesaria por lo que el modelo 2 es el más simple que se ajusta de la mejor manera.
moraleja: “no solo son predicciones perfectas, son modelos que se acerquen a los valores observados pero que las variables que lo componen sean útiles”