Se trabajará con una encuesta en la que se preguntó a los estudiantes del primer semestre de cierta facultad si habían consumido alcohol, cigarrillos o marihuana.
La respuesta fueron las siguientes:
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")))Usamos la función “ftable” para que la tabla apareza.
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 20Usamos la función “prop.table” para calcular las proporciones de las celdas. La opción margin permite especificar en qué margen queremos trabajar. Por lo cual, se va a calcular las proporciones a través de las columnas (2) a lo largo de las filas (1) para cada capa (3). De ahí el argumento margin =c(1,3)
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.909En la tabla podemos observar que los estudiantes que probaron los cigarrillos y el alcohol, el 62,5% también probó la marihuana.
Y aquellos que no probaron alcohol pero sí probaron cigarillo y marihuana, fue del 80%.
La función addmargins nos proporciona los totales marginales.
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 472Se usará la función glm de los modelos lineales generalizados. La función loglin trabaja directamente con las tablas de contigencia.
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 20mod0 <- glm(Freq ~ cig + Mar + Alc, data = tabcon.df, family = poisson)
summary.glm(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: 5Podemos observar que el modelo parece tener coeficientes altamente significativos y con valores p cercanos a 0.
Pero para modelos loglineales se pretende comprobar la desviación residual (DR). El valor de DR fue de 133.63 en 4 grados de libertad. Esto indica un mal ajuste.
Por lo tanto, se utilizará la función pchisq para poder calcular un p con una distribución aproximada de chi-cuadrado. La función df.residual extrae los grados de libertad. Si establecemos la función lower.tail = F, obtenemos la probabilidad de ver una desviación.
pchisq(deviance(mod0), df = df.residual(mod0), lower.tail = F)## [1] 6.53e-28La Ho de esta prueba es que las frecuencias esperadas satisfacen el modelo loglineal dado.
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.6Para saber el ajuste del modelo, es necesario comparar los valores ajustados con los valores observados. Por eso, se combinó los datos originales.
Finalmente, con el modelo ajustado, se observó este se encuentra muy lejos de los datos observados. Por ejemplo, observamos 10 esutidantes que probaron solo marihuana y alcohol y nuestro modelo predice unos 51.
exp(coef(mod0)[3])## MarS
## 1.05Se exponenció el coeficiente de la marihuana obteniendo un valor de 1.05.
Se calculará entonces las probabilidades directamente mediante la siguiente función:
prueba_odd = margin.table(tabcon, margin = 2)/sum(margin.table(tabcon, margin = 2))
prueba_odd[1]/prueba_odd[2]## S
## 1.05Según el modelo, las posbilidades (odd) de usar marihuana son de aproximadamente 1,05 a 1, independientemente de que se haya probado alcohol o los cigarrillos.
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: 4Este modelo parece ajustarse mucho mejor que el anterior. La desviación residual en este caso es (1,51) comparada con los grados de libertad (1). Podemos calcular un valor p si queremos:
pchisq(deviance(mod1), df = df.residual(mod1), lower.tail = F)## [1] 0.219El valor alto de p dice que no tenemos pruebas suficientes para rechazar la hipótesis nula de que las frecuencias esperadas satisfacen nuestro modelo. Pero, podemos comparar los valores ajustados y los observados y ver si coinciden:
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.44exp(coef(mod1)["MarS:AlcS"])## MarS:AlcS
## 0.477La asociación entre las variables se usó los coeficientes de la interacciónes exponenciándolos.
Por lo cual, se puede decir entonces, que los estudiantes que probaron la marihuana tienen posibilidades estimadas de haber probado el alcohol que es aproximadamente 0,5 veces la posibilidad estimada para los estudiantes que no probaron marihuana.
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.739Estos coeficientes nos pueden ayudar para calcular un intervalo de confianza para estimaciones de odds ratios (posibilidades). Usamos la función de confin
exp(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.844Vemos que las posibilidades de probar la marihuana si se probó los cigarillos son de al menos 10 veces más altas que las probabilidades de probar la marihuana si no se probó los cigarillos.