Modelos Log-linear
Cristian David Muñoz Mayorga
Generacion de Datos
tabcon <- array(data = c(180,10, 108,90,50, 2, 12, 20),
dim = c(2,2,2),
dimnames = list("Cig" = c("S","N"),
"Mar" = c("S","N"),
"Alc" = c("S","N")))Tabla de contingencia
ftable(tabcon, row.vars = c("Alc","Cig"))## Mar S N
## Alc Cig
## S S 180 108
## N 10 90
## N S 50 12
## N 2 20Analisis de los datos:
Margenes
addmargins(tabcon)## , , Alc = S
##
## Mar
## Cig S N Sum
## S 180 108 288
## N 10 90 100
## Sum 190 198 388
##
## , , Alc = N
##
## Mar
## Cig S N Sum
## S 50 12 62
## N 2 20 22
## Sum 52 32 84
##
## , , Alc = Sum
##
## Mar
## Cig S N Sum
## S 230 120 350
## N 12 110 122
## Sum 242 230 472Tabla de proporciones
options(digits=3)
prop.table(tabcon, margin = c(1,3))## , , Alc = S
##
## Mar
## Cig S N
## S 0.625 0.375
## N 0.100 0.900
##
## , , Alc = N
##
## Mar
## Cig S N
## S 0.8065 0.194
## N 0.0909 0.909Analisis grafico
Diagrama de Sieve
library(vcd)## Warning: package 'vcd' was built under R version 4.0.3## Loading required package: gridvcd::mosaic( ~ Mar + Cig + Alc, data = tabcon)
sieve(tabcon, shade=TRUE)
Odds ratios para Cig y Mar por Alc:
oddsratio(tabcon, log=FALSE)## odds ratios for Cig and Mar by Alc
##
## S N
## 15.0 41.7odds_ratio <- oddsratio(tabcon) # capture log odds ratios
plot(odds_ratio) 
Preparacipon de los datos en dataframe para Modelado:
x = as.table(tabcon)
tabcon.df <- as.data.frame(x)
tabcon.df[,-4] <- lapply(tabcon.df[,-4], relevel, ref = "N")
tabcon.df## Cig Mar Alc 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 20Modelo Log Lineal de independencia:
Este de modelo de independencia completa establece que las filas, columnas y capas son todas independientes las unas de las otras en la distribución de Poissosn, y por lo tanto no hay interacciones.
\[\text{log}(\mu_{ijk})=\pi+\pi_i^{Cig}+\pi_j^{Alc}+\pi_k^{Mar}\]
log_lineal_ind <- glm(Freq ~ Cig + Mar + Alc,
data = tabcon.df, family = poisson)
summary(log_lineal_ind)##
## Call:
## glm(formula = Freq ~ Cig + Mar + Alc, 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 ***
## CigS 1.0539 0.1051 10.02 <2e-16 ***
## MarS 0.0509 0.0921 0.55 0.58
## AlcS 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: 5plot(log_lineal_ind,1)
El modelo log lineal no es bueno pues el AIC es muy alto y la deviancia residual de 133.63 y los grados de libertad de 4 difieren bastante entre si.
Mediante el grafico se puede observar que el modelo no se ajsuta a los valores esperados.
Probabilidad de encontrar deviancia igual o mayor a la del modelo de independencia:
pchisq(deviance(log_lineal_ind), df = df.residual(log_lineal_ind), lower.tail = F)## [1] 6.53e-28Se extraen los grados de libertad y se observa que la probabilidad calculada es muy baja, por lo que se rechaza la hipotesis de encontrar una deviancia igual o mayor a la del modelo lo cual indica que las frecuencias esperadas no satisfacen a las del modelo
Coeficiente para "Marihuana Sí"" del modelo de independencia:
coef(log_lineal_ind)[3]## MarS
## 0.0509Exponenciacion o Antilogaritimo del coeficiente para hallar el Odds:
exp(coef(log_lineal_ind)[3])## MarS
## 1.05Según este modelo, las posibilidades (odd) de probar marihuana son de aproximadamente 1,05 a 1, independientemente de que hayas probado el alcohol o los cigarrillos. Pero se sabe que el modelo no es acertado segun el diagrama, la tabla y la prueba pchisq, y ademas se puede confirmar que las tres variables no son independientes entre si.
Modelo Saturado de interaccion triple:
Los modelos Log-linear son jerárquicos, por lo que siempre es bueno comenzar con los modelos mas complejos y ver que tanto se puede reducir.
\[\text{log}(\mu_{ijk})=\pi+\pi_i^{Cig}+\pi_j^{Alc}+\pi_k^{Mar}+\pi_{ij}^{CigAlc}+\pi_{jk}^{AlcMar}+\pi_{ik}^{CigMar}+\pi_{ijk}^{CigAlcMar}\]
int_triple <- glm(Freq ~ Cig * Mar * Alc,
data = tabcon.df, family = poisson)
summary(int_triple)##
## Call:
## glm(formula = Freq ~ Cig * Mar * Alc, 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) 2.996 0.224 13.40 < 2e-16 ***
## CigS -0.511 0.365 -1.40 0.1618
## MarS -2.303 0.742 -3.10 0.0019 **
## AlcS 1.504 0.247 6.08 1.2e-09 ***
## CigS:MarS 3.730 0.808 4.61 3.9e-06 ***
## CigS:AlcS 0.693 0.392 1.77 0.0771 .
## MarS:AlcS 0.105 0.813 0.13 0.8969
## CigS:MarS:AlcS -1.022 0.883 -1.16 0.2471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.6107e+02 on 7 degrees of freedom
## Residual deviance: 1.2879e-14 on 0 degrees of freedom
## AIC: 57.6
##
## Number of Fisher Scoring iterations: 3plot(int_triple,1)
Este modelo parece ser mejor al de independencia por su AIC menor y por poseer unos residuales similares a los grados de libertad. Aunque no se puede asegurar que es el mejor pues solo una de sus interacciones es significativa lo cual no tendria sentido. No es un modelo preferible, pero con el se verifica que la interaccion de orden superior no es estadisticamente significativa.
Mediante el grafico se puede observar que el modelo se ajusta bastante a los valores esperados.
Comparación de frecuencias
cbind(tabcon, fitted(int_triple))## tabcon
## 1 180 180
## 2 10 10
## 3 108 108
## 4 90 90
## 5 50 50
## 6 2 2
## 7 12 12
## 8 20 20La comparacion entre las frecuencias esperadas y las ajustadas son exactas.
Modelo de asociación homogéneo:
\[\text{log}(\mu_{ijk})=\pi+\pi_i^{Cig}+\pi_j^{Alc}+\pi_k^{Mar}+\pi_{ij}^{CigAlc}+\pi_{jk}^{AlcMar}+\pi_{ik}^{CigMar}\]
aso_hom <- glm(Freq ~ (Cig + Mar + Alc)^2,
data = tabcon.df, family = poisson)
summary(aso_hom)##
## Call:
## glm(formula = Freq ~ (Cig + Mar + Alc)^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 ***
## CigS -0.307 0.315 -0.97 0.330
## MarS -1.645 0.373 -4.41 1.1e-05 ***
## AlcS 1.602 0.240 6.67 2.6e-11 ***
## CigS:MarS 2.918 0.329 8.88 < 2e-16 ***
## CigS:AlcS 0.458 0.332 1.38 0.167
## MarS:AlcS -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: 4plot(aso_hom,1)
Este modelo tiene un AIC bajo al igual que el de interacción triple, lo cual es bueno, y, ademas, se encontró una interaccion mas que el modelo de interacción triple que indica que sería un mejor modelo. Por otro, lado la desviancia residual y los grados de libertad son silimares, aunque no tanto como el modelo saturado.
Mediante el grafico se puede observar que el modelo no se ajusta perfectamente a los valores esperados pero lo hace de manera muy cercana.
cbind(aso_hom$data, fitted(aso_hom))## Cig Mar Alc Freq fitted(aso_hom)
## 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.44Probabilidad de encontrar deviancia igual o mayor a la del modelo de independencia:
pchisq(deviance(aso_hom), df = df.residual(aso_hom), lower.tail = F)## [1] 0.219El p-valor obtenido indica que no hay suficiente evidencia que puede rechazar la hipotesis de que las frecuencias esperadas no satisfagan el modelo.
Comparacion entre Modelos de Asociacíon homogenea y interacción triple:
anova(aso_hom,int_triple)## Analysis of Deviance Table
##
## Model 1: Freq ~ (Cig + Mar + Alc)^2
## Model 2: Freq ~ Cig * Mar * Alc
## Resid. Df Resid. Dev Df Deviance
## 1 1 1.51
## 2 0 0.00 1 1.51El analisis de desviancia arroja un resultado de 1.51 que se usa en la prueba pchisq.
pchisq(1.51, df = 1, lower.tail = F)## [1] 0.219Con un p-valor de 0.22 no se rechaza la hipotesis de que el modelo de asociación homogeneo ajusta tan bien como el modelo saturado de interacción de tres vias. Esto quiere decir que a pesar de que el modelo de asociación homogeneo no ajusta sus valores tan bien como el de saturado, es capaz de arrojar mejores resultados por ser un modelo mas simple y que las variables que estan en el modelo evidencien relacion estadistica con la respuesta por medio de la significancia en las interacciones.
Extracción de coeficientes y aplicacion de antilogaritmo:
exp(coef(aso_hom)["MarS:AlcS"])## MarS:AlcS
## 0.477Este Odd se puede interpretar como que de los estudiantes que probaron alcohol por cada estudiante que probo marihuana dos estudiantes no probaron marihuana.
exp(confint(aso_hom, parm = c("CigS:MarS","CigS:AlcS","MarS:AlcS")))## Waiting for profiling to be done...## 2.5 % 97.5 %
## CigS:MarS 10.092 36.965
## CigS:AlcS 0.822 3.037
## MarS:AlcS 0.262 0.844Un ODDS de 10 en Cig:Mar significa que son 10 veces mas grandes la posibilidaes de si probar marihuana si se probo cigarillos que de haber probado marihuana no habiendo probado cigarillos.
Modelo de Independencia condicional:
A diferencia del modelo de independencia completa, este modelo asume que hay independecia para una variable, mientras que en las demas se dan interacciones por pares. En este caso, las asociaciones entre cigarrillo y marihuana no depende de los niveles de alcohol, y las asociaciones entre alcohol y marihuana no dependen de los niveles del cigarrillo, por lo que la interaccion cigarrillo alcohol queda fuera del modelo al ser consideradas independientes la una de la otra.
\[\text{log}(\mu_{ijk})=\pi+\pi_i^{Cig}+\pi_j^{Alc}+\pi_k^{Mar}+\pi_{ij}^{CigAlc}+\pi_{jk}^{AlcMar}\]
indep_cond = glm(Freq ~ (Cig * Mar) + (Alc * Mar),
data = tabcon.df, family = poisson)
summary(indep_cond)##
## Call:
## glm(formula = Freq ~ (Cig * Mar) + (Alc * Mar), 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 ***
## CigS 0.087 0.132 0.66 0.510
## MarS -1.781 0.367 -4.86 1.2e-06 ***
## AlcS 1.823 0.190 9.57 < 2e-16 ***
## CigS:MarS 2.866 0.324 8.84 < 2e-16 ***
## MarS:AlcS -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: 4Este modelo arroja un AIC bajo aceptable pero la deviancia residual es casi el doble a la de los grados de libertad. Tambien se encuentran p-valores significativos para las interacciones.
plot(indep_cond, 1)
Mediante el grafico se puede observar que el modelo ajusta un poco menos que el de asociación homogenea en los valores esperados.
Comparacion entre Modelos de Asociacíon homogenea e independencia condicionada:
anova(indep_cond,aso_hom)## Analysis of Deviance Table
##
## Model 1: Freq ~ (Cig * Mar) + (Alc * Mar)
## Model 2: Freq ~ (Cig + Mar + Alc)^2
## Resid. Df Resid. Dev Df Deviance
## 1 2 3.41
## 2 1 1.51 1 1.9pchisq(1.9, df = 1, lower.tail = F)## [1] 0.168Con un p-valor de 0.16 no se rechaza la hipotesis de que el modelo de independencia condicionada ajusta tan bien como el modelo asociación homogeneo.
exp(confint(indep_cond, parm = c("CigS:MarS","MarS:AlcS")))## Waiting for profiling to be done...## 2.5 % 97.5 %
## CigS:MarS 9.662 34.790
## MarS:AlcS 0.361 0.953La estimación de los odds ratio son muy parecidos respecto al del modelo de asociación homogeneo.
Si se quisiera cortar por Cigarillo:
por_cig <- array(data = c(180,50, 108,12,10, 2, 90, 20), dim = c(2,2,2),
dimnames = list("Alc" = c("S","N"),
"Mar" = c("S","N"),
"Cig" = c("S","N"))); por_cig## , , Cig = S
##
## Mar
## Alc S N
## S 180 108
## N 50 12
##
## , , Cig = N
##
## Mar
## Alc S N
## S 10 90
## N 2 20ftable(por_cig, row.vars = c("Cig","Alc"))## Mar S N
## Cig Alc
## S S 180 108
## N 50 12
## N S 10 90
## N 2 20addmargins(por_cig)## , , Cig = S
##
## Mar
## Alc S N Sum
## S 180 108 288
## N 50 12 62
## Sum 230 120 350
##
## , , Cig = N
##
## Mar
## Alc S N Sum
## S 10 90 100
## N 2 20 22
## Sum 12 110 122
##
## , , Cig = Sum
##
## Mar
## Alc S N Sum
## S 190 198 388
## N 52 32 84
## Sum 242 230 472options(digits=3)
prop.table(por_cig, margin = c(1,3))## , , Cig = S
##
## Mar
## Alc S N
## S 0.625 0.375
## N 0.806 0.194
##
## , , Cig = N
##
## Mar
## Alc S N
## S 0.1000 0.900
## N 0.0909 0.909- ¿Aquellos que fuman cigarros son más propensos al uso de la marihuana?
De acuerdo a las proporciones obtenidas, aquellos que consumen cigarillo son mas propensos al consumo de marihuana puesto que en aquellos que consumen alcohol el 62.5% consumen la marihuana y en el grupo de los que no consumen alcohol el 81% si consume marihuana. Por otro lado, los que no consumen cigarrillo tiene un comportamiento contrario: el 90% aproximadamente no son consumidores de marihuana idependiente de si consume o no alcohol.
- ¿Aquellos que toman alcohol son más propensos al uso de la cigarrillos?
Por medio de los margenes se puede observar que de las 350 personas que consumen cigarillo, 288 personas consumen alcohol, mientras que 62 no lo consumen, mientras que en el grupo que no consume cigarrillo 100 consumen alcohol y 22 no, por lo que en ambos grupos hay una proporcion muy similar indicando que no existe una relacion entre el consumo de alcohol y el cigarrillo, haciendolos, efectivamente, ni mas ni menos propensos.
Selección del modelo:
Con el modelo de asociación homogeneo no se encontró una interacción significativa en cigarrillos:alcohol y por esta razon se decide probar con el modelo de independencia condicionada. El modelo de independencia condicionada es un buen modelo, tan bueno como el de asociación homogeneo y al realizarce una analisis descriptivo se concluye que no es neceasario dejar esta interacción en el modelo, pues no es significativa la interacción. De este modo, el modelo de independencia condicionada, en este caso, se ajusta mejor que el de asociación homogeneo pues es mas simple y obtiene buenas estimaciones en los odds ratio.