Los datos y códigos explorados en esta ocasión son los de Agresti (1996, Tabla 6.3), código de la Univeridad de Virginia y apoyo en las cifras de clase (diapositivas).
| Uso de alcohol | Uso de cigarrillo | Uso de Marihuana | ||
|---|---|---|---|---|
| Si | No | |||
| Si | Si | 180 | 108 | |
| No | 10 | 90 | ||
| No | Si | 50 | 12 | |
| No | 2 | 20 |
La tabla presentada se basa en una encuesta a estudiantes de último año de secundaria donde se quizo evaluar el consumo de alcohol, cigarrillos y/o marihuana.
Al basarse en un conteo con ciertos niveles de variables categóricas (cualitativas) se propone un modelo log-lineal para esta tabla de múltiples entradas. Se busca evaluar la interacción o relación entre las variables sin tomar una en especifico como variable respuesta.
tabcon <- array(data = c(180,10, 108,90,50, 2, 12, 20),
dim = c(2,2,2),
dimnames = list("Cigarrillo" = c("Si","No"),
"Marihuana" = c("Si","No"),
"Alcohol" = c("Si","No"))) ftable(tabcon, row.vars = c("Alcohol","Cigarrillo"))## Marihuana Si No
## Alcohol Cigarrillo
## Si Si 180 108
## No 10 90
## No Si 50 12
## No 2 20#-La función “ftable” permite organizar en una tabla las 2 creadas anteriormente y row.vars especifica las variables y orden en las filas.Antes de modelar, se desea explorar los datos. La función \(addmargins\) proporciona totales marginales. Vemos que 472 estudiantes participaron en la encuesta y la mayoría de ellos probaron el alcohol.
addmargins(tabcon)## , , Alcohol = Si
##
## Marihuana
## Cigarrillo Si No Sum
## Si 180 108 288
## No 10 90 100
## Sum 190 198 388
##
## , , Alcohol = No
##
## Marihuana
## Cigarrillo Si No Sum
## Si 50 12 62
## No 2 20 22
## Sum 52 32 84
##
## , , Alcohol = Sum
##
## Marihuana
## Cigarrillo Si No Sum
## Si 230 120 350
## No 12 110 122
## Sum 242 230 472options(digits=3)
prop.table(tabcon, margin = c(1,3))## , , Alcohol = Si
##
## Marihuana
## Cigarrillo Si No
## Si 0.625 0.375
## No 0.100 0.900
##
## , , Alcohol = No
##
## Marihuana
## Cigarrillo Si No
## Si 0.8065 0.194
## No 0.0909 0.909Vemos que de los estudiantes encuestados que probaron los cigarrillos y el alcohol, el 62.5% también probó la marihuana. Por otro lado, de aquellos estudiantes que no probaron alcohol, pero si cigarrillo tienen un 80.65% de probabilidad de usar marihuana, mientras que el 90% que no usa alcohol ni cigarrillo tampoco probó la marihuana. Esto sugiere relación más amplia entre consumir cigarrillo y marihuana vs alcohol y marihuana.
Con esta aproximación consideramos usar la función \(glm\) para trabajar mediante transformación de los datos. También es posible considerar \(loglin\), que es especifica para tablas de contingencia, pero glm permite usar fórmulas.
Se transforma en un marco de datos. Usamos \(“as.data.frame”\) para convertir el objeto tabcon en una tabla con las frecuencias de la combinación de variables (2 niveles, 3 variables = 2²= 8).
\(as.table\) permite en conjunto con as.data.frame crear la tabla con la combinación de las variables y sus frecuencias. Luego, con \(lapply\), y los argumentos \(relevel\) y \(ref\) se asigna una cuarta columna que tomara como referencia el estado, nivel o categoria no, de tal manera que se expresa todo en terminos del estado contrastante “si”.
tabcon.df <- as.data.frame(as.table(tabcon))
tabcon.df[,-4] <- lapply(tabcon.df[,-4], relevel, ref = "No")
tabcon.df## Cigarrillo Marihuana Alcohol Freq
## 1 Si Si Si 180
## 2 No Si Si 10
## 3 Si No Si 108
## 4 No No Si 90
## 5 Si Si No 50
## 6 No Si No 2
## 7 Si No No 12
## 8 No No No 20“loglineal”.Como ejemplo, se voltea una moneda de 200, 100 veces. También una moneda de 1000, 100 veces, y se crea una tabla de contingencia de 2 x 2. Aquí hay una simulación rápida:
set.seed(1)
m1000 <- sample(c("numero","simbolo"),40,replace = T)
m200 <- sample(c("numero","simbolo"),40,replace = T)
table(m1000,m200)## m200
## m1000 numero simbolo
## numero 9 15
## simbolo 5 11Como las monedas son independientes, supondríamos que los recuentos de celdas serían aproximadamente iguales (alrededor de 10 en cada celda).La independencia permite calcular a aprtir de sus marginales obtener sus probabilidades conjuntas.
En nuestro ejemplo, las probabilidades marginales son de 0,5 para la cara y la cruz de cada moneda. Las probabilidades conjuntas son las probabilidades de un resultado particular en la tabla. Por ejemplo, la probabilidad de que ambas monedas caigan cara es de 0,5 * 0,5 = 0,25. Para obtener un recuento esperado, multiplicamos esa probabilidad por el número de pruebas. En este caso es 40 * 0,25 = 10. En la simulación se encontró 11,15,9 y 5.
Podríamos escribirlo como: \[Conteo~celdas_{i,j}= nπ_i*nπ_j\]
Conteo de celdas en la fila i y la columna j es igual al producto de las probabilidades marginales y el número de observaciones por lo que depende del tamaño de la tabla (variables de fila y columna), siempre que sean independientes. Que si usamos logaritmos y propiedades se obtiene una adición (modelo lineal)
MODELO LOG-LINEAL DE INDEPENDENCIA PARA TABLAS DE CONTINGENCIA DE DOS VIAS
\[logμ_{ij}=logn+logπ_i+logπ_j\] Dos variables independentes, de lo contrario se debe añadir un parametro que explique o permita la interacción al modelo (término de interacción o de orden superior). La interacciones entre variables permiten determinan si dos variables están asociadas y de qué manera.
Volviendo a la encuesta, modelamos la frecuencia como función de las tres variables usando \(glm\). Asumimos una distrbución Poisson, aunque se puedan usar otras útiles para conteo (Modelo de independencia). Con los datos que tenemos parece muy poco probable.
mod1 <- glm(Freq ~ Cigarrillo + Marihuana + Alcohol,
data = tabcon.df, family = poisson)
summary(mod1)##
## Call:
## glm(formula = Freq ~ Cigarrillo + Marihuana + Alcohol, 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 ***
## CigarrilloSi 1.0539 0.1051 10.02 <2e-16 ***
## MarihuanaSi 0.0509 0.0921 0.55 0.58
## AlcoholSi 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: 5Mirando coeficientes y valores p cercanos a 0, excepto para Marihuana, el modelo parece explicar decentemente las variables. Pero para los modelos log-lineales queremos comprobar la desviación residual. Como regla general, nos gustaría que se acercara en valor a los grados de libertad.Se obtuvo 133.63 en 4 grados de libertad (mal ajuste). La estadística de desviación tiene una distribución aproximada de chi-cuadrado, así que usamos la función \(pchisq\). La función \(df.residual\) extrae los grados de libertad. Si estableciendo el parametro lower.tail = F, obtenemos la probabilidad de ver una desviación tan grande o más grande que la que observamos.
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 6.53e-28\[ Ho: Las~frecuencias~satisfacen~el~modelo~log-lineal~independiente\]. ¡Claramente no lo hacen!. Otra posibildad es comparar los valores ajustados con los valores observados. Combinandolos:
cbind(mod1$data, fitted(mod1))## Cigarrillo Marihuana Alcohol Freq fitted(mod1)
## 1 Si Si Si 180 147.5
## 2 No Si Si 10 51.4
## 3 Si No Si 108 140.2
## 4 No No Si 90 48.9
## 5 Si Si No 50 31.9
## 6 No Si No 2 11.1
## 7 Si No No 12 30.4
## 8 No No No 20 10.6El modelo ajustado está muy lejos de los datos observados, en su gran mayoria son muy diferentes y disperos al comparlos directamente. Los coeficientes pueden ser interpretados como probabilidades. Por ejemplo, para marihuana, obtenemos lo siguiente:
exp(coef(mod1)[3])## MarihuanaSi
## 1.05Esas son las probabilidades de usar marihuana porque “no” es la referencia, y al calcularlas probabilidades directamente:
prueba_odd=margin.table(tabcon, margin = 2)/sum(margin.table(tabcon, margin = 2))
prueba_odd[1]/prueba_odd[2]## Si
## 1.05Ajustemos un modelo más complejo que permita que las variables se asocien entre sí, pero que mantenga la misma independemcia de tercera variable. Llamamos a esto asociación homogénea. Esto dice que, por ejemplo, el consumo de alcohol y de marihuana tienen algún tipo de relación, pero independientemente de si probaron o no los cigarrillos. Ajustar este modelo significa que añadimos interacciones para cada combinación de variables por pares. R lo hace fácil con su notación de fórmula. Simplemente ponga paréntesis alrededor de sus variables y añada un exponente de 2. Esto se traduce en “todas las variables y todas las interacciones por pares”
mod2 <- glm(Freq ~ (Cigarrillo + Marihuana + Alcohol)^2,
data = tabcon.df, family = poisson)
summary(mod2)##
## Call:
## glm(formula = Freq ~ (Cigarrillo + Marihuana + Alcohol)^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 ***
## CigarrilloSi -0.307 0.315 -0.97 0.330
## MarihuanaSi -1.645 0.373 -4.41 1.1e-05 ***
## AlcoholSi 1.602 0.240 6.67 2.6e-11 ***
## CigarrilloSi:MarihuanaSi 2.918 0.329 8.88 < 2e-16 ***
## CigarrilloSi:AlcoholSi 0.458 0.332 1.38 0.167
## MarihuanaSi:AlcoholSi -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: 4Este modelo encaja mucho mejor. Fíjese en la desviación residual (1,5132 comparada con un grado de libertad (1). Podemos calcular un valor p si queremos:
\[ Ho: Las~frecuencias~satisfacen~el~modelo~2\].
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 6.53e-28El alto valor p dice que no hay evidencia para rechazar la hipótesis nula . Una vez más podemos comparar los valores ajustados y los observados y ver si coinciden:
cbind(mod1$data, fitted(mod1))## Cigarrillo Marihuana Alcohol Freq fitted(mod1)
## 1 Si Si Si 180 147.5
## 2 No Si Si 10 51.4
## 3 Si No Si 108 140.2
## 4 No No Si 90 48.9
## 5 Si Si No 50 31.9
## 6 No Si No 2 11.1
## 7 Si No No 12 30.4
## 8 No No No 20 10.6Así que el modelo parece encajar (muy bien), pero ¿cómo describimos la asociación entre las variables? ¿Qué efecto tienen las variables entre sí? Para responder a esto miramos los coeficientes de las interacciones. Exponiendo los coeficientes, obtenemos los odds ratios. Por ejemplo, veamos el coeficiente para “marihuanasi:alcoholsi”.
exp(coef(mod1)["MarihuanaSi:AlcoholSi"])## <NA>
## NALos estudiantes que probaron marihuana tienen probabilidades bajas de haber probado el alcohol que son 0.477 veces las probabilidades estimadas para los estudiantes que no probaron la marihuana. Debido a que encajamos en un modelo de asociación homogéneo, la proporción de probabilidades es la misma independientemente de si probaron cigarrillos. Es una buena idea calcular un intervalo de confianza para estas estimaciones de la proporción de probabilidades. Podemos hacerlo con la función \(“confint”\):
mod2$coefficients## (Intercept) CigarrilloSi MarihuanaSi
## 2.915 -0.307 -1.645
## AlcoholSi CigarrilloSi:MarihuanaSi CigarrilloSi:AlcoholSi
## 1.602 2.918 0.458
## MarihuanaSi:AlcoholSi
## -0.739exp(confint(mod2, parm = c("CigarrilloSi:MarihuanaSi","CigarrilloSi:AlcoholSi","MarihuanaSi:AlcoholSi")))## Waiting for profiling to be done...## 2.5 % 97.5 %
## CigarrilloSi:MarihuanaSi 10.092 36.965
## CigarrilloSi:AlcoholSi 0.822 3.037
## MarihuanaSi:AlcoholSi 0.262 0.844Las asociaciones son claras. Las probabilidades de probar la marihuana si probaste los cigarrillos son al alreddor de 12 veces más altas que las probabilidades de probar la marihuana si no probaste los cigarrillos, y viceversa, como se discutia en los primeros pasos.
¿Qué pasa si encajamos un modelo con una interacción de tres vías? -Esto permite que las asociaciones por pares cambien según el nivel de la tercera variable. Este es el modelo saturado ya que tiene tantos coeficientes como celdas en nuestra tabla: 8. Se ajusta perfectamente a los datos. La fórmula “Cigarrillo * Marihuana * Alcohol” significa que encaja en todas las interacciones.
mod3 <- glm(Freq ~ Cigarrillo * Marihuana * Alcohol,
data = tabcon.df, family = poisson)
summary(mod3)##
## Call:
## glm(formula = Freq ~ Cigarrillo * Marihuana * Alcohol, 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 ***
## CigarrilloSi -0.511 0.365 -1.40 0.1618
## MarihuanaSi -2.303 0.742 -3.10 0.0019 **
## AlcoholSi 1.504 0.247 6.08 1.2e-09 ***
## CigarrilloSi:MarihuanaSi 3.730 0.808 4.61 3.9e-06 ***
## CigarrilloSi:AlcoholSi 0.693 0.392 1.77 0.0771 .
## MarihuanaSi:AlcoholSi 0.105 0.813 0.13 0.8969
## CigarrilloSi:MarihuanaSi:AlcoholSi -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: 4.6185e-14 on 0 degrees of freedom
## AIC: 57.6
##
## Number of Fisher Scoring iterations: 3La desviación residual del modelo es básicamente 0 sobre 0 grados de libertad. Además, los recuentos ajustados coinciden con los observados:
cbind(mod3$data, fitted(mod3))## Cigarrillo Marihuana Alcohol Freq fitted(mod3)
## 1 Si Si Si 180 180
## 2 No Si Si 10 10
## 3 Si No Si 108 108
## 4 No No Si 90 90
## 5 Si Si No 50 50
## 6 No Si No 2 2
## 7 Si No No 12 12
## 8 No No No 20 20En igualdad de condiciones, preferimos un modelo más simple. Normalmente no queremos 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.
\[El~modelo~2~encaja~tan~bien~como~el~modelo~3\]
anova(mod2, mod3)## Analysis of Deviance Table
##
## Model 1: Freq ~ (Cigarrillo + Marihuana + Alcohol)^2
## Model 2: Freq ~ Cigarrillo * Marihuana * Alcohol
## Resid. Df Resid. Dev Df Deviance
## 1 1 1.51
## 2 0 0.00 1 1.51pchisq(0.373, df = 1, lower.tail = F)## [1] 0.541No rechazamos la hipótesis nula. Podríamos probar modelos con sólo ciertas interacciones. Por ejemplo, el siguiente modelo:
mod4 <- glm(Freq ~ (Cigarrillo * Marihuana) + (Alcohol * Marihuana), data = tabcon.df, family = poisson)
summary(mod4)##
## Call:
## glm(formula = Freq ~ (Cigarrillo * Marihuana) + (Alcohol * Marihuana),
## 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 ***
## CigarrilloSi 0.087 0.132 0.66 0.510
## MarihuanaSi -1.781 0.367 -4.86 1.2e-06 ***
## AlcoholSi 1.823 0.190 9.57 < 2e-16 ***
## CigarrilloSi:MarihuanaSi 2.866 0.324 8.84 < 2e-16 ***
## MarihuanaSi:AlcoholSi -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: 4Esto encaja con las interacciones para el cigarrillo y la marihuana, y el alcohol y la marihuana, pero no para el cigarrillo * el alcohol. La implicación es que el uso de alcohol y cigarrillos es independiente del otro, controlando el uso de la marihuana. ¿Es este un modelo adecuado? Una vez más podemos realizar una prueba de proporción de probabilidad:
anova(mod4, mod2)## Analysis of Deviance Table
##
## Model 1: Freq ~ (Cigarrillo * Marihuana) + (Alcohol * Marihuana)
## Model 2: Freq ~ (Cigarrillo + Marihuana + Alcohol)^2
## Resid. Df Resid. Dev Df Deviance
## 1 2 3.41
## 2 1 1.51 1 1.9La probabilidad de ver un cambio tan grande en la desviación (3.41) si los modelos no fueran realmente diferentes es remota. Parece haber buenas pruebas de que el modelo de asociación homogénea encaja mucho mejor que el modelo que asume la independencia condicional entre el alcohol y el uso de cigarrillos. Los modelos log-lineales funcionan para tablas más grandes que se extienden en 4 o más dimensiones. Obviamente la interpretación de las interacciones se vuelve mucho más complicada cuando involucran 3 o más variables.
tabcon2 <- array(data = c(180,10, 108,90,50, 2, 12, 20), dim = c(2,2,2), dimnames = list("Marihuana" = c("si","no"), "Alcohol" = c("si","no"),"Cigarrillo" = c("si","no"))); tabcon2## , , Cigarrillo = si
##
## Alcohol
## Marihuana si no
## si 180 108
## no 10 90
##
## , , Cigarrillo = no
##
## Alcohol
## Marihuana si no
## si 50 12
## no 2 20