Generando arreglo
tabcon_hum <- 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")))
tabcon_hum## , , Alc = S
##
## Mar
## Cig S N
## S 180 108
## N 10 90
##
## , , Alc = N
##
## Mar
## Cig S N
## S 50 12
## N 2 20Aplanando la tabla
ftable(tabcon_hum, row.vars = c("Alc","Cig"))## Mar S N
## Alc Cig
## S S 180 108
## N 10 90
## N S 50 12
## N 2 20Totales marginales
addmargins(tabcon_hum)## , , 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 472Proporciones entre variables
options(digits=3)
prop_alc <- prop.table(tabcon_hum, margin = c(1,3))
prop_alc## , , 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.909Creando una tabla que rompe por consumo de cigarrillo y que se ve las tablas 2*2 para las otras dos variables
tabcon_hum_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")))
tabcon_hum_cig## , , Cig = S
##
## Mar
## Alc S N
## S 180 108
## N 50 12
##
## , , Cig = N
##
## Mar
## Alc S N
## S 10 90
## N 2 20Proporciones
options(digits=3)
prop_cig <- prop.table(tabcon_hum_cig, margin = c(1,3))
prop_cig## , , 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.909tabcon_hum_alc1 <- array(data = c(180,108, 10,90,50, 12, 2, 20),
dim = c(2,2,2),
dimnames = list("Mar" = c("S","N"),
"Cig" = c("S","N"),
"Alc" = c("S","N")))
tabcon_hum_alc1## , , Alc = S
##
## Cig
## Mar S N
## S 180 10
## N 108 90
##
## , , Alc = N
##
## Cig
## Mar S N
## S 50 2
## N 12 20Proporciones
options(digits=3)
prop_alc.cig <- prop.table(tabcon_hum_alc1, margin = c(1,3))
prop_alc.cig## , , Alc = S
##
## Cig
## Mar S N
## S 0.947 0.0526
## N 0.545 0.4545
##
## , , Alc = N
##
## Cig
## Mar S N
## S 0.962 0.0385
## N 0.375 0.6250Comparando todas las probabilidades de interés
prop_alc## , , 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.909prop_cig## , , 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.909prop_alc.cig## , , Alc = S
##
## Cig
## Mar S N
## S 0.947 0.0526
## N 0.545 0.4545
##
## , , Alc = N
##
## Cig
## Mar S N
## S 0.962 0.0385
## N 0.375 0.6250¿Percibe alguna relación?
Se percibe que el alcohol no tiene relación con el uso de la marihuana ya que las probabilidades del uso de la marihuana son muy parecidas cuando consume o no alcohol, mientras que si se ve una clara interacción entre el uso de cigarrillo y el uso de marihuana ya que las probabilidades son contrastantes, es decir, que las probabilidades más altas de que fume marihuana están en la tabla que corresponde al uso de cigarrillos y del mismo modo las probabilidades más altas de que no use marihuana se encuentran en la tabla correspondiente al no uso de cigarrillos, esta es otra razón más de que el alcohol no tiene influencia en el uso de marihuana.
¿Aquellos que fuman cigarrillos son más propensos al uso de la marihuana?
Si, como se explicó anteriormente el uso de cigarrillos tiene una fuerte relación con el uso de marihuana, vemos que la probabilidad de que consuman marihuana cuando si usan cigarrillo es del 62,5% para los que consumen alcohol y 80,6% para los que no consumen alcohol, y la probabilidad de que no usen marihuana si no usan cigarrillos es del 90% sin importar el consuman de alcohol.
¿Aquellos que toman alcohol son más propensos al uso de cigarrillos?
No, analizando la última tabla de proporciones confirmamos lo que habíamos percibido con las dos primeras, no hay ninguna relación fuerte entre el consumo de alcohol y el uso de cigarrillos, tampoco entre el consumo de alcohol y el uso de marihuana. En ninguna de las tablas se puede determinar un patrón ya que las probabilidades varían mucho sin importar que si consuma alcohol o que no lo haga.
tabcon_hum.df <- as.data.frame(as.table(tabcon_hum))
tabcon_hum.df[,-4] <- lapply(tabcon_hum.df[,-4], relevel, ref = "N")
tabcon_hum.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_hum.df, family = poisson)
summary(mod0)##
## Call:
## glm(formula = Freq ~ Cig + Mar + Alc, family = poisson, data = tabcon_hum.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: 5Modelo desajustado a causa de la deviancia y el valor de AIC muy alto.
pchisq(deviance(mod0), df = df.residual(mod0), lower.tail = F)## [1] 6.53e-28P-valor muy bajo lo que quiere decir que se rechaza \(H_0\) por lo que las frecuencias observadas no se ajustan las frecuencias estimadas por el modelo.
Comparando valores estimados con los 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.6Se evidencia gran diferencia entre los valores. Suponiendo que el modelo sirviera se realiza un primer análisis del “odd” para la marihuana
exp(coef(mod0)[3])## MarS
## 1.05prueba_odd = margin.table(tabcon_hum, margin = 2)/sum(margin.table(tabcon_hum, margin = 2))
prueba_odd[1]/prueba_odd[2]## S
## 1.05Se genera un segundo modelo denominado modelo de asociación homogénea, en donde se tiene en cuenta las interacciones para cada combinación de variables por pares.
mod1_hum <- glm(Freq ~ (Cig + Mar + Alc)^2,
data = tabcon_hum.df, family = poisson)
summary(mod1_hum)##
## Call:
## glm(formula = Freq ~ (Cig + Mar + Alc)^2, family = poisson, data = tabcon_hum.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: 4Vemos que mejoro bastante ya que la devianza residual se acerca mucho más al número de grados de libertad, además el valor del AIC es mucho más bajo.
pchisq(deviance(mod1_hum), df = df.residual(mod1_hum), lower.tail = F)## [1] 0.219Con este p-valor no se puede rechazar \(H_0\) lo que quiere decir que las frecuencias estimadas se logran ajustar a las frecuencias esperadas.
cbind(mod1_hum$data, fitted(mod1_hum))## Cig Mar Alc Freq fitted(mod1_hum)
## 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.44Ahora coinciden mucho más los valores observados con los estimados por el modelo, por lo que ya es más viable hacer el análsis con este modelo.
Se analiza el “odd ratio” para la interacción “Mar_S:Alc_S”
exp(coef(mod1_hum)["MarS:AlcS"])## MarS:AlcS
## 0.477Es buena idea calcular un intervalo de confianza para la estimación de los odds ratios
mod1_hum$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_hum, 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.844Se evidencias fuertes interaccione entre las variables, y se determina que hay 10 veces más probabilidad que los que fuman cigarrillo fumen marihuana a que los que no fumen cigarrillo fumen marihuana.
Es un modelo que tiene en cuenta todas las interacciones entre las variables.
mod2_hum <- glm(Freq ~ Cig*Mar*Alc,
data = tabcon_hum.df, family = poisson)
summary(mod2_hum)##
## Call:
## glm(formula = Freq ~ Cig * Mar * Alc, family = poisson, data = tabcon_hum.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: 3No se evidencia que la interacción triple sea necesaria.
Comparando valores estimados con los observados
cbind(mod2_hum$data, fitted(mod2_hum))## Cig Mar Alc Freq fitted(mod2_hum)
## 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 20Los valores estimados coinciden perfectamente con los valores observados, esto es característico del modelo saturado.Se prefiere un modelo más simple, sin embargo es útil ajustar un modelo saturado para verificar que la interacción de orden superior no es significativa.
Para verificar que un modelo se ajusta tan bien como el modelo saturado se realiza una prueba de la razón de verosimilidad.
anova(mod1_hum,mod2_hum)## 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.51pchisq(1.51, df = 1, lower.tail = F)## [1] 0.219Con ese p-valor no se rechaza \(H_0\) lo que quiere decir que aunque el modelo saturado coincida perfectamente con los datos observados es mejor el modelo mas sencillo.