Trabajo Final Estadística 2
Lucas Rodriguez Saa
5/12/2022
Ejercicio 1: ANOVA
El Departamento de nutrición de la Ciudad desea analizar los efectos del tipo de harina y del porcentaje de endulzamiento sobre los atributos físicos de un determinado tipo de torta, a la venta actualmente. Para ello, se diseñó un experimento utilizando dos tipos de harina (multipropósito y harina para tortas), y distintos porcentajes de endulzamiento, en cuatro niveles diferentes. Los siguientes datos corresponden a la información de la densidad específica de muestras de tortas preparadas con las distintas combinaciones posibles. Se prepararon 3 tortas con cada una de las combinaciones posibles
- ¿Existen diferencias estadísticamente significativas debidas al tipo de harina utilizada?
- ¿Existen diferencias estadísticamente significativas debidas al nivel de endulzamiento utilizado?
- ¿Existe interacción estadísticamente significativa entre los factores considerados?
- Realizar los gráficos de medias e interpretar.
Resolución Ejercicio 1
En primer lugar, generamos los datos que se muestran en la tabla de la consigna modelados para los fines de este trabajo
densidad <- c(.90, .87, .90, .86, .89, .91,
.93, .88, .87, .79, .82, .80, # general
.91, .90, .80, .88, .82, .83,
.86, .85, .80, .86, .85, .85) # torta
tipo_harina <- c(rep("general", 12),
rep("para_tortas", 12))
conc_azucar <- c(rep("0", 3), rep("50", 3),
rep("75", 3), rep("100", 3),
rep("0", 3), rep("50", 3),
rep("75", 3), rep("100", 3))
data_tortas <- data.frame(densidad, tipo_harina, conc_azucar)
data_tortas$tipo_harina <- as.factor(data_tortas$tipo_harina)
data_tortas$conc_azucar <- factor(data_tortas$conc_azucar,
levels = c("0","50","75","100"))A continuación, para responder a las primeras consignas se realizan
test de ANOVA para las variables tipo_harina y
conc_azucar.
Consigna 1: ¿Existen diferencias estadísticamente significativas debidas al tipo de harina utilizada?
Para estos tests vamos a trabajar con un nivel de confianza de \(95\%\), es decir con un \(\alpha=0.05\). Para realizar esta prueba tenemos que tener en claro cuáles son las hipótesis:
- \(H_0: \mu_1 = \mu_2 = \mu_n\)
- \(H_1:\) al menos una de las medias es diferente
model_anova <- aov(densidad ~ tipo_harina, data = data_tortas)
summary(model_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## tipo_harina 1 0.00184 0.001837 1.16 0.293
## Residuals 22 0.03486 0.001584Obtenemos un valor P de \(0.293\), siendo mayor a nuestro nivel de significancia de \(\alpha=0.05\) no rechazamos la hipótesis nula de igualdad de medias.
Consigna 2: ¿Existen diferencias estadísticamente significativas debidas al nivel de endulzamiento utilizado?
model_anova1 <- aov(densidad ~ conc_azucar, data = data_tortas)
summary(model_anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## conc_azucar 3 0.008713 0.002904 2.076 0.136
## Residuals 20 0.027983 0.001399Obtenemos un valor P de \(0.136\), siendo mayor a nuestro nivel de significancia de \(\alpha=0.05\) no rechazamos la hipótesis nula de igualdad de medias.
No parece haber evidencia estadísticamente significativa para rechazar la hipótesis de igualdad de medias. Tanto el Tipo de Harina como la Concentración de azúcar no parecen tener efecto sobre las densidades de las tortas.
Consigna 3: ¿Existe interacción estadísticamente significativa entre los factores considerados?
El siguiente paso es realizar un test de ANOVA con interacción para las variables en estudio.
model_anova2 <- aov(densidad ~ tipo_harina*conc_azucar, data = data_tortas)
summary(model_anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## tipo_harina 1 0.001838 0.001838 1.838 0.1941
## conc_azucar 3 0.008713 0.002904 2.904 0.0670 .
## tipo_harina:conc_azucar 3 0.010146 0.003382 3.382 0.0442 *
## Residuals 16 0.016000 0.001000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En este caso, obtenemos un valor P de \(0.0442\), siendo menor a nuestro nivel de significancia de \(\alpha=0.05\) rechazamos la hipótesis nula de no interacción entre los factores. Pareciera que la densidad de las tortas se ven afectadas por la interacción entre el Tipo de Harina y la Concentración de Azúcar.
Consigna 4: Realizar los gráficos de medias e interpretar.
En este sentido, también resulta útil explorar los datos con un Gráfico de Interacción:
with(data_tortas, interaction.plot(conc_azucar, tipo_harina, densidad,
fun = mean, col = c("red2", "blue4"),
legend = TRUE,
main = "Interaction Plot"))
Este gráfico nos sugiere que hay un efecto de interacción dado que los segmentos de cada línea no son paralelos y se cruzan. También creamos una tabla de datos para observar las medias (como en el gráfico previo) y los desvíos.
data_tortas %>%
group_by(tipo_harina, conc_azucar) %>%
get_summary_stats(densidad, type = "mean_sd")| tipo_harina | conc_azucar | variable | n | mean | sd |
|---|---|---|---|---|---|
| general | 0 | densidad | 3 | 0.890 | 0.017 |
| general | 50 | densidad | 3 | 0.887 | 0.025 |
| general | 75 | densidad | 3 | 0.893 | 0.032 |
| general | 100 | densidad | 3 | 0.803 | 0.015 |
| para_tortas | 0 | densidad | 3 | 0.870 | 0.061 |
| para_tortas | 50 | densidad | 3 | 0.843 | 0.032 |
| para_tortas | 75 | densidad | 3 | 0.837 | 0.032 |
| para_tortas | 100 | densidad | 3 | 0.853 | 0.006 |
Además, utilizamos el siguiente gráfico de medias que resume los elementos anteriores.
ggline(data_tortas, x = "conc_azucar", y = "densidad",
color = "tipo_harina",
add = c("mean_se", "jitter"),
order = c(0, 50, 75, 100),
palette = c("#00AFBB", "#E7B800"),
ylab = "Densidad", xlab = "Concentración de Azúcar")
Parece haber diferencias entre las medias pero solamente en un subgrupo particular (Harina general y 100 de concentración de azúcar).
Finalmente, hacemos el chequeo de supuestos sobre nuestros datos para comprobar la validez de los resultados y la aplicación del test de ANOVA.
El shapiro_test() nos permite verificar el primer
supuesto de Normalidad con las siguientes hipótesis:
- \(H_0:\) los datos provienen de una distribución Normal
- \(H_1:\) los datos no provienen de una distribución Normal
## Shapiro test para normalidad
data_tortas %>%
group_by(tipo_harina, conc_azucar) %>%
shapiro_test(densidad)| tipo_harina | conc_azucar | variable | statistic | p |
|---|---|---|---|---|
| general | 0 | densidad | 0.7500000 | 0.0000000 |
| general | 50 | densidad | 0.9868421 | 0.7804408 |
| general | 75 | densidad | 0.8709677 | 0.2982759 |
| general | 100 | densidad | 0.9642857 | 0.6368868 |
| para_tortas | 0 | densidad | 0.8175676 | 0.1571668 |
| para_tortas | 50 | densidad | 0.8709677 | 0.2982759 |
| para_tortas | 75 | densidad | 0.8709677 | 0.2982759 |
| para_tortas | 100 | densidad | 0.7500000 | 0.0000000 |
Con un P-value mayor a \(\alpha =
0.05%\) no se rechaza \(H_0\),
es decir que se acepta que los datos provienen de una distribución
Normal. Esto ocurre para la mayoría de los grupos, con lo cual aceptamos
ese resultado. Aunque también podemos testear la Normalidad de forma
gráfica a partir de un `ggqqplot() que, en este caso, nos
servirá de apoyo para validar la decisión.
ggqqplot(data_tortas, "densidad", facet.by = "tipo_harina*conc_azucar")
El qqplot permite comparar gráficamente los cuantiles de una Normal con los cuantiles de una muestra, si la muestra sigue una distribución Normal, los cuantiles deberían coincidir, es decir, los puntitos deberían caer sobre la línea o en la zona sombreada.
Se observa para todos los grupos que los datos no parecen alejarse de los cuantiles teóricos de la Normal y se encuentran dentro del área sombreada. Esto nos lleva a pensar que el resultado del Test de Shapiro esté relacionado con la poca cantidad de casos (3) en cada subgrupo.
El siguiente supuesto a chequear es el de Homocedasticidad, es decir que las varianza son constantes. Para evaluar este comportamiento trabajaremos con las siguientes hipótesis:
- \(H_0:\) las varianzas son iguales
- \(H_1:\) las varianzas no son iguales
## levene test para homocedasticidad
data_tortas %>% levene_test(densidad ~ tipo_harina*conc_azucar)| df1 | df2 | statistic | p |
|---|---|---|---|
| 7 | 16 | 0.4129721 | 0.8803457 |
El resultado nos muestra un valor \(P =0.880\), mayor a nuestro \(\alpha=0.05\), por lo que no rechazamos la hipótesis nula de igualdad de las varianzas
Ejercicio 2: Regresión Lineal
En un estudio que investiga el papel de las armas de fuego en el aumento de la tasa de homicidios de Detroit, se recogieron datos para los años 1961-1973. La variable de respuesta (la tasa de homicidios) y las variables potencialmente predictoras del comportamiento de ésta, se describen a continuación:
| Variable | Descripción |
|---|---|
| FTP | Cantidad de policias full time cada 100.000 habitantes |
| UEMP | Porcentaje de desempleados |
| M | Cantidad de trabajadores en la Industria Manufacturera (en miles) |
| LIC | Cantidad de licencias de portación de armas cada 100.000 habitantes |
| GR | Cantidad de armas de fuego registradas cada 100.000 habitantes |
| CLEAR | Porcentaje de homicidios resueltos con arresto del responsable |
| W | Cantidad de hombres blancos en la población |
| NMAN | Cantidad de trabajadores fuera de la industria manufacturera |
| G | Cantidad de trabajadores del gobierno |
| HE | Ingreso promedio por hora |
| WE | Ingreso promedio semanal |
| H | Tasa de homicidios cada 100.000 habitantes |
Se requiere construir una ecuación de regresión que relacione la Tasa de Homicidios (variable H), con el resto de las variables disponibles, y determinar si estas variables son útiles para predecir la Tasa de Homicidios.
- Cuál es el nivel de asociación lineal de las variables predictoras con la variable H? Comentar.
- Realizar una regresión lineal múltiple, seleccionando los mejores predictores entre las variables independientes disponibles, utilizando un método de selección automática. Describir el proceso de selección automática utilizado. (Sug. Considerar como probabilidad de entrada 0.10 y de salida del modelo 0.15.)
- ¿Qué información da el coeficiente de determinación?
- ¿Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámetros?
- Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
- Realizar un análisis de los residuos del modelo para evaluar el cumplimiento de los supuestos. Para esto, realizar gráficos de los residuos con el valor predicho.
- Analizar la colinealidad de las variables predictoras presentes en la ecuación.
- Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso.
Resolución Ejercicio 2
Para comenzar, levantamos e imprimimos los datos junto a una primera vista de los datos
detroit <- read_sav("2_REG LINEAL/4_p301.sav")
head(detroit)| YEAR | FTP | UNEMP | M | LIC | GR | CLEAR | W | NMAN | G | HE | WE | H |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1961 | 260.35 | 11.0 | 455.5 | 178.15 | 215.98 | 93.4 | 558724 | 538.1 | 133.9 | 2.98 | 117.18 | 8.600000 |
| 1962 | 269.80 | 7.0 | 480.2 | 156.41 | 180.48 | 88.5 | 538584 | 547.6 | 137.6 | 3.09 | 134.02 | 8.900000 |
| 1963 | 272.04 | 5.2 | 506.1 | 198.02 | 209.57 | 94.4 | 519171 | 562.8 | 143.6 | 3.23 | 141.68 | 8.520001 |
| 1964 | 272.96 | 4.3 | 535.8 | 222.10 | 231.67 | 92.0 | 500457 | 591.0 | 150.3 | 3.33 | 147.98 | 8.890000 |
| 1965 | 272.51 | 3.5 | 576.0 | 301.92 | 297.65 | 91.0 | 482418 | 626.1 | 164.3 | 3.46 | 159.85 | 13.070000 |
| 1966 | 261.34 | 3.2 | 601.7 | 391.22 | 367.62 | 87.4 | 465029 | 659.8 | 179.5 | 3.60 | 157.19 | 14.570000 |
summary(detroit)## YEAR FTP UNEMP M
## Min. :1961 Min. :260.4 Min. : 3.200 Min. :455.5
## 1st Qu.:1964 1st Qu.:269.8 1st Qu.: 3.900 1st Qu.:535.8
## Median :1967 Median :273.0 Median : 5.200 Median :569.3
## Mean :1967 Mean :304.5 Mean : 5.792 Mean :556.4
## 3rd Qu.:1970 3rd Qu.:341.4 3rd Qu.: 7.100 3rd Qu.:596.9
## Max. :1973 Max. :390.2 Max. :11.000 Max. :613.5
## LIC GR CLEAR W
## Min. : 156.4 Min. : 180.5 Min. :58.90 Min. :359647
## 1st Qu.: 222.1 1st Qu.: 231.7 1st Qu.:73.90 1st Qu.:401518
## Median : 583.2 Median : 616.5 Median :87.40 Median :448267
## Mean : 537.5 Mean : 545.7 Mean :81.45 Mean :453354
## 3rd Qu.: 794.9 3rd Qu.: 750.4 3rd Qu.:91.00 3rd Qu.:500457
## Max. :1131.2 Max. :1029.8 Max. :94.40 Max. :558724
## NMAN G HE WE
## Min. :538.1 Min. :133.9 Min. :2.910 Min. :117.2
## 1st Qu.:591.0 1st Qu.:150.3 1st Qu.:3.230 1st Qu.:141.7
## Median :686.2 Median :187.5 Median :3.600 Median :157.2
## Mean :673.9 Mean :185.8 Mean :3.948 Mean :170.0
## 3rd Qu.:755.3 3rd Qu.:223.8 3rd Qu.:4.470 3rd Qu.:178.7
## Max. :819.8 Max. :230.9 Max. :5.760 Max. :258.0
## H
## Min. : 8.52
## 1st Qu.: 8.90
## Median :21.36
## Mean :25.13
## 3rd Qu.:37.39
## Max. :52.33str(detroit)## tibble [13 × 13] (S3: tbl_df/tbl/data.frame)
## $ YEAR : num [1:13] 1961 1962 1963 1964 1965 ...
## ..- attr(*, "label")= chr "Year"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ FTP : num [1:13] 260 270 272 273 273 ...
## ..- attr(*, "label")= chr "FTP"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ UNEMP: num [1:13] 11 7 5.2 4.3 3.5 ...
## ..- attr(*, "label")= chr "UNEMP"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ M : num [1:13] 456 480 506 536 576 ...
## ..- attr(*, "label")= chr "M"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ LIC : num [1:13] 178 156 198 222 302 ...
## ..- attr(*, "label")= chr "LIC"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ GR : num [1:13] 216 180 210 232 298 ...
## ..- attr(*, "label")= chr "GR"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ CLEAR: num [1:13] 93.4 88.5 94.4 92 91 ...
## ..- attr(*, "label")= chr "CLEAR"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ W : num [1:13] 558724 538584 519171 500457 482418 ...
## ..- attr(*, "label")= chr "W"
## ..- attr(*, "format.spss")= chr "F12.0"
## ..- attr(*, "display_width")= int 12
## $ NMAN : num [1:13] 538 548 563 591 626 ...
## ..- attr(*, "label")= chr "NMAN"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ G : num [1:13] 134 138 144 150 164 ...
## ..- attr(*, "label")= chr "G"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ HE : num [1:13] 2.98 3.09 3.23 3.33 3.46 ...
## ..- attr(*, "label")= chr "HE"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ WE : num [1:13] 117 134 142 148 160 ...
## ..- attr(*, "label")= chr "WE"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9
## $ H : num [1:13] 8.6 8.9 8.52 8.89 13.07 ...
## ..- attr(*, "label")= chr "H"
## ..- attr(*, "format.spss")= chr "F9.2"
## ..- attr(*, "display_width")= int 9Consigna 1: Cuál es el nivel de asociación lineal de las variables predictoras con la variable H? Comentar.
Vemos que todas las variables son númericas, así que procedemos a realizar un análisis acerca de la correlación entre las variables. Descartamos de nuestro análisis la variable YEAR.
M <-cor(detroit[2:13])
M## FTP UNEMP M LIC GR CLEAR
## FTP 1.0000000 0.292549420 0.4182908 0.5685212 0.70236680 -0.9742599
## UNEMP 0.2925494 1.000000000 -0.6517190 -0.1671449 -0.03287467 -0.3056631
## M 0.4182908 -0.651718959 1.0000000 0.6983145 0.62778120 -0.4292949
## LIC 0.5685212 -0.167144946 0.6983145 1.0000000 0.90379302 -0.5548320
## GR 0.7023668 -0.032874670 0.6277812 0.9037930 1.00000000 -0.6831611
## CLEAR -0.9742599 -0.305663099 -0.4292949 -0.5548320 -0.68316112 1.0000000
## W -0.8799107 0.090013214 -0.7652255 -0.7815460 -0.84782368 0.8821079
## NMAN 0.8818185 -0.038947279 0.7503422 0.7848606 0.84329238 -0.8915961
## G 0.8793086 0.007512467 0.7097715 0.8037000 0.86242889 -0.8934356
## HE 0.9365055 0.230809707 0.4536647 0.4216629 0.57352189 -0.9574269
## WE 0.9222516 0.131281687 0.5023350 0.3910047 0.55589549 -0.9362843
## H 0.9640610 0.210141728 0.5464227 0.7262896 0.81628728 -0.9684603
## W NMAN G HE WE H
## FTP -0.87991074 0.88181847 0.879308571 0.9365055 0.9222516 0.9640610
## UNEMP 0.09001321 -0.03894728 0.007512467 0.2308097 0.1312817 0.2101417
## M -0.76522553 0.75034215 0.709771506 0.4536647 0.5023350 0.5464227
## LIC -0.78154603 0.78486061 0.803700034 0.4216629 0.3910047 0.7262896
## GR -0.84782368 0.84329238 0.862428887 0.5735219 0.5558955 0.8162873
## CLEAR 0.88210786 -0.89159606 -0.893435576 -0.9574269 -0.9362843 -0.9684603
## W 1.00000000 -0.99537633 -0.985670491 -0.8622352 -0.8585848 -0.9469213
## NMAN -0.99537633 1.00000000 0.990035941 0.8695054 0.8555694 0.9559347
## G -0.98567049 0.99003594 1.000000000 0.8572888 0.8264266 0.9580545
## HE -0.86223522 0.86950540 0.857288766 1.0000000 0.9827635 0.9134063
## WE -0.85858484 0.85556944 0.826426612 0.9827635 1.0000000 0.8881526
## H -0.94692129 0.95593473 0.958054543 0.9134063 0.8881526 1.0000000Para un análisis más claro, graficaremos las correlaciones
corrplot::corrplot(M, type="upper", order="original",
addCoef.col = "black", method = "color",
col=brewer.pal(n=8, name="RdYlBu"))
Este gráfico nos muestra con números en cada celda el Coeficiente de Correlación R pareando todas las variables en estudio junto a un color de referencia.
Se observa que para la variable de respuesta H (Tasa de homicidios cada 100.000 habitantes) la mayoría de las variables tienen una correlación cercana a 1 o -1, entendemos que estaríamos en condiciones de descartar para el próximo análisis de regresión múltiple a las variables con un R más bajo que no serían optimas candidatas como variables predictoras, éstas son: UNEMP (Porcentaje de desempleados) con un \(R=0.210\) y M (Cantidad de trabajadores en la Industria Manufacturera) con un \(R=0.546\).
De forma similar, analizaremos el siguiente gráfico de dispersión para todas las variables en análisis:
pairs(detroit[2:13], pch = 19, lower.panel = NULL)
Para un análisis más preciso, utilizaremos los diagramas de dispersión de H con las variables predictoras
pivot_detroit <- detroit[2:13] %>%
pivot_longer(-H, names_to = "var", values_to = "value")
ggplot(pivot_detroit, aes(x = value, y = H)) +
geom_point() +
stat_smooth(method = "lm")+
facet_wrap(~ var, scales = "free")
Realizando un examen visual de los diagramas de dispersión, pareciera que las variables UNEMP y M no seguirían un patrón lineal de relación respecto a H. Debido a esto y al resultado del coeficiente R, podríamos pensar en excluir las variables UNEMP y M de nuestro modelo.
Consigna 2: Realizar una regresión lineal múltiple, seleccionando los mejores predictores entre las variables independientes disponibles, utilizando un método de selección automática. Describir el proceso de selección automática utilizado. (Sug. Considerar como probabilidad de entrada 0.10 y de salida del modelo 0.15.)
Siguiendo con las consignas, y en busca de determinar la ecuación de regresión múltiple seleccionando los mejores predictores, generamos un modelo que incluya todas las variables en estudio.
model_detroit <- lm(H ~ FTP + UNEMP + M + LIC + GR +
CLEAR + W + NMAN + G + HE + WE,
data = detroit)
summary(model_detroit)##
## Call:
## lm(formula = H ~ FTP + UNEMP + M + LIC + GR + CLEAR + W + NMAN +
## G + HE + WE, data = detroit)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.004254 -0.037388 0.252667 -0.174622 -0.169070 0.170520 -0.043951 0.036741
## 9 10 11 12 13
## -0.041456 -0.003869 -0.044030 0.003702 0.055010
## attr(,"label")
## [1] "H"
## attr(,"format.spss")
## [1] "F9.2"
## attr(,"display_width")
## [1] 9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.227e+02 6.681e+01 -1.837 0.3173
## FTP 4.011e-02 1.608e-02 2.495 0.2427
## UNEMP 3.980e-01 3.643e-01 1.093 0.4719
## M -6.186e-02 2.472e-02 -2.502 0.2420
## LIC 1.916e-02 1.856e-03 10.326 0.0615 .
## GR -2.063e-03 1.494e-03 -1.381 0.3989
## CLEAR -6.854e-02 8.864e-02 -0.773 0.5810
## W 1.191e-04 7.825e-05 1.522 0.3700
## NMAN 8.545e-02 4.327e-02 1.975 0.2984
## G 1.495e-01 7.022e-02 2.129 0.2796
## HE -5.890e+00 2.571e+00 -2.291 0.2620
## WE 2.829e-01 6.329e-02 4.470 0.1401
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4042 on 1 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9994
## F-statistic: 1792 on 11 and 1 DF, p-value: 0.01842De este modelo vemos que los valores P de las variables predictoras son mayores a nuestro \(\alpha=0.05\), por lo que no son estadísticamente significativas. Sin embargo, nuestro \(R^2=0.9999\) nos indica que la variación de H es explicada o determinada por nuestro modelo de regresión, aunque al agregar más variables al modelo de regresión el \(R^2\) aumentará siempre, incluso si esas nuevas variables son sólo ruido. Para subsanar este inconveniente tenemos un Coficiente de Determinación Ajustado \(R^2=0.9994\) que nos indicaría la proporción de la variación de H estaría explicada por el modelo, ajustando para la cantidad de variables que se utilizan.
Por otro lado, se nos presenta el resultado del Test de ANOVA para la Regresión Lineal con un estadístico de prueba F y un \(p.value=0.01842\). Consideremos las hipótesis de prueba:
- \(H_0:\) la regresión no es significativa
- \(H_1:\) la regresión es significativa
Con un valor \(P < \alpha = 0.05\) rechazamos la hipótesis nula, es decir que la regresión resulta significativa y continuamos con el proceso.
Esto nos lleva a pensar que hay variables intervinientes que no permiten avanzar u ocultan otras relaciones posibles entre las variables, por ello debemos volver a un punto anterior y descartar como variables predictoras a UNEMP y M. Generamos un nuevo modelo:
model_detroit2 <- lm(H ~ FTP + LIC + GR +
CLEAR + W + NMAN + G + HE + WE,
data = detroit)
summary(model_detroit2)##
## Call:
## lm(formula = H ~ FTP + LIC + GR + CLEAR + W + NMAN + G + HE +
## WE, data = detroit)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 0.31251 -0.80560 0.71173 -0.56435 0.64110 -0.43260 0.78349 0.35282
## 9 10 11 12 13
## -2.08603 0.28411 0.66355 0.07720 0.06208
## attr(,"label")
## [1] "H"
## attr(,"format.spss")
## [1] "F9.2"
## attr(,"display_width")
## [1] 9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.052e+02 1.556e+02 -0.676 0.548
## FTP 8.401e-02 5.382e-02 1.561 0.216
## LIC 1.505e-02 6.992e-03 2.152 0.120
## GR 4.475e-03 4.550e-03 0.984 0.398
## CLEAR -3.810e-01 2.823e-01 -1.349 0.270
## W 1.271e-04 1.537e-04 0.827 0.469
## NMAN 6.842e-02 7.044e-02 0.971 0.403
## G 3.731e-03 2.079e-01 0.018 0.987
## HE 3.255e+00 6.299e+00 0.517 0.641
## WE 4.647e-02 1.596e-01 0.291 0.790
##
## Residual standard error: 1.611 on 3 degrees of freedom
## Multiple R-squared: 0.9976, Adjusted R-squared: 0.9903
## F-statistic: 137.5 on 9 and 3 DF, p-value: 0.0009181Vemos que el R^2 ajustado arroja un valor de \(0.9903\), si bien es un poco menor al
modelo anterior sigue siendo significativo en su nivel de determinación.
Por otro lado, obtemos un valor \(P =
0.0009181\) que resulta considerablemente menor que antes.
Tampoco vemos que los valores P de las variables predictoras nos
indiquen que éstas son estadísticamente significativas, pero decidimos
continuar igual con la selección automática de variables para luego
estudiar los residuos, la multicolinealidad y finalmente tomar alguna
decisión respecto al modelo.
Utilizaremos el método de selección automática Stepwise sobre un modelo que incluya todas las variables, utilizando las probabilidades de entrada y salida del modelo sugeridas en la consigna (Entrada: 0.10 y Salida: 0.15):
(stepwise <- ols_step_both_p(model_detroit, pent = 0.10, prem= 0.15))##
## Stepwise Selection Summary
## --------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## --------------------------------------------------------------------------------------
## 1 CLEAR addition 0.938 0.932 1215.0960 78.4277 4.2643
## 2 LIC addition 0.989 0.987 200.0290 57.3254 1.8393
## 3 HE addition 0.993 0.991 124.7580 53.2520 1.5349
## 4 FTP addition 0.996 0.994 78.5760 49.2182 1.2908
## 5 UNEMP addition 0.997 0.996 49.2760 44.9261 1.0833
## 6 WE addition 0.999 0.999 15.5890 30.8418 0.6303
## 7 HE removal 0.999 0.999 15.1810 30.1883 0.6146
## 8 W addition 1.000 0.999 8.9610 22.9673 0.4656
## --------------------------------------------------------------------------------------Según este método, el mejor modelo debería tener siete variables predictoras. También utilizaremos otro método de selección automática que nos mostrará los mejores modelos para cada cantidad de variables predictoras, esto nos permitirá comparar con el Stepwise.
ols_step_best_subset(model_detroit)| mindex | n | predictors | rsquare | adjr | predrsq | cp | aic | sbic | sbc | msep | fpe | apc | hsp | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6 | 1 | 1 | CLEAR | 0.9379154 | 0.9322714 | 0.9224734 | 1215.095947 | 78.427706 | 35.6200395 | 80.12255 | 236.950785 | 20.9814745 | 0.0846608 | 1.8183945 |
| 40 | 2 | 2 | LIC CLEAR | 0.9894998 | 0.9873997 | 0.9830643 | 200.029013 | 57.325420 | 13.0530640 | 59.58522 | 44.527816 | 4.1636399 | 0.0168004 | 0.3758842 |
| 126 | 3 | 3 | UNEMP LIC WE | 0.9978960 | 0.9971946 | 0.9931895 | 36.484532 | 38.427422 | -4.9009429 | 41.25217 | 10.037808 | 0.9849509 | 0.0039743 | 0.0941497 |
| 390 | 4 | 4 | UNEMP LIC CLEAR WE | 0.9988213 | 0.9982320 | 0.9956383 | 20.239564 | 32.894385 | -8.7923815 | 36.28408 | 6.426472 | 0.6572528 | 0.0026520 | 0.0678118 |
| 846 | 5 | 5 | UNEMP LIC CLEAR W WE | 0.9992298 | 0.9986797 | 0.9956439 | 14.185748 | 29.362900 | -9.2981461 | 33.31755 | 4.899228 | 0.5181001 | 0.0020905 | 0.0590816 |
| 1098 | 6 | 6 | FTP UNEMP LIC CLEAR W WE | 0.9995962 | 0.9991925 | 0.9976599 | 8.960916 | 22.967327 | -5.8647061 | 27.48692 | 3.082022 | 0.3335522 | 0.0013459 | 0.0433618 |
| 1507 | 7 | 7 | FTP UNEMP M LIC CLEAR NMAN WE | 0.9996709 | 0.9992102 | 0.9878065 | 9.488024 | 22.307709 | -0.5404201 | 27.39225 | 3.139751 | 0.3425182 | 0.0013821 | 0.0530088 |
| 1920 | 8 | 8 | FTP M LIC W NMAN G HE WE | 0.9997219 | 0.9991658 | 0.9893626 | 10.482430 | 22.118383 | 6.1473460 | 27.76788 | 3.537484 | 0.3790162 | 0.0015293 | 0.0746547 |
| 2022 | 9 | 9 | FTP M LIC GR W NMAN G HE WE | 0.9998738 | 0.9994950 | 0.9718084 | 9.489177 | 13.853632 | 25.7525589 | 20.06807 | 2.409176 | 0.2398746 | 0.0009679 | 0.0677907 |
| 2041 | 10 | 10 | FTP UNEMP M LIC GR W NMAN G HE WE | 0.9999190 | 0.9995137 | 0.6715221 | 10.597927 | 10.091448 | 28.3490625 | 16.87084 | 3.093141 | 0.2410239 | 0.0009725 | 0.1305546 |
| 2047 | 11 | 11 | FTP UNEMP M LIC GR CLEAR W NMAN G HE WE | 0.9999493 | 0.9993914 | -0.6458061 | 12.000000 | 5.998254 | -30.8941476 | 13.34260 | Inf | 0.3142404 | 0.0012680 | Inf |
Observamos que en ambos métodos de selección, las variables predictores CLEAR y LIC parecerían ser las que aportan mayor significancía. Creamos un objeto y luego un plot para ver de forma más clara el impacto de la incorporación de cada variable al modelo.
(step_sum <- data_frame(
Step = 1:8,
Order = stepwise$orders,
Method = stepwise$method,
Adj_R_2 = stepwise$adjr))| Step | Order | Method | Adj_R_2 |
|---|---|---|---|
| 1 | CLEAR | addition | 0.9322714 |
| 2 | LIC | addition | 0.9873997 |
| 3 | HE | addition | 0.9912251 |
| 4 | FTP | addition | 0.9937939 |
| 5 | UNEMP | addition | 0.9956287 |
| 6 | WE | addition | 0.9985201 |
| 7 | HE | removal | 0.9985931 |
| 8 | W | addition | 0.9991925 |
plot(x = step_sum$Step, y = step_sum$Adj_R_2, type = "o",
xlab = "Modelo de Stepwise",
ylab = "R2 ajustado")
Vemos que si bien el \(R^2\) Ajustado crece, no creemos que aporte valor explicativo de forma significativa.
modelo_0 <- lm(H ~ CLEAR,
data = detroit)
modelo_1 <- lm(H ~ CLEAR + LIC,
data = detroit)
modelo_2 <- lm(H ~ CLEAR + LIC + HE,
data = detroit)
modelo_3 <- lm(H ~ CLEAR + LIC + HE + FTP,
data = detroit)
anova(modelo_0, modelo_1, modelo_2, modelo_3)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11 | 200.02339 | NA | NA | NA | NA |
| 10 | 33.82957 | 1 | 166.193816 | 99.741999 | 0.0000086 |
| 9 | 21.20306 | 1 | 12.626514 | 7.577862 | 0.0249511 |
| 8 | 13.32990 | 1 | 7.873164 | 4.725117 | 0.0614599 |
Realizamos el test de ANOVA para comparar distintos modelos de
regresión que serían competitivos en cuánto a su valor explicativo.
Llegamos a la misma conclusión que observando el \(R^2\) Ajustado. El ANOVA nos indica que
incorporar la variable LIC es estadísticamente significativo, y que al
seguir incorporando variables la significancia estadística va
disminuyendo comparativamente. Considerando lo expuesto, y con la
finalidad de obtener un modelo parsimonioso, nos quedamos con el modelo
H ~ CLEAR + LIC
Debemos considerar que, no siempre hay un mejor modelo; no siempre el mejor modelo tiene más predictores; pueden existir varios modelos que resulten competitivos. Vemos que los modelos son competitivos, entonces para elegir el modelo más adecuado tenemos que tener en cuenta la multicolinealidad. El siguiente paso es analizar la multicolinealidad del modelo, esto coincide con la consigna 7 del ejercicio: Analizar la colinealidad de las variables predictoras presentes en la ecuación
model_detroit <- lm(H ~ CLEAR + LIC, data = detroit)
ols_vif_tol(model_detroit)| Variables | Tolerance | VIF |
|---|---|---|
| CLEAR | 0.6921615 | 1.44475 |
| LIC | 0.6921615 | 1.44475 |
Este fenómeno ocurre cuando las variables predictoras están relacionadas entre sí, lo que implicaría diversos efectos sobre nuestro modelo: Limita el tamaño del \(R^2\) por que los predictores correlacionados “explican” la misma variación; Dificulta la determinación de la importancia de los predictores; Incrementa la varianza de los estimadores de los parámetros.
Vemos que para las variables predictoras no se obtienen valores de VIF > 10 y de Tolerancia < 0.1 que nos indicaría que estamos en presencia de un fenómeno de multicolinealidad.
Consigna 3: Qué información da el coeficiente de determinación?
El coeficiente de determinación es la proporción de la variación en
y que se explica por la línea de regresión. Veamoslo en nuestro
ejercicio, recordamos la formula de trabajo:
lm(H ~ CLEAR + LIC, data = detroit)
summary(model_detroit)##
## Call:
## lm(formula = H ~ CLEAR + LIC, data = detroit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3802 -0.6176 0.6449 1.3765 1.8907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 103.655542 4.820176 21.505 1.05e-09 ***
## CLEAR -1.057474 0.050413 -20.976 1.35e-09 ***
## LIC 0.014136 0.002017 7.009 3.68e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.839 on 10 degrees of freedom
## Multiple R-squared: 0.9895, Adjusted R-squared: 0.9874
## F-statistic: 471.2 on 2 and 10 DF, p-value: 1.276e-10Tenemos un \(R^2=0.9895\) y un valor ajustado de \(R^2_A=0.9874\). El primer valor es el coeficiente de determinación múltiple, nos indica qué tan bien se ajusta la ecuación de regresión múltiple a los datos muestrales. El Coeficiente de determinación ajustado, modifica este valor para tener en cuenta el número de variables y el tamaño de la muestra. Esto resulta relevante ya que el Coeficiente de determinación múltiple \(R^2\) aumenta a medida que se incluyan más variables, incluso aunque éstas sean simple ruido estadístico sin correlación.
En síntesis, con un \(R^2_A=0.9874\) el modelo explica el \(98.74\%\) de la variación de la variable respuesta H.
Consigna 4: Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámeros?
Los modelos de regresión lineal tienen los siguientes requisitos:
La muestra de datos es una muestra aleatoria de datos cuantitativos
Los datos tienen una distribución Normal. Es decir, para cada valor fijo de \(x\) los valores correspondientes de \(y\) siguen una distribución Normal
Linearidad, debe haber correlación lineal entre las variables. Se puede comprobar con un análisis exploratorio de la correlación con el Coeficiente \(R\) y gráficos de dispersión, como se vió al inicio del ejercicio en el análisis exploratorio. También se puede comprobar con un gráfico de residuos.
Homocedasticidad, es decir constancia o igualdad de varianzas. Se puede comprobar con un gráfico de residuos que se verá más adelante.
Independencia, las muestras deben ser independientes. Se puede comprobar mediante el análisis de los residuos.
Ausencia de valores atípicos
Consigna 5: Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera
En un modelo de Regresión Lineal Múltiple, la bondad del ajuste puede aproximarse a través de \(R^2\). En general, un modelo se ajusta bien a los datos a partir de la diferencia entre los valores observados y los valores predichos. Un \(R^2\) elevado, cercano a uno, implica el porcentaje de variabilidad explicado o determinado por el modelo. Esto es así a partir de la misma definición del coeficiente:
- \(SSR/SST\)
SSR se llama suma de cuadrados de la regresión y refleja la cantidad de variación de los valores de y explicados por el modelo.
SST es la suma de cuadrados de y y refleja su variabilidad en torno a su media. Es decir, la variabilidad total
Consigna 6: Realizar un análisis de los residuos del modelo para evaluar el cumplimiento de los supuestos. Para esto, realizar gráficos de los residuos con el valor predicho.
El gráfico de residuos nos permite constatar la homocedasticidad, linealidad y normalidad. Se observará un gráfico de dispersión donde se reemplazan los valores de \(y\) por el residual \(y-\hat{y}\), entonces los puntos de la gráfica son \((x, y-\hat{y})\). Retomamos la ecuación de trabajo y graficamos:
model_detroit$call## lm(formula = H ~ CLEAR + LIC, data = detroit)ols_plot_resid_stud_fit(model_detroit, print_plot = TRUE)
Vemos que sólo uno de los valores se encuentra fuera del rango (2;-2). El problema surge aquí al ver que los Residuos siguen un patrón determinado cuando idealmente no debería seguir ningún patrón obvio para confirmar que los datos muestrales seguirían un patrón de línea recta. Vemos que la gráfica no se ensancha, esto confirmaría la homocedasticidad.
Consigna 7: Analizar la colinealidad de las variables predictoras presentes en la ecuación.
La colinealidad ya ha sido analizada en el punto 2 para la elección del modelo.
Consigna 8: Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso.
Continuamos con nuestro análisis haciendo un gráfico que detecta valores atípicos (outliers) y puntos influyentes (leverage)
ols_plot_resid_lev(model_detroit, print_plot = TRUE)
Detectamos el caso 8 como punto influyente de nuestro modelo que podría alterar la línea de regresión y por lo tanto los resultados. En este sentido, hacemos un gráfico de la Distancia de Cook para establecer la importancia de este leverage.
ols_plot_cooksd_bar(model_detroit, print_plot = TRUE)
Nuevamente, vemos que el caso 2 representa un outlier significativo. Cabe destacar también la poca cantidad de casos de la muestra para mantener este caso y este modelo.
Continuamos con un análisis de los gráficos dfBetas que permitirá observar la influencia de los leverage sobre la estimación de los parámetros para las variables predictoras.
ols_plot_dfbetas(model_detroit, print_plot = TRUE)
Al parecer, el caso 2 tendría una influencia importante es la estimación de los parámetros. Sin embargo, también se sugiere que el caso 12 también tiene una influencia ligeramente importante. Esto nos lleva a pensar que eliminar el caso 2 llevaría a sobreajuste innecesario que luego podría determinar el caso 12 como punto influyente. Por lo tanto, decidimos mantener el modelo actual.
Ecuación
La ecuación de regresión lineal sería:
model_detroit$call## lm(formula = H ~ CLEAR + LIC, data = detroit)model_detroit$coefficients## (Intercept) CLEAR LIC
## 103.65554203 -1.05747371 0.01413625
- \(Y = 103.655\:-\: 1..057*X_C\:+\:0.014*X_L\)
Siendo \(X_C\) = CLEAR y \(X_L\) = LIC
Las variables predictoras de H (Tasa de homicidios cada 100.000 habitantes) en nuestro modelo son CLEAR (Porcentaje de homicidios resueltos con arresto del responsable) y LIC (Cantidad de licencias de portación de armas cada 100.000 habitantes).
Ejercicio 3: Regresión Logística
El bajo peso al nacer, definido por un peso al nacer inferior a 2500 gr., ha sido una preocupación de los médicos durante años debido a que tanto las tasas de mortalidad como la de nacimientos defectuosos son muy altas para los niños con bajo peso al nacer. El comportamiento de la mujer durante el embarazo (incluyendo la dieta, los hábitos tabáquicos y los cuidados prenatales) pueden alterar las chances de un parto de un niño con bajo peso.
Los datos que se presentan en este ejercicio corresponden a 189 nacimientos de los cuales 59 han resultado en niños con bajo peso. El objetivo de este ejercicio es determinar cuáles de las variables presentes en la base de datos que se adjunta son factores de riesgo de bajo peso al nacer.
La base de datos que se presenta (archivo LOWBWT.sav) contiene las siguientes variables:
• LOW: Bajo peso al nacer. (0: \(\geq\) 2500 g; 1: \(<\) 2500 g) (variable dependiente)
• AGE: Edad de la madre
• LWT: Peso de la madre el momento de la última menstruación (en libras)
• RACE: Raza (1 = White; 2 = Black; 3 = Other)
• SMOKE: Fumó durante el embarazo (0 = No 1 = Yes)
• PTL: Antecedentes de embarazos prematuros (0 = None; 1 = One; 2 = Two, etc).
• HT: Antecedentes de hipertensión arterial (0 = No; 1 = Yes)
• UI: Irritabilidad uterina (0 = No; 1 = Yes)
• FTV: Cantidad de consultas obstétricas durante el primer trimestre (0 = None; 1 = One; 2 = Two, etc.)
• BWT: Peso al nacer del bebé en gramos
Se requiere construir una ecuación de regresión logística que relacione la variable dicotómica que indica si se trata de un nacimiento con bajo peso al nacer con el resto de las variables que corresponda, y determinar si estas variables son útiles para predecir la variable dependiente
- Calcular el riesgo relativo y los odds ratio de la variable dependiente con todas las variables dicotómicas. Analizar los resultados.
- Cuál es la definición de odds ratio? Qué información suministra y de qué manera puede calcularse utilizando la regresión logística?
- Calcule los odds ratio de cada una de las variables predictoras con la variable dependiente? Comentar.
- Realizar una regresión logística múltiple, seleccionando los mejores predictores entre las variables independientes disponibles, utilizando un método de selección automática. Describir el proceso de selección automática utilizado. (Sug. Considerar como probabilidad de entrada 0.10 y de salida del modelo 0.15.)
- Según el modelo obtenido, cuáles son los principales factores de riesgo del bajo peso y cuál es la magnitud de su efecto?
- Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámetros?
- Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
- Indicar porcentaje de casos bien predichos por el modelo.
- Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso.
Resolución Ejercicio 3
Cargamos los datos:
LOWBWT <- read_sav("3_REG LOG/7_LOWBWT.sav")
head(LOWBWT)| ID | LOW | AGE | LWT | RACE | SMOKE | PTL | HT | UI | FTV | BWT |
|---|---|---|---|---|---|---|---|---|---|---|
| 4 | 1 | 28 | 120 | 3 | 1 | 1 | 0 | 1 | 0 | 709 |
| 10 | 1 | 29 | 130 | 1 | 0 | 0 | 0 | 1 | 2 | 1021 |
| 11 | 1 | 34 | 187 | 2 | 1 | 0 | 1 | 0 | 0 | 1135 |
| 13 | 1 | 25 | 105 | 3 | 0 | 1 | 1 | 0 | 0 | 1330 |
| 15 | 1 | 25 | 85 | 3 | 0 | 0 | 0 | 1 | 0 | 1474 |
| 16 | 1 | 27 | 150 | 3 | 0 | 0 | 0 | 0 | 0 | 1588 |
Consigna 1: Calcular el riesgo relativo y los odds ratio de la variable dependiente con todas las variables dicotómicas. Analizar los resultados.
Laa variable dependiente es LOW, mientras que las variables dicotómicas para esta consigna son: SMOKE, HT y UI.
smoke_l <- glm(LOW ~ SMOKE,
data = LOWBWT,
family = "binomial")
summary(smoke_l)##
## Call:
## glm(formula = LOW ~ SMOKE, family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0197 -0.7623 -0.7623 1.3438 1.6599
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
## SMOKE 0.7041 0.3196 2.203 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.80 on 187 degrees of freedom
## AIC: 233.8
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = smoke_l,
incr = list(SMOKE = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| SMOKE | 2.022 | 1.082 | 3.801 | 1 |
Esta función nos devuelve el Odds Ratio (OR) de 2.022, y una estimación por intervalos de confianza con un \(L_i=1.082\) y \(L_s=3.081\). El valor de Odds Ratio expresa la posibilidad de ocurrencia de un evento, en este caso la posibilidad de LOW = 1 (bajo peso al nacer) prácticamente se duplica para SMOKE = 1 (Fumar durante el embarazo) respecto a SMOKE = 0 (no fumar durante el embarazo). El IC nos indica que la relación es estadísticamente significativa, dado que el Límite inferior sigue siendo mayor a 1, es decir que la posibilidad se incrementa de la misma forma que vimos recién. Si el OR fuese igual a 1 implicaría que no hay efecto de la variable independiente sobre la dependiente.
Repetimos el proceso para la variable predictora HT: Antecedentes de hipertensión arterial.
ht_l <- glm(LOW ~ HT,
data = LOWBWT,
family = "binomial")
summary(ht_l)##
## Call:
## glm(formula = LOW ~ HT, family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3232 -0.8341 -0.8341 1.5652 1.5652
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8771 0.1650 -5.315 1.07e-07 ***
## HT 1.2135 0.6083 1.995 0.0461 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 230.65 on 187 degrees of freedom
## AIC: 234.65
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = ht_l,
incr = list(HT = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| HT | 3.365 | 1.028 | 11.829 | 1 |
Nuevamente el límite inferior es mayor a 1, por lo que la relación es estadísticamente significativa. Entonces, pareciera que la posibilidad de LOW = 1 es 3.365 veces más grande si HT = 1 respecto a HT=0.
Finalmente, replicamos para la variable UI: Irritabilidad uterina
ui_l <- glm(LOW ~ UI,
data = LOWBWT,
family = "binomial")
summary(ui_l)##
## Call:
## glm(formula = LOW ~ UI, family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1774 -0.8097 -0.8097 1.1774 1.5967
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9469 0.1756 -5.392 6.97e-08 ***
## UI 0.9469 0.4168 2.272 0.0231 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.60 on 187 degrees of freedom
## AIC: 233.6
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = ui_l,
incr = list(UI = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| UI | 2.578 | 1.133 | 5.881 | 1 |
Llegamos a la misma conclusión que para los OR anteriores.
Consigna 2: Cuál es la definición de odds ratio? Qué información suministra y de qué manera puede calcularse utilizando la regresión logística?
Los Odds Ratio se definen como el cociente que surge entre las posibilidades de ocurrencia de un evento, y sus posibilidades de no ocurrencia.
La interpretación de este cociente es que si se obtiene un valor mayor a 1, indica que si aumenta el predictor los odds de la variable dependiente crecen. Si es igual a 1 el efecto sobre la dependiente sería nulo. Un valor menor a 1 indicaría que la variable dependiente decrece. Por este motivo, el Odds Ratio es una buena medida del tamaño del efecto.
Consigna 3: Calcule los odds ratio de cada una de las variables predictoras con la variable dependiente? Comentar.
Las tres variables dicotómicas ya fueron analizadas en la Consigna 1, ahora nos quedaría calcular los OR para las otras variables cualitativas no dicotómicas (RACE + PTL + FTV) y para las variables cuantitativas (AGE + LWT).
race_l <- glm(LOW ~ factor(RACE),
data = LOWBWT,
family = "binomial")
summary(race_l)##
## Call:
## glm(formula = LOW ~ factor(RACE), family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0489 -0.9665 -0.7401 1.4042 1.6905
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1550 0.2391 -4.830 1.36e-06 ***
## factor(RACE)2 0.8448 0.4634 1.823 0.0683 .
## factor(RACE)3 0.6362 0.3478 1.829 0.0674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = race_l,
incr = list(RACE = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| factor(RACE)2 | 2.328 | 0.926 | 5.775 | Indicator variable |
| factor(RACE)3 | 1.889 | 0.957 | 3.758 | Indicator variable |
Para esta variable vemos un OR mayor a 1 en las categorías, pero un límite inferior menor a 1 que indicaría que no es estadísticamente significativo.
ptl_l <- glm(LOW ~ factor(PTL),
data = LOWBWT,
family = "binomial")
summary(ptl_l)##
## Call:
## glm(formula = LOW ~ factor(PTL), family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4823 -0.7723 -0.7723 0.9005 1.6464
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0571 0.1813 -5.831 5.5e-09 ***
## factor(PTL)1 1.7503 0.4694 3.728 0.000193 ***
## factor(PTL)2 0.6516 0.9307 0.700 0.483821
## factor(PTL)3 -13.5090 882.7434 -0.015 0.987790
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 218.80 on 185 degrees of freedom
## AIC: 226.8
##
## Number of Fisher Scoring iterations: 13or_glm(data = LOWBWT,
model = ptl_l,
incr = list(PTL = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| factor(PTL)1 | 5.756 | 2.353 | 1.514200e+01 | Indicator variable |
| factor(PTL)2 | 1.919 | 0.246 | 1.197000e+01 | Indicator variable |
| factor(PTL)3 | 0.000 | NA | 1.813066e+72 | Indicator variable |
Para la variable PTL (Antecedentes de embarazos prematuros) vemos que resulta significativa para sólo un valor de la variable.
ftv_l <- glm(LOW ~ factor(FTV),
data = LOWBWT,
family = "binomial")
summary(ftv_l)##
## Call:
## glm(formula = LOW ~ factor(FTV), family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3018 -0.9448 -0.7302 1.4294 1.7060
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5754 0.2083 -2.762 0.00575 **
## factor(FTV)1 -0.6103 0.4026 -1.516 0.12957
## factor(FTV)2 -0.6142 0.4793 -1.281 0.20003
## factor(FTV)3 0.8630 0.7917 1.090 0.27564
## factor(FTV)4 -0.5232 1.1733 -0.446 0.65564
## factor(FTV)6 -13.9907 882.7434 -0.016 0.98735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 228.49 on 183 degrees of freedom
## AIC: 240.49
##
## Number of Fisher Scoring iterations: 13or_glm(data = LOWBWT,
model = ftv_l,
incr = list(FTV = 1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| factor(FTV)1 | 0.543 | 0.239 | 1.16900e+00 | Indicator variable |
| factor(FTV)2 | 0.541 | 0.198 | 1.33100e+00 | Indicator variable |
| factor(FTV)3 | 2.370 | 0.496 | 1.25950e+01 | Indicator variable |
| factor(FTV)4 | 0.593 | 0.029 | 4.82200e+00 | Indicator variable |
| factor(FTV)6 | 0.000 | NA | 1.11173e+72 | Indicator variable |
Para la variable FTV (Cantidad de consultas obstétricas durante el primer trimestre) no parecía haber relación significativa con la variable dependiente.
A continuación, seguiremos con el mismo análisis para las variables cuantitativas
age_l <- glm(LOW~AGE,data=LOWBWT,family="binomial")
summary(age_l)##
## Call:
## glm(formula = LOW ~ AGE, family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0402 -0.9018 -0.7754 1.4119 1.7800
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.38458 0.73212 0.525 0.599
## AGE -0.05115 0.03151 -1.623 0.105
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 231.91 on 187 degrees of freedom
## AIC: 235.91
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = age_l,
incr = list(AGE=5))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| AGE | 0.774 | 0.562 | 1.046 | 5 |
lwt_l <- glm(LOW~LWT,data=LOWBWT,family="binomial")
summary(lwt_l)##
## Call:
## glm(formula = LOW ~ LWT, family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0951 -0.9022 -0.8018 1.3609 1.9821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99831 0.78529 1.271 0.2036
## LWT -0.01406 0.00617 -2.279 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 228.69 on 187 degrees of freedom
## AIC: 232.69
##
## Number of Fisher Scoring iterations: 4or_glm(data = LOWBWT,
model = lwt_l,
incr = list(LWT=15))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| LWT | 0.81 | 0.667 | 0.961 | 15 |
En ninguno de los dos casos el OR pareciera ser significativo.
En definitiva, las variables independientes que mostraron ser significativas en este primer análisis son:
- SMOKE
- HT
- UI
- PTL = 1
Consigna 4: Realizar una regresión logística múltiple, seleccionando los mejores predictores entre las variables independientes disponibles, utilizando un método de selección automática. Describir el proceso de selección automática utilizado. (Sug. Considerar como probabilidad de entrada 0.10 y de salida del modelo 0.15.)
Para comenzar con el proceso, vamos a correr un modelo que incluya todas las variables independientes:
model_low <- glm(LOW ~ AGE + LWT +
factor(RACE) + factor(PTL) + factor(FTV)+
SMOKE + HT + UI,
data = LOWBWT,
family = binomial)
or_glm(data = LOWBWT,
model = model_low,
incr = list(AGE=5,LWT=15,
RACE=1, PTL=1,FTV=1,
SMOKE=1, HT=1, UI=1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| AGE | 0.823 | 0.553 | 1.205000e+00 | 5 |
| LWT | 0.794 | 0.631 | 9.820000e-01 | 15 |
| factor(RACE)2 | 3.036 | 1.033 | 9.114000e+00 | Indicator variable |
| factor(RACE)3 | 1.954 | 0.772 | 5.077000e+00 | Indicator variable |
| factor(PTL)1 | 6.456 | 2.190 | 2.118800e+01 | Indicator variable |
| factor(PTL)2 | 1.628 | 0.184 | 1.166300e+01 | Indicator variable |
| factor(PTL)3 | 0.000 | NA | 8.053712e+121 | Indicator variable |
| factor(FTV)1 | 0.571 | 0.207 | 1.470000e+00 | Indicator variable |
| factor(FTV)2 | 0.922 | 0.301 | 2.610000e+00 | Indicator variable |
| factor(FTV)3 | 3.013 | 0.549 | 1.823500e+01 | Indicator variable |
| factor(FTV)4 | 0.398 | 0.015 | 5.032000e+00 | Indicator variable |
| factor(FTV)6 | 0.000 | NA | 3.051192e+122 | Indicator variable |
| SMOKE | 2.029 | 0.859 | 4.877000e+00 | 1 |
| HT | 6.078 | 1.453 | 2.875600e+01 | 1 |
| UI | 2.210 | 0.860 | 5.684000e+00 | 1 |
Al realizar el modelo con todas las variables, sólo algunas variables se muestran significativas en cuanto al valor del OR:
- RACE = 2
- PTL = 1
- HT
Utilizamos la función summary() para visualizar los
estimadores y el valor P de las variables predictoras:
summary(model_low)##
## Call:
## glm(formula = LOW ~ AGE + LWT + factor(RACE) + factor(PTL) +
## factor(FTV) + SMOKE + HT + UI, family = binomial, data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7713 -0.7839 -0.4859 0.8305 2.1983
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.986e-01 1.299e+00 0.692 0.48908
## AGE -3.889e-02 3.952e-02 -0.984 0.32508
## LWT -1.536e-02 7.483e-03 -2.052 0.04017 *
## factor(RACE)2 1.110e+00 5.509e-01 2.016 0.04382 *
## factor(RACE)3 6.698e-01 4.774e-01 1.403 0.16061
## factor(PTL)1 1.865e+00 5.722e-01 3.259 0.00112 **
## factor(PTL)2 4.875e-01 1.008e+00 0.484 0.62858
## factor(PTL)3 -1.553e+01 1.455e+03 -0.011 0.99148
## factor(FTV)1 -5.598e-01 4.958e-01 -1.129 0.25879
## factor(FTV)2 -8.141e-02 5.447e-01 -0.149 0.88120
## factor(FTV)3 1.103e+00 8.648e-01 1.276 0.20213
## factor(FTV)4 -9.225e-01 1.384e+00 -0.667 0.50496
## factor(FTV)6 -1.291e+01 1.455e+03 -0.009 0.99292
## SMOKE 7.073e-01 4.399e-01 1.608 0.10784
## HT 1.805e+00 7.488e-01 2.410 0.01595 *
## UI 7.930e-01 4.780e-01 1.659 0.09712 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 188.60 on 173 degrees of freedom
## AIC: 220.6
##
## Number of Fisher Scoring iterations: 14Vemos que habría una variable más que resulta estadísticamente significativa. El listado completo quedaría:
- RACE = 2
- PTL = 1
- HT
- LWT
Más allá de este análisis, sabemos que correr un modelo con todas las variables independientes puede mostrarnos relaciones significativas que no lo son y ocultar otras que sí lo son. Por ello, utilizaremos una herramienta de selección automática Stepwise, similar al utilizado en la Regresión Lineal pero basado en el AIC y en la Deviance para realizar dicha selección.
## Selección automática de modelos
model_sel <- glm (LOW ~ 1,
family = "binomial",
data =LOWBWT)
step_sel <- stepAIC(model_sel, scope = list(
upper = ~AGE + LWT +
factor(RACE) + factor(PTL) + factor(FTV)+
SMOKE + HT + UI),
direction = c("both"),
trace = F)
step_sel$anova| Step | Df | Deviance | Resid. Df | Resid. Dev | AIC |
|---|---|---|---|---|---|
| NA | NA | 188 | 234.6720 | 236.6720 | |
| + factor(PTL) | 3 | 15.872143 | 185 | 218.7999 | 226.7999 |
| + LWT | 1 | 5.264712 | 184 | 213.5351 | 223.5351 |
| + HT | 1 | 7.293221 | 183 | 206.2419 | 218.2419 |
| + UI | 1 | 3.646880 | 182 | 202.5950 | 216.5950 |
| + AGE | 1 | 2.497925 | 181 | 200.0971 | 216.0971 |
| + SMOKE | 1 | 2.190360 | 180 | 197.9068 | 215.9068 |
| + factor(RACE) | 2 | 5.452610 | 178 | 192.4541 | 214.4541 |
| - AGE | 1 | 1.136371 | 179 | 193.5905 | 213.5905 |
El modelo final quedaría determinado por el modelo con la siguiente formula:
final_model <- glm(LOW ~ factor(PTL) + LWT +
HT + UI + SMOKE + factor(RACE),
family = "binomial",
data = LOWBWT)
final_model$call## glm(formula = LOW ~ factor(PTL) + LWT + HT + UI + SMOKE + factor(RACE),
## family = "binomial", data = LOWBWT)El modelo seleccionado coincide con lo expuesto sobre los OR en el primer modelo múltiple con todas las variables independientes, sumando las variables SMOKE y UI que también se mostraron significativas en el análisis de regresión simple.
Consigna 5: Según el modelo obtenido, cuáles son los principales factores de riesgo del bajo peso y cuál es la magnitud de su efecto?
En primera instancia, para responder a esta consigna debemos analizar los OR de nuestro modelo:
or_glm(data = LOWBWT,
model = final_model,
incr = list(LWT=15,RACE=1, PTL=1,
SMOKE=1, HT=1, UI=1))| predictor | oddsratio | ci_low (2.5) | ci_high (97.5) | increment |
|---|---|---|---|---|
| factor(PTL)1 | 4.297 | 1.625 | 1.212800e+01 | Indicator variable |
| factor(PTL)2 | 1.315 | 0.155 | 8.937000e+00 | Indicator variable |
| factor(PTL)3 | 0.000 | NA | 4.446159e+71 | Indicator variable |
| LWT | 0.773 | 0.619 | 9.420000e-01 | 15 |
| HT | 6.674 | 1.693 | 2.981600e+01 | 1 |
| UI | 2.445 | 0.968 | 6.192000e+00 | 1 |
| SMOKE | 2.424 | 1.098 | 5.521000e+00 | 1 |
| factor(RACE)2 | 3.486 | 1.222 | 1.013100e+01 | Indicator variable |
| factor(RACE)3 | 2.218 | 0.933 | 5.452000e+00 | Indicator variable |
En cuanto a los OR, los resultados son similares a los obtenidos previamente. La variable con mayor efecto pareciera ser HT = 1, es decir que los antecedentes de Hipertensión arterial implicarían una posibilidad 6.674 veces más alta de un Bajo peso al nacer. Asímismo, para PTL = 1 (1-un- Antecedente de embarazo prematuro) se presenta un OR de 4.297 sobre la posibilidad de LOW = 1.
También resultan útiles los valores P del modelo para validar lo expuesto:
summary(final_model)##
## Call:
## glm(formula = LOW ~ factor(PTL) + LWT + HT + UI + SMOKE + factor(RACE),
## family = "binomial", data = LOWBWT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8644 -0.7707 -0.5171 0.9271 2.2084
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.030369 0.986202 0.031 0.97543
## factor(PTL)1 1.457868 0.507406 2.873 0.00406 **
## factor(PTL)2 0.273850 0.980762 0.279 0.78007
## factor(PTL)3 -14.744564 882.743533 -0.017 0.98667
## LWT -0.017173 0.007121 -2.412 0.01588 *
## HT 1.898206 0.717535 2.645 0.00816 **
## UI 0.894205 0.469649 1.904 0.05691 .
## SMOKE 0.885373 0.409389 2.163 0.03057 *
## factor(RACE)2 1.248872 0.535197 2.333 0.01962 *
## factor(RACE)3 0.796707 0.447359 1.781 0.07493 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 193.59 on 179 degrees of freedom
## AIC: 213.59
##
## Number of Fisher Scoring iterations: 13Vemos que las variables con mayor OR, analizados recientemente, también tienen un menor valor de P.
Consigna 6: Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámetros?
La Regresión Logística comparte algunos supuestos con la Regresión Lineal:
- Linealidad: esto se logra ya que se aplica un logarítmo natural a la variable respuesta, y se espera una relación lineal con las variables predictoras.
- Distribución: se trataría de una distribución de probabilidad Binomial, ya que la variable respuesta de dicotómica
Consigna 7: Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
En principio, recordamos que en Regresión Logística predecir la Probabilidad de Y, en nuestro caso \(P(LOW)\); en realidad, la medida que se utiliza es el Odds Ratio. Además, para realizar la regresión se utiliza la Función de Máxima Verosimilitud (a diferencia del Método de Mínimos Cuadrados de la Regresión Lineal).
Ahora bien, el logarítimo de máxima verosimilitud del modelo se
multiplica por \(-2\) y se obtiene la
deviance que se representa como \(D =
-2LL\). La devianza tiene una distribución de probabilidad
conocida chi-cuadrado y compara los valores predichos con
los observados. El valor de la deviance tendería a ser más bajo
si el modelo es bueno.
Para conocer la eficacia del modelo para predecir la variable de respuesta, utilizamos el estadístico de prueba chi-cuadrado.
dev <- final_model$deviance # deviance del modelo
nullDev <- final_model$null.deviance # deviance nula
modelChi <- nullDev - dev # diferencia
chigl <- final_model$df.null - final_model$df.residual # grados de libertad
chisq.prob <- 1 - pchisq(modelChi, chigl)
format(chisq.prob, scientific = FALSE)## [1] "0.000004834"La probabilidad es menor que nuestro \(\alpha=0.05\) rechazamos la hipótesis nula y concluimos que el modelo es estadísticamente significativo para predecir la variable respuesta.
Consigna 8: Indicar porcentaje de casos bien predichos por el modelo.
Es posible clasificar a los individuos en función de los estimadores del modelo. Como es posible calcular la probabilidad de cada sujeto, se puede “predecir” un valor de Y en términos de esa probabilidad, y calcular los aciertos. Para esto, es necesario decidir a partir de qué probabilidad se considera al individuo como que tuvo el evento
predic_final <- plogis(predict(final_model, LOWBWT))
optCutOff_final <- optimalCutoff(LOWBWT$LOW, predic_final)
optCutOff_final## [1] 0.4352869Se define el punto óptimo de corte en 0.4352869; es decir que si la probabilidad estimada para cada caso es mayor que el Corte Óptimo se clasificará como LOW = 1.
En este sentido, se presentan dos indicadores de Clasificación: sensitividad y especificidad
sensitivity(LOWBWT$LOW, predic_final, threshold = optCutOff_final)## [1] 0.5254237La Sensitividad mide la proporción de casos positivos correctamente clasificados como tales. Para nuestro modelo se estima que la proporción de positivos verdaderos es 0.5254237, por lo que para la sensitividad se presenta un Error de Clasificación de 1 - 0.5254237 = 0.4745763
Este error de clasificación implicaría un \(47.45\%\) de Falsos negativos.
specificity(LOWBWT$LOW, predic_final, threshold = optCutOff_final)## [1] 0.8769231La Especificidad mide la proporción de casos negativos correctamente clasificados como tales. Para nuestro modelo se estima un Error de Clasificación de 1 - 0.8769231 = 0.1230769
Este error de clasificación implicaría un \(12.3\%\) de Falsos positivos.
Esto también se puede formular a través de una Matriz de Confusión
confusionMatrix(LOWBWT$LOW, predic_final, threshold = optCutOff_final)| 0 | 1 | |
|---|---|---|
| 0 | 114 | 28 |
| 1 | 16 | 31 |
Finalmente, la Curva de ROC nos permite visualizar el porcentaje bien predicho por el modelo al cruzar la Especificidad y la Sensitividad, y representar dicho porcentaje en el área bajo la curva
plotROC(LOWBWT$LOW, predic_final)
Mientras mayor es el área cubierta, mejor es el carácter predictor del modelo. También se establecen una serie de niveles para concluir qué tan bueno es el test:
- (0.5): Es como lanzar una moneda.
- (0.5; 0.6): Test malo.
- (0.6; 0.75): Test regular.
- (0.75; 0.9): Test bueno.
- (0.9; 0.97): Test muy bueno.
- (0.97; 1): Test excelente.
Con estos criterios de calidad, se concluye que el modelo es bueno y clasifica correctamente un \(76.98\%\) de los casos.