Modelo log lineal

Nombres:

• Kahomy Lorena Avila Robayo - 1000854408

• Carlos Manuel Burgos De La Cruz - 1001089617

El siguiente ejercicio es didáctico y no tiene en cuenta las caracteristicas, ni el comportamiento de cada enfermedad.

Se desea estudiar la correlacion de estas tres enfermedades en una variedad de papa

1. tizon tardio

2. marchitez bacteriana

3. pierna negra

Se quiere saber si hay propensión a presentar la enfermedas piel negra, como consecuencia de tener tizón tardio y/o marchitez bacteriana

• Se inicia generando 8 datos, con una media de 120 y una varianza de 40

set.seed(2020)
datos = (round(rnorm(8, 120, 40)))
datos

## [1] 135 132  76  75   8 149 158 111

tabcon = array(data=datos, 
                dim = c(2,2,2), 
                dimnames = list("tis"=c("si","no"),
                                "mar"=c("si","no"), 
                                "pier"=c("si","no")))
tabcon

## , , pier = si
## 
##     mar
## tis   si no
##   si 135 76
##   no 132 75
## 
## , , pier = no
## 
##     mar
## tis   si  no
##   si   8 158
##   no 149 111

• Usamos la función ftable para aplanar la tabla

ftable(tabcon, row.vars=c("tis","mar"))

##         pier  si  no
## tis mar             
## si  si       135   8
##     no        76 158
## no  si       132 149
##     no        75 111

• Con addmargis calculamos totales marginales de la tabla

addmargins(tabcon)

## , , pier = si
## 
##      mar
## tis    si  no Sum
##   si  135  76 211
##   no  132  75 207
##   Sum 267 151 418
## 
## , , pier = no
## 
##      mar
## tis    si  no Sum
##   si    8 158 166
##   no  149 111 260
##   Sum 157 269 426
## 
## , , pier = Sum
## 
##      mar
## tis    si  no Sum
##   si  143 234 377
##   no  281 186 467
##   Sum 424 420 844

• Proporciones

options(digit=3)
prop.table(tabcon,margin = c(1,3))

## , , pier = si
## 
##     mar
## tis         si        no
##   si 0.6398104 0.3601896
##   no 0.6376812 0.3623188
## 
## , , pier = no
## 
##     mar
## tis          si        no
##   si 0.04819277 0.9518072
##   no 0.57307692 0.4269231

Mediante las tablas anteriores, podemos observar que al tener la enfermedad de piel negra, hay una 64% de probabilidad de tener tanto marchitez bacteriana, como tizón tardio, mientras que al no presentar piel negra hay un 48% de probabilidad.

Modelo log lineal

Se genera el dataframe a partir de la tabla de contingencia, con tres columnas asociadas a las enfermedades y sus respectivas frecuencias.

tabcon.df <- as.data.frame(as.table(tabcon))
tabcon.df[,-4] <- lapply(tabcon.df[,-4], relevel, ref= "no")
tabcon.df

##   tis mar pier Freq
## 1  si  si   si  135
## 2  no  si   si  132
## 3  si  no   si   76
## 4  no  no   si   75
## 5  si  si   no    8
## 6  no  si   no  149
## 7  si  no   no  158
## 8  no  no   no  111

Modelo glm

Se usará la función glm de los modelos lineales generalizado.

Modelo log lineal de independencia para tabla de contingencia de dos vías.

\[ log(\mu_{ij}) = log~n + log~ \pi_i + log~\pi_j \]

Se procede con el modelado:

mod0 <- glm(Freq ~ mar + pier + tis,
            data = tabcon.df,family = poisson)
summary(mod0)

## 
## Call:
## glm(formula = Freq ~ mar + pier + tis, family = poisson, data = tabcon.df)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
##   3.9888    1.4351   -1.8125   -3.9937  -11.6404    2.7010    5.9300   -0.5868  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.764718   0.067082  71.028  < 2e-16 ***
## marsi        0.009479   0.068844   0.138  0.89049    
## piersi      -0.018958   0.068846  -0.275  0.78303    
## tissi       -0.214084   0.069238  -3.092  0.00199 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 225.22  on 7  degrees of freedom
## Residual deviance: 215.51  on 4  degrees of freedom
## AIC: 273.53
## 
## Number of Fisher Scoring iterations: 5

El modelo (mod0), no ajusta para estos datos con un modelo de Poisson, ya que la deviancia resiudal de 215.5 con 4 grados de libertad, da como resultado un número muy grande, mostrando así que los datos tienen más variabilidad de lo esperado para un modelo de Poisson.

• Prueba pchisq

\[H_o = Las~frecuencias~esperadas~satisfacen~el~modelo~loglineal\]

pchisq(deviance(mod0), df=df.residual(mod0), lower.tail = F)

## [1] 1.736506e-45

Con la prueba pchisq comprobamos al rechazarse la \(H_0\), que el modelo no es confianble y no es util al momento de concluir si existe una relación entre las variables

cbind(mod0$data, fitted(mod0))

##   tis mar pier Freq fitted(mod0)
## 1  si  si   si  135     93.79908
## 2  no  si   si  132    116.19144
## 3  si  no   si   76     92.91419
## 4  no  no   si   75    115.09529
## 5  si  si   no    8     95.59428
## 6  no  si   no  149    118.41520
## 7  si  no   no  158     94.69245
## 8  no  no   no  111    117.29807

En esta tabla podemos ver que los datos observados se alejan mucho de los datos esperados, evidenciando que le modelo no ajusta

exp(coef(mod0)[3])

##    piersi 
## 0.9812207

• calculo del odd, coeficiente de dispariedad

prueba_odd=margin.table(tabcon,margin = 2)/sum(margin.table(tabcon,margin = 2))

prueba_odd[1]/prueba_odd[2]

##       si 
## 1.009524

Modelo de asociación homogénea con interacciones dobles

A continuación correremos el modelo de asociación homogénea con interacciones dobles, buscando un modelo que se ajuste mejor a los datos observados.

\[ log\mu_{ikj} = \lambda+\lambda^X_i+\lambda^Y_j+\lambda^Z_k+\lambda^{XY}_{ij}+\lambda^{YZ}_{jk}+\lambda^{XZ}_{ik} \]

mod1 <- glm(Freq ~ (mar + pier + tis)^2,
            data=tabcon.df, family = poisson)
summary(mod1)

## 
## Call:
## glm(formula = Freq ~ (mar + pier + tis)^2, family = poisson, 
##     data = tabcon.df)
## 
## Deviance Residuals: 
##      1       2       3       4       5       6       7       8  
##  2.547  -2.243  -2.836   3.672  -5.578   2.406   2.327  -2.419  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.93068    0.08177  60.296  < 2e-16 ***
## marsi        -0.13056    0.11403  -1.145    0.252    
## piersi       -1.06941    0.14169  -7.547 4.44e-14 ***
## tissi        -0.05911    0.11231  -0.526    0.599    
## marsi:piersi  1.34117    0.15499   8.653  < 2e-16 ***
## marsi:tissi  -1.17220    0.15598  -7.515 5.68e-14 ***
## piersi:tissi  0.83717    0.15585   5.372 7.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 225.217  on 7  degrees of freedom
## Residual deviance:  81.223  on 1  degrees of freedom
## AIC: 145.25
## 
## Number of Fisher Scoring iterations: 5

El modelo (mod1), no ajusta para estos datos con un modelo de asociación homogenea, ya que la deviancia resiudal de 81.22 con 1 grado de libertad, da como resultado un número grande, aunque más pequeño que el mod0, mostrando así que los datos tienen más variabilidad de lo esperad.

\[H_o = Las~frecuencias~esperadas~satisfacen~el~modelo~loglineal\]

pchisq(deviance(mod1), df=df.residual(mod1), lower.tail = F)

## [1] 2.016544e-19

Con la prueba pchisq comprobamos al rechazarse la \(H_0\), que el modelo no es confianble y no es util al momento de concluir si existe una relación entre las variables

cbind(mod1$data, fitted(mod1))

##   tis mar pier Freq fitted(mod1)
## 1  si  si   si  135    107.52573
## 2  no  si   si  132    159.47427
## 3  si  no   si   76    103.47427
## 4  no  no   si   75     47.52573
## 5  si  si   no    8     35.47427
## 6  no  si   no  149    121.52573
## 7  si  no   no  158    130.52573
## 8  no  no   no  111    138.47427

En esta tabla podemos ver que los datos observados se alejan de los datos esperados, evidenciando que le modelo no ajusta

exp(coef(mod1)["piersi:tissi"])

## piersi:tissi 
##     2.309811

mod1$coefficients

##  (Intercept)        marsi       piersi        tissi marsi:piersi  marsi:tissi 
##   4.93068454  -0.13055853  -1.06941330  -0.05911417   1.34116989  -1.17220418 
## piersi:tissi 
##   0.83716592

exp(confint(mod1,parm = c("marsi:piersi", "marsi:tissi", "piersi:tissi")))

## Waiting for profiling to be done...

##                  2.5 %    97.5 %
## marsi:piersi 2.8315121 5.2009817
## marsi:tissi  0.2272457 0.4190405
## piersi:tissi 1.7069365 3.1460574

• Las probabilidades que presente piel negra, si tiene marchitez son aprox 3 veces más altas

• Las probabilidades de que presenten tizón tardio, si tiene marchitez, son de aprox 0

• La probabilidad de que presente tizón tardio, si tiene piel negra con de aprox 2 vesces más altas

Gráfico mod1

# install.packages("plotly")
library(plotly)

## Loading required package: ggplot2

## 
## 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':
## 
##     layout

fig_exp <- plot_ly(data=tabcon.df,y=~tabcon.df$Freq,type='bar',name = "Observado")
fig_exp <- fig_exp %>% add_trace(y=~mod1$fitted.values,name = "Esperado")
fig_exp <- fig_exp %>% layout(xaxis=list(title="Enfermedad"),
                              yaxis=list(title="Frecuencias"))

fig_exp

## 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.

En este gráfico evidenciamos lo que anteriormente se menciona, acerca de las fercuencias esperadas y observadas.

Modelo saturado (todas las interacciones, simples dobles y triples)

Seguimos buscando un modelo que se ajuste mejor a las frecuencias observadas

\[ log\mu_{ikj} = \lambda+\lambda^X_i+\lambda^Y_j+\lambda^Z_k+\lambda^{XY}_{ij}+\lambda^{YZ}_{jk}+\lambda^{XZ}_{ik} +\lambda^{XYZ}_{ijk} \]

mod2 <- glm(Freq ~ mar*pier*tis, data = tabcon.df, family = poisson)
summary(mod2)

## 
## Call:
## glm(formula = Freq ~ mar * pier * tis, family = poisson, data = tabcon.df)
## 
## Deviance Residuals: 
## [1]  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         4.70953    0.09492  49.618  < 2e-16 ***
## marsi               0.29442    0.12538   2.348  0.01887 *  
## piersi             -0.39204    0.14947  -2.623  0.00872 ** 
## tissi               0.35306    0.12385   2.851  0.00436 ** 
## marsi:piersi        0.27090    0.19139   1.415  0.15694    
## marsi:tissi        -3.27757    0.38347  -8.547  < 2e-16 ***
## piersi:tissi       -0.33982    0.20452  -1.662  0.09661 .  
## marsi:piersi:tissi  3.28680    0.43419   7.570 3.74e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.2522e+02  on 7  degrees of freedom
## Residual deviance: 7.5495e-15  on 0  degrees of freedom
## AIC: 66.026
## 
## Number of Fisher Scoring iterations: 3

Observando lo deviancia residual de 7.5495e-15, con o grados de libertad, podemos concluir que este es modelo el que mejor se ajusta a la frecuecia de los datos obtenidos

cbind(mod1$data, fitted(mod2))

##   tis mar pier Freq fitted(mod2)
## 1  si  si   si  135          135
## 2  no  si   si  132          132
## 3  si  no   si   76           76
## 4  no  no   si   75           75
## 5  si  si   no    8            8
## 6  no  si   no  149          149
## 7  si  no   no  158          158
## 8  no  no   no  111          111

Aquí observamos que los datos esperados y observados concuerdan.

\[H_o = Las~frecuencias~esperadas~satisfacen~el~modelo~loglineal\]

pchisq(deviance(mod2), df=df.residual(mod2), lower.tail = F)

## [1] 0

El valor de p = 0, no rechaza la \(H_0\), por lo que este modelo ajusta de mejor manera

Grafico mod2

#install.packages("plotly")
library(plotly)
fig_exp <- plot_ly(data=tabcon.df,y=~tabcon.df$Freq,type='bar',name = "Observado")
fig_exp <- fig_exp %>% add_trace(y=~mod2$fitted.values,name = "Esperado")
fig_exp <- fig_exp %>% layout(xaxis=list(title="Enfermedad"),
                              yaxis=list(title="Frecuencias"))

fig_exp

Con este gráfico podemos apoyar la conclusión anteriormente planteada acerca del ajuste del modelo y la concordancia entre los valores observados y esperados.