Se quiere saber si hay propensión a usar la marihuana como consecuencia de usar alcohol y cigarillo. Nos interesa saber cómo los conteos en esta tabla depende de los niveles de las variables categóricas.
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")))ftable(tabcon, row.vars = c("Alc","Cig")) # Razon## Mar S N
## Alc Cig
## S S 180 108
## N 10 90
## N S 50 12
## N 2 20addmargins(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 472options(digits=3)
prop.table(tabcon, margin = c(1,3)) # Proporcion## , , 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.909Con las tablas anteriores primero se observa la relación entre los totales de consumir alcohol y de no consumirlo (388/84) y esto da más o menos una relación de 5/1; lo que significa que por cada 5 que consuman alcohol 1 no lo hace. Esto es una razón.
Cuando se convierte a proporciones los margenes de la tabla se observa que el porcentaje mas alto de los que si consumen alcohol es NN (0.900) los que no usan cigarillos ni marihuana; y de los que no consumen alcohol el porcentaje mas alto es también NN (0.909) los que no usan cigarrillo ni marihuana.
Se ve que de los estudiantes que probaron los cigarillos y el alcohol, el 62,5% también probó marihuana.
tabconC<-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"))); tabconC## , , Cig = S
##
## Mar
## Alc S N
## S 180 108
## N 50 12
##
## , , Cig = N
##
## Mar
## Alc S N
## S 10 90
## N 2 20ftable(tabconC, row.vars = c("Cig", "Alc")) # Razon## Mar S N
## Cig Alc
## S S 180 108
## N 50 12
## N S 10 90
## N 2 20addmargins(tabconC)## , , 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(tabconC, margin = c(1,3)) # Proporciones ## , , 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.909Sí, que los que fuman cigarillo son más propensos a consumir alcohol y/o marihuana, mientras que los que no fuman cigarrilo son menos propensos a consumir alcohol y/o marihuana.
Sí, porque se observa que los que fuman cigarros y consumen marihuana son 230, mientras que los que no fuman cigarros, pero consumen marihuana son 12. La relación es 230:12 aproximadamente 19:1
Si, porque 288 que consumen alcohol fuman cigarillos, mientras que solo 100 consumen alcohol pero no fuman cigarillos. Por otro lado, en cuanto a proporciones también se puede observar que las personas que consumen cigarrillo, marihuana y alcohol es del 63% aproximadamente, mientras que las personas que consumen cigarrillo, marihuana y no consumen alcohol es aproximadamente del 81%. De ahí, que el consumo de cigarrillo está más relacionado al consumo de marihuana sin necesidad de consumir alcohol.
Se usará la funció glm de los modelos lineales generalizados
Tabla de contingencia ahora es un data frame con las tres columnas asociadas a las categorías y sus respectivas frecuencias.
tabcon.df <- as.data.frame(as.table(tabcon))
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 20¿Por qué modelo loglineal?
# Ejemplo de las monedas y su probabilidad...
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 11Acá no hay interacción \[\log \ (\mu_{ij}) = \log \ n + \log \ \pi_{i} + \log \ \pi_{j}\] # valores esperados son igual a 0.
Los modelos loglineales permiten determinar si dos variables están asociadas y de qué manera.
Modelo de independencia
mod0 <- glm(Freq ~ Cig + Mar + Alc,
data = tabcon.df, family = poisson) #poisson quiere decir que estoy modelando conteos
summary(mod0)##
## 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: 5Este modelo está mal,acá se puede observar que hay un problema de sobredispersión, ya que al mirar la deviancia residual y comparar 133.63 con 4 grados de confianza da un número muy grande indicando que los datos tienen más variabilidad de los que se esperan en un modelado de Poisson.Esto es un mal ajuste.
pchisq(deviance(mod0), df = df.residual(mod0), lower.tail = F)## [1] 6.53e-28Primero la Hipótesis nula de esta prueba es que las frecuencias esperadas satisfacen el modelo loglineal. Por tanto, el pchisq al dar un valor muy cercano a 0 ( 6.53e-28) indica que se rechaza la H0 y de ahí que el modelo es malo y no confiable. Entonces este modelo no me sirve para concluir si hay relación entre las tres variables de estudio.
Hallando los valores estimados del modelo. Una forma más intuitiva de investigar el ajuste es comparar los valores ajustados con los valores observados.
cbind(mod0$data, fitted(mod0))## Cig Mar Alc 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.6Acá se observa el problema de sobredispersión. Se observa que el modelo ajustado está muy lejos de los datos observados.
exp(coef(mod0)[3]) #se escogio el estimado de la posicion 3## MarS
## 1.05Este es el odd (cociente de disparidad). Es el odd de usar marihuana (1.05), la probabilidad de usar y no usar marihuana es practimante la misma.
prueba_odd=margin.table(tabcon, margin = 2)/sum(margin.table(tabcon, margin = 2))
prueba_odd[1]/prueba_odd[2]## S
## 1.05\[\log \ (\mu_{ijk}) = \lambda + \lambda^X_{i} + \lambda^Y_{j} + \lambda^Z_{k} + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} + \lambda^{XZ}_{ik}\]
mod1 <- glm(Freq ~ (Cig + Mar + Alc)^2,
data = tabcon.df, family = poisson)
summary(mod1)##
## 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: 4Ahora al mirar la deviancia residual (1.5132) y los grados de libertad (1) son muy parecidos.
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 0.219Con este pchisq ya no rechazo la hipótesis nula, de ahí que este modelo se ajusta y no tiene problema de sobredispersión, predice mejor las cosas a comparación del anterior.Esto se puede observar en la siguiente tabla en donde los valores tienden a ser muy parecidos.
cbind(mod1$data, fitted(mod1))## Cig Mar Alc Freq fitted(mod1)
## 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.44ODD de marihuana Sí y alcohol Sí
exp(coef(mod1)["MarS:AlcS"]) # este es el ODD## MarS:AlcS
## 0.477La relación es aproximadamente 0.5 o 1:2, lo que quiere decir que por cada uno de Sí marihuana 2 Sí alcohol.
mod1$coefficients## (Intercept) CigS MarS AlcS CigS:MarS CigS:AlcS
## 2.915 -0.307 -1.645 1.602 2.918 0.458
## MarS:AlcS
## -0.739exp(confint(mod1, 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.844Las asociaciones son CASI todas muy fuertes, Se observa que las posibilidades de probar la marihuana si se probó el cigarrillo son de al menos 10 veces más altas que las probabilidades de probar marihuana si no se probó el cigarrillo. La puerta para la marihuana es el cigarrillo… El 10 significa que por cada 10 cigarrillo sí 1 marihuana sí.
\[\log \ (\mu_{ijk}) = \lambda + \lambda^X_{i} + \lambda^Y_{j} + \lambda^Z_{k} + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk}\] Acá se quitó la interacción que no generó un p significativo en el anterior:
mod2 <- glm(Freq ~ (Cig * Mar) + (Mar * Alc),
data = tabcon.df, family = poisson)
summary(mod2)##
## Call:
## glm(formula = Freq ~ (Cig * Mar) + (Mar * Alc), 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: 4pchisq(deviance(mod2), df = df.residual(mod2), lower.tail = F)## [1] 0.182Primero la Hipótesis nula de esta prueba es que las frecuencias esperadas satisfacen el modelo loglineal. Por tanto, el pchisq al dar un valor mayor a 0.05 (0.182) indica que NO se rechaza la H0, por tanto, el modelo es bueno y confiable.
Acá se puede observar también que a pesar de que se quitó la interacción entre cigarrillo y alcohol el AIC disminuyó a 57, pero la deviancia residual aumentó un poco ahora al comparar sería 3.4093 con 2 grados de libertad, sin embargo, el valor es cercano.
cbind(mod2$data, fitted(mod2))## Cig Mar Alc Freq fitted(mod2)
## 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.30Al comparar los valores ajustados con los valores observados se observa que son muy parecidos, sin embargo, hay algunas diferencias.
exp(confint(mod2, 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.953Se puede observar por el límite inferior que por cada 10, aproximadamente, que fuman cigarrillo 1 fuma marihuana, y por el límite superior por cada 35, aproximadamente, que fuman cigarrillo 1 fuma marihuana. Por otro lado, se puede observar por el límite inferior (~0.4) que por cada 2, aproximadamente, que fuman marihuana 5 consumen alcohol, y por el límite superior (~1) por cada uno, aproximadamente, que fuma marihuana 1 consume alcohol.
Para finalizar, escojo el 2 modelo que presenta las interacciones entre variables que dieron significativas.
\[\log \ (\mu_{ijk}) = \lambda + \lambda^X_{i} + \lambda^Y_{j} + \lambda^Z_{k} + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk}\]
Ahora corriendo el MODELO SATURADO, el cual tiene la interacción triple
\[\log \ (\mu_{ijk}) = \lambda + \lambda^X_{i} + \lambda^Y_{j} + \lambda^Z_{k} + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} + \lambda^{XZ}_{ik} + \lambda^{XYZ}_{ijk}\]
mod3 <- glm(Freq ~ (Cig * Mar * Alc)^2,
data = tabcon.df, family = poisson)
summary(mod3)##
## Call:
## glm(formula = Freq ~ (Cig * Mar * Alc)^2, 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: 3pchisq(deviance(mod3), df = df.residual(mod3), lower.tail = F)## [1] 0El valor p dio mayor a 0.05 por tanto no se rechaza la hipótesis nula y este modelo se ajusta bien. Al mirar la siguiente tabla los valores esperados y observados son los mismos.
cbind(mod2$data, fitted(mod3))## Cig Mar Alc Freq fitted(mod3)
## 1 S S S 180 180
## 2 N S S 10 10
## 3 S N S 108 108
## 4 N N S 90 90
## 5 S S N 50 50
## 6 N S N 2 2
## 7 S N N 12 12
## 8 N N N 20 20Sin embargo, se espera escoger un modelo más simple, en donde se eliminen las interacciones dobles que no sean significativas.