Chapter 7 REGRESIÓN LOGÍSTICA MÚLTIPLE
Veamos cómo ajustar varios modelos de regresión logística a un conjunto de datos. Tomemos el conjunto de datos Pima.tr y cree una variable categórica Mom donde una mujer con npreg> 0 se clasifica como madre (sin tener en cuenta la posibilidad de que algunos los embarazos pueden no haber resultado en un nacimiento vivo).
Pima.tr <- Pima.tr %>%
mutate(Mom=ifelse(npreg==0,"No","Yes"))
xtabs(~type+Mom,data=Pima.tr)## Mom
## type No Yes
## No 16 116
## Yes 12 56Calculemos la razón de posibilidades(odss-ratio) de la tabla :
\[ OR=\dfrac{\frac{a}{b}}{\frac{c}{d}}=\dfrac{\frac{16}{116}}{\frac{12}{256}}=0.6437\] Tomaremos el recíproco para facilitar la interpretación, \(\frac{1}{OR}=\frac{1}{0.6437}=1.5536\)
Las probabilidades
de tener diabetes tipo II son 1.55 veces mayores para las no madres que para las madres.
Observe que aproximadamente el \(43\%\) de las no madres y el
\(33\%\) de las madres tienen el Tipo II diabetes.
Primero, encuentre un intervalo de confianza para el logaritmo de la razón de posibilidades(logit)
\[ln(OR)\pm z\sqrt{\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d} }\] \[ln(1.5536)\pm (1.96)\sqrt{\frac{1}{16}+\frac{1}{116}+\frac{1}{12}+\frac{1}{56}} \]
\[\left( 0.4406 \pm 0.8316\right) \] \[(-0.3730, 1.2542)\]
Este intervalo de confianza contiene 0, por lo que la variable Mom no es un predictor significativo de la diabetes tipo II. Si prefiere el IC en términos de razón de posibilidades, exponencial.
\[ \left( e^{-0.3730} , e^{1.2542}\right) \]
\[\left( 0.6887, 3.5050\right) \]
El IC incluye el valor 1 (que indica que no hay efecto para la razón de posibilidades).
Usemos R para ajustar el tipo de type~Mom
mod0<- glm(type~1,data=Pima.tr,family="binomial"(link=logit))
mod1<- glm(type~Mom,data=Pima.tr,family="binomial"(link=logit))
summary(mod1)##
## Call:
## glm(formula = type ~ Mom, family = binomial(link = logit), data = Pima.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0579 -0.8876 -0.8876 1.4981 1.4981
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2877 0.3819 -0.753 0.451
## MomYes -0.4406 0.4151 -1.061 0.289
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 256.41 on 199 degrees of freedom
## Residual deviance: 255.31 on 198 degrees of freedom
## AIC: 259.31
##
## Number of Fisher Scoring iterations: 4Anova(mod1,test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: type
## LR Chisq Df Pr(>Chisq)
## Mom 1.1056 1 0.2931confint(mod1)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -1.058545 0.4561131
## MomYes -1.250578 0.3907974Vemos que Mom no es significativa con la prueba de Wald o la razón de verosimilitud prueba. El intervalo de confianza se calcula con una fórmula más compleja que dado aquí, por lo que los resultados no son idénticos. Los signos son opuestos porque mi La tabla usó No como un éxito y el ajuste glm de R usó Sí como un éxito. Hagamos una regresión logística múltiple y hagamos una predicción. Usaré 3 predictores bmi, age y Mom.
mod2 <- glm(type~bmi+age+Mom,data=Pima.tr,family="binomial"(link=logit))
summary(mod2)##
## Call:
## glm(formula = type ~ bmi + age + Mom, family = binomial(link = logit),
## data = Pima.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8477 -0.8045 -0.4907 0.9983 2.3009
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.76192 1.23679 -4.659 3.18e-06 ***
## bmi 0.09703 0.03011 3.223 0.00127 **
## age 0.07830 0.01628 4.809 1.52e-06 ***
## MomYes -0.83471 0.48052 -1.737 0.08237 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 256.41 on 199 degrees of freedom
## Residual deviance: 212.96 on 196 degrees of freedom
## AIC: 220.96
##
## Number of Fisher Scoring iterations: 4Observe que el bmi y la edad son predictores significativos de bmi con valores pendiente positivas. Lo que indica un mayor riesgo a medida que aumenta el bmi o la edad. Mom no es significativo en \(\alpha=0.05\), pero tiene una pendiente negativa que indica que las madres fueron menos probabilidades de ser diabéticas que las no madres.
Si una mujer tiene 40 años, un bmi de 28 y es madre, calculemos su probabilidad de diabetes tipo II.
\[ln\left(\frac{\hat{\pi}}{1-\hat{\pi}}\right)=-5.76192 + 0.09703(28) + 0.07830(40)-0.83471(1) = -0.74779 \]
Su logit negativo indica menos del \(50\%\) de posibilidades de diabetes. Tomando el logit inverso:
\[\frac{e^{-0.74779}}{1+e^{-0.74779}}=0.321 \]
La probabilidad es de aproximadamente el \(32\%\). Ahora hazlo para una persona con las mismas estadísticas, excepto que no sea madre.
\[ ln\left(\frac{\hat{\pi}}{1-\hat{\pi}}\right)=-5.76192 + 0.09703(28) + 0.07830(40)-0.83471(0) = 0.08692\]
Ahora el logit es positivo, por lo que la probabilidad será superior al \(50\%\).
\[\frac{e^{0.08692}}{1+e^{0.08692}}=0.522 \]
Una prueba de razón de verosimilitud que compara el modelo 1 (solo con mamá) y el modelo 2 (con bmi, age y Mom) tendrá \(df=2\) con dos parámetros adicionales, y tendríamos esperar que el segundo modelo sea una mejora significativa.
anova(mod1,mod2,test="LR")## Analysis of Deviance Table
##
## Model 1: type ~ Mom
## Model 2: type ~ bmi + age + Mom
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 198 255.31
## 2 196 212.96 2 42.353 6.355e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que hay una mejora significativa, con \(\chi^2=42.353\), \(df=2\) y \(p=0.0001\).
Tal vez no deberíamos haber creado una variable categórica como Mom, pero solo usar npreg. Encajaré un tercer tipo de type ~ bmi + age + npreg.
mod3 <- glm(type~bmi+age+npreg,data=Pima.tr,family="binomial"(link=logit))
summary(mod3)##
## Call:
## glm(formula = type ~ bmi + age + npreg, family = binomial(link = logit),
## data = Pima.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7413 -0.8235 -0.4918 0.9773 2.2382
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.44462 1.18360 -5.445 5.18e-08 ***
## bmi 0.10761 0.02986 3.604 0.000313 ***
## age 0.05937 0.01817 3.267 0.001086 **
## npreg 0.06508 0.05720 1.138 0.255226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 256.41 on 199 degrees of freedom
## Residual deviance: 214.62 on 196 degrees of freedom
## AIC: 222.62
##
## Number of Fisher Scoring iterations: 4Algo extraño, ya que npreg no es significativo, pero el signo de su estimación es positivo, en lugar de negativo para Mom. Puede haber alguna explicación médica. No soy conciente de.
Por último, suponga que quisiera comparar los modelos 2 y 3. No están anidados, por lo que necesita usar AIC en su lugar. Crearé una tabla AIC para los cuatro modelos (incluidos el modelo nulo).
require(bbmle)## Loading required package: bbmle## Loading required package: stats4##
## Attaching package: 'bbmle'## The following object is masked from 'package:dplyr':
##
## sliceAICtab(mod0,mod1,mod2,mod3,base=TRUE,delta=TRUE,weights=TRUE,sort=TRUE)## AIC dAIC df weight
## mod2 221.0 0.0 4 0.7
## mod3 222.6 1.7 4 0.3
## mod0 258.4 37.5 1 <0.001
## mod1 259.3 38.4 2 <0.001Al AIC parece gustarle un poco más el Model 2 que el Model 3, aunque no hasta cierto punto eso se consideraría sustancial. El modelo 0 y el modelo 1 son muy débiles, con \(\Delta_i>10\) y pesos diminutos Akaike \(w_i<0.001\).
predict(object = mod2, newdata = data.frame(bmi =30.3 ,age=27,Mom="No"))## 1
## -0.7077035