Modelo Lineal

“El análisis de la varianza, la regresión lineal y la regresión logística tiene la misma lógica matemática pero resuelven problemas distintos y son todos modelos lineales” (Apuntes de clase).

Veremos a continuación una aplicación de los distintos modelos lineales mencionados anteriores en cada uno de los ejercicios.

Consigna del Ejercicio 1

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.

resultados

1) ¿Existen diferencias estadísticamente significativas debidas al tipo de harina utilizada?

Para resolver este ejercicio necesitamos comparar medias entre grupos, es decir, vamos a realizar pruebas de hipótesis. Lo primero que debemos hacer es preparar nuestra base de datos para trabajarla en R.

Al realizar las pruebas de hipótesis con R, obtendremos el p-value, que es el valor que nos permitirá decidir si aceptamos o rechazamos la Hipotesis Nula (Ho), en relación al alfa, que en definitiva es lo que lxs investigadores/as podemos definir o controlar. El p-value nos indica la probabilidad de equivocarnos al rechazar Ho.

Estadistica_Mario_Triola

La Hipótesis Nula (Ho) enuncia que no hay diferencias o cambios entre los grupos que vamos a comparar. La Hipótesis alternativa que si hay diferencias entre los grupos.

Para este ejercicio nuestro alfa será de 0.05

Ahora que hemos aclarado el método que necesitamos utilizar, procedemos a cargar las librerías que vamos a requerir para el procesamiento de los datos en R:

library(tidyverse)
library(ggpubr)
library(rstatix)
library(plyr)
library(dplyr)
library(apa)
library(descr)
library(openxlsx)
library(haven)
library(ggplot2)
library(olsrr)
library(oddsratio)
library(MASS)
library(InformationValue)
library(plotROC)

Podemos leer con R una base de datos en distintos formatos o generarla escribiendo los vectores y uniéndoles en un data frame o base de datos. En este caso procedemos a leer con R la base de datos que construimos en formato excel y vemos los primeros 6 datos de la base:

base_torta <- read.xlsx("C:/Users/Lorena/Documents/UNTREF/5.ESTADISTICA 2/2021/tp2021/anova/base/base_torta.xlsx")
head(base_torta)
##     Torta Azucar Densidad      Harina
## 1 Primera      0     0.90     General
## 2 Primera     50     0.86     General
## 3 Primera     75     0.93     General
## 4 Primera    100     0.79     General
## 5 Primera      0     0.91 Para tortas
## 6 Primera     50     0.88 Para tortas

La base tiene la misma estructura que en formato excel.

bd

A continuación vamos a analizar que tipo de variables tiene la base de datos en R (numérica, factor o character). Ya que para que el ANOVA funcione correctamente debe relacionar una variable categórica o cualitativa (factor) con otra cuantitativa.

str(base_torta)
## 'data.frame':    24 obs. of  4 variables:
##  $ Torta   : chr  "Primera" "Primera" "Primera" "Primera" ...
##  $ Azucar  : num  0 50 75 100 0 50 75 100 0 50 ...
##  $ Densidad: num  0.9 0.86 0.93 0.79 0.91 0.88 0.86 0.86 0.87 0.89 ...
##  $ Harina  : chr  "General" "General" "General" "General" ...

Creamos otras variables dentro de la base que contengan la misma información pero que este definidas como factor para R, esto lo hacemos de la siguiente manera:

base_torta <- base_torta%>%
  mutate(fharina=factor(Harina, labels=(c("General","Para tortas"))))

base_torta <- base_torta%>%
  mutate(fazucar=factor(Azucar, labels=(c("0","50","75","100"))))

Con la función str comprobamos que haya creado las variables factor con 4 niveles para azúcar y 2 niveles para harina:

str(base_torta)
## 'data.frame':    24 obs. of  6 variables:
##  $ Torta   : chr  "Primera" "Primera" "Primera" "Primera" ...
##  $ Azucar  : num  0 50 75 100 0 50 75 100 0 50 ...
##  $ Densidad: num  0.9 0.86 0.93 0.79 0.91 0.88 0.86 0.86 0.87 0.89 ...
##  $ Harina  : chr  "General" "General" "General" "General" ...
##  $ fharina : Factor w/ 2 levels "General","Para tortas": 1 1 1 1 2 2 2 2 1 1 ...
##  $ fazucar : Factor w/ 4 levels "0","50","75",..: 1 2 3 4 1 2 3 4 1 2 ...

Podemos corrobar que ahora tenemos las variables categóricas como factor. Son las variables que utilizaremos para comparar los grupos dada la consigna del ejercicio, y dado que nuestra variable cuantitativa es “densidad”.

En primer lugar, realizaremos un análisis descriptivo de la variable densidad por tipo de harina:

base_torta %>%
  group_by(fharina) %>%
  get_summary_stats(Densidad, type = "mean_sd")
## # A tibble: 2 x 5
##   fharina     variable     n  mean    sd
##   <fct>       <chr>    <dbl> <dbl> <dbl>
## 1 General     Densidad    12 0.868 0.044
## 2 Para tortas Densidad    12 0.851 0.035

Para los dos tipo de harina la media y el desvio standar tiene valores muy similares, por lo que hasta aquí no se observan diferencias en la densidad según el tipo de Harina.

Podemos observalo gráficamente a través de:

ggboxplot(base_torta, x = "fharina", y = "Densidad")

Al parecer tiene una variabilidad similar aunque el tipo de harina “general” parece obtener niveles más altos de densidad en relación a la harina “para tortas”.

Para realizar el test ANOVA y comparar la densidad según los grupos de tipo de harina, necesitamos chequear previamente que se cumplan tres supuestos:

  1. que los datos tenga una distribución normal

  2. que las varianzas sean iguales, y

  3. que las muestras sean independientes.

En este caso, damos por descontado que estamos trabajando con muestras independientes y procedemos a verificar los supuestos a través de pruebas de hipotesis.

1)Test de normalidad:

A continuación realizamos el test de normalidad (Kolmogorov-Smirnov para observaciones mayores a 50 y el test de Shapiro para observaciones menores a 50) en cada grupo de tipo de harina.

Ho: los datos corresponden a una distribución normal

base_torta %>%
  group_by(fharina) %>%
  shapiro_test(Densidad)
## # A tibble: 2 x 4
##   fharina     variable statistic     p
##   <fct>       <chr>        <dbl> <dbl>
## 1 General     Densidad     0.924 0.319
## 2 Para tortas Densidad     0.950 0.643

Observamos el gráfico de normalidad con:

ggqqplot(base_torta, "Densidad", facet.by = "fharina")

Podemos observar que los grupos tienen una distribución normal como se visualiza en el gráfico, y en los dos p-value del test de shapiro resulta mayor que 0,05, por lo que NO rechazamos nuestra hipótesis nula.

En consecuencia, podemos decir que con 0.05 nivel de significancia se observa evidencia estadística de que la variable tipo de harina tiene una distribución normal.

2)Test de Homocedasticidad

Realizamos el test de homocedasticidad de igualdad de varianzas Ho: No hay diferencias entre las varianzas de los grupos

base_torta %>% levene_test(Densidad ~ fharina)
## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    22     0.528 0.475

En este sentido, dado que el p-value del test de levene es de 0.475 podemos decir que existe evidencia estadística de que las varianzas entre los grupos son iguales.

Verificados los supuestos, realizamos el análisis ANOVA Ho: no hay diferencias entre las medias de los grupos por tipo de harina

summary(aov(Densidad~fharina, data=base_torta))
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## fharina      1 0.00184 0.001837    1.16  0.293
## Residuals   22 0.03486 0.001584

Dado que el p-value es de 0.293 no rechazamos Ho.En consecuencia podemos afirmar que con una significancia de alfa de 0.05 no hay diferencias significativas entre las medias de los grupos por tipo de harina utilizado.

2) ¿Existen diferencias estadísticamente significativas debidas al nivel de endulzamiento utilizado?

En este caso, aplicamos el mismo procedimiento utilizado en el punto 1 pero con la variable nivel de azúcar

Análisis descriptivo de densidad de las tortas según nivel de azúcar

base_torta %>%
  group_by(fazucar) %>%
  get_summary_stats(Densidad, type = "mean_sd")
## # A tibble: 4 x 5
##   fazucar variable     n  mean    sd
##   <fct>   <chr>    <dbl> <dbl> <dbl>
## 1 0       Densidad     6 0.88  0.041
## 2 50      Densidad     6 0.865 0.035
## 3 75      Densidad     6 0.865 0.042
## 4 100     Densidad     6 0.828 0.029
ggboxplot(base_torta, x = "fazucar", y = "Densidad")

Observamos que los valores de las medias para los distintos grupos son muy parecidos y que se observan algunas diferencias en la desviación standar. En el gráfico podemos observar que la variabilidad de los datos en el grupo de 0 y de 75 tiene incluso valores atipicos para una distribución normal. Procedemos a chequear los supuestos para aplicar ANOVA:

1)Test de Normalidad para nivel de azúcar

base_torta %>%
  group_by(fazucar) %>%
  shapiro_test(Densidad)
## # A tibble: 4 x 4
##   fazucar variable statistic      p
##   <fct>   <chr>        <dbl>  <dbl>
## 1 0       Densidad     0.728 0.0122
## 2 50      Densidad     0.951 0.748 
## 3 75      Densidad     0.964 0.851 
## 4 100     Densidad     0.889 0.310
ggqqplot(base_torta, "Densidad", facet.by = "fazucar")

Nuevamene, el p-value es mayor en la mayoría de los grupos de nivel de azúcar, excepto para el nivel 0 de azúcar. Por otro lado, se visualiza en el gráfico, valores por fuera del área sombreada o de lo esperable para una distribución normal. Vamos a suponer que es un caso atípico y que no se trata de un error en la carga de datos. Este supuesto no se cumple completamente para todos los grupos de la variable nivel de azúcar.

  1. Test de homocedasticidad o de igualdad de varianzas entre los grupos:
base_torta %>% levene_test(Densidad ~ fazucar)
## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    20    0.0530 0.983

El p-value es 0.983 muy superior a nuestro alfa de 0.05. por lo que No rechazamos Ho.

Apesar de que no se cumple totalmente el test de normalidad, realizamos de todas formas el test ANOVA para nivel de azúcar

summary(aov(Densidad~fazucar, data=base_torta))
##             Df   Sum Sq  Mean Sq F value Pr(>F)
## fazucar      3 0.008713 0.002904   2.076  0.136
## Residuals   20 0.027983 0.001399

Como conclusión en relación a las pruebas realizadas, podemos afirmar que dado el p-value de 0.136 no existe evidencia estadistica para rechazar Ho. Siendo Ho no hay diferencias entre las medias de los grupos por nivel de azúcar.

3) ¿Existe interacción estadísticamente significativa entre los factores considerados?

En este caso, se aplica el ANOVA de dos factores, de la siguiente manera:

ANOVA DE DOS FACTORES SIN INTERACCIÓN, es decir, agregando variables que no tienen interacción entre sí:

Ho: No existen diferencias en la densidad sin interacción entre el nivel de azúcar y el tipo de harina utilizado

dd0 <- aov(Densidad ~ fharina+fazucar, data=base_torta)
anova_apa(dd0)
##        Effect                                                
## 1 (Intercept) F(1, 19) = 12886.60, p < .001, petasq > .99 ***
## 2     fharina F(1, 19) =     1.34, p = .262, petasq = .07    
## 3     fazucar F(3, 19) =     2.11, p = .133, petasq = .25
summary(dd0)
##             Df   Sum Sq  Mean Sq F value Pr(>F)
## fharina      1 0.001838 0.001837   1.335  0.262
## fazucar      3 0.008713 0.002904   2.110  0.133
## Residuals   19 0.026146 0.001376

Para ambos p-value de cada variable no da un valor superior a 0.05 por lo que No rechazamos Ho, pero nos está mostrando que si es significativo el intercept, por ello realizamos el siguiente análisis:

ANOVA DE DOS FACTORES CON INTERACCIÓN, es decir, agregando variables que tienen interacción entre sí:

Ho: No existen diferencias en la densidad con interacción entre el nivel de azúcar y el tipo de harina utilizado

dd <- aov(Densidad ~ fazucar*fharina, data=base_torta)
anova_apa(dd)
##            Effect                                                
## 1     (Intercept) F(1, 16) = 17733.20, p < .001, petasq > .99 ***
## 2         fazucar F(3, 16) =     2.90, p = .067, petasq = .35 .  
## 3         fharina F(1, 16) =     1.84, p = .194, petasq = .10    
## 4 fazucar:fharina F(3, 16) =     3.38, p = .044, petasq = .39 *
summary(dd)
##                 Df   Sum Sq  Mean Sq F value Pr(>F)  
## fazucar          3 0.008713 0.002904   2.904 0.0670 .
## fharina          1 0.001838 0.001838   1.838 0.1941  
## fazucar:fharina  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 ' ' 1

Los p-value obtenidos aqui nos vuelve a indicar que las variables por separado no son significativas en la densidad de las tortas pero si son significativas en su interacción y el tamaño del efecto de la interacción es del 0.39. Lo que indica que el tamaño del efecto es fuerte.

4) Realizar los gráficos de medias e interpretar

with(base_torta, interaction.plot(fazucar, fharina, Densidad, fun = mean,
                                  main = "Interaction Plot"))

En el gráfico observamos claramente la interacción entre las variables y las diferencias en la densidad de las tortas debido a la influencia de los ditintos tipos de harina y los distintos niveles de azúcar.

Como conclusión podemos decir que se visualiza un comportamiento de la variable tipo de torta que cambia radicalmente cuando se trata de niveles de azúcar superior al 75%, que mientras que para las harinas comunes baja considerablemente la densidad para las harinas para torta sube. En consecuencia la tendencia a disminuir o aumentar la densidad de los tipos de harina se ve condicionado por los niveles de azúcar.

Por otro lado, es llamativo observar que para el tipo de harina para tortas tiene una tendencia a bajar la densidad mientras aumenta el nivel azúcar, cambia la tendencia y sube cuando el nivel de azucar es del 100%.

Consigna del Ejercicio 2

Tasa de homicidios en Detroit

En un estudio que investiga el papel de las armas de fuego en el aumento de la tasa de homocidios de Detroit, se recogieron datos para los años 1961-1973. La variable de respuesta (la tasa de homicidios) y las variales potencialmente predictoras del comportamiento de ésta, se describen a continuación:

detroit

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.

  1. Cuál es el nivel de asociación lineal de las variables predictoras con la variable H? Comentar.
  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.)
  3. Qué información da el coeficiente de determinación?
  4. Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámeros?
  5. Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
  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.
  7. Analizar la colinealidad de las variables predictoras presentes en la ecuación.
  8. Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso

Como en esta consigna nos pide encontrar variables dentro de la base de datos disponible que puedan predecir la tasa de homicidios cada 100.000 habitantes. Para ello, debemos aplicar el método de regresión lineal, vamos a empezar a analizar la relación que tiene cada variable potencialmente predictora de la variable H.

detroit <- read_sav("p301.sav")
head(detroit)
## # A tibble: 6 x 13
##    YEAR   FTP UNEMP     M   LIC    GR CLEAR      W  NMAN     G    HE    WE     H
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1  1961  260. 11     456.  178.  216.  93.4 558724  538.  134.  2.98  117.  8.60
## 2  1962  270.  7     480.  156.  180.  88.5 538584  548.  138.  3.09  134.  8.90
## 3  1963  272.  5.20  506.  198.  210.  94.4 519171  563.  144.  3.23  142.  8.52
## 4  1964  273.  4.30  536.  222.  232.  92   500457  591   150.  3.33  148.  8.89
## 5  1965  273.  3.5   576   302.  298.  91   482418  626.  164.  3.46  160. 13.1 
## 6  1966  261.  3.20  602.  391.  368.  87.4 465029  660.  180.  3.60  157. 14.6

Vemos su composición y tipo de datos:

str(detroit)
## tibble [13 x 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"
##   ..- attr(*, "display_width")= int 0
##  $ FTP  : num [1:13] 260 270 272 273 273 ...
##   ..- attr(*, "label")= chr "FTP"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ 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 0
##  $ M    : num [1:13] 456 480 506 536 576 ...
##   ..- attr(*, "label")= chr "M"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ LIC  : num [1:13] 178 156 198 222 302 ...
##   ..- attr(*, "label")= chr "LIC"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ GR   : num [1:13] 216 180 210 232 298 ...
##   ..- attr(*, "label")= chr "GR"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ 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 0
##  $ W    : num [1:13] 558724 538584 519171 500457 482418 ...
##   ..- attr(*, "label")= chr "W"
##   ..- attr(*, "format.spss")= chr "F12.0"
##   ..- attr(*, "display_width")= int 0
##  $ NMAN : num [1:13] 538 548 563 591 626 ...
##   ..- attr(*, "label")= chr "NMAN"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ G    : num [1:13] 134 138 144 150 164 ...
##   ..- attr(*, "label")= chr "G"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ 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 0
##  $ WE   : num [1:13] 117 134 142 148 160 ...
##   ..- attr(*, "label")= chr "WE"
##   ..- attr(*, "format.spss")= chr "F9.2"
##   ..- attr(*, "display_width")= int 0
##  $ 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 0

Notamos que la variable años no nos será de utilidad ya que no es una variable que vaya a servir como predictora, y tambien vemos que todas las variables son numéricas. Por otro lado, vale aclarar que la base de datos contiene datos agrupado por años que no es lo recomendable para realizar trabajos de regresión. Además la base de datos detroit esta compuesta por 13 variables, una es año y no la tendremos en cuenta y 13 observaciones. Lo que nos indica que tenemos una número excesivo de variables y muy pocas observaciones para construir un modelo.En general y como recomendación brindada en clase, se trabajan con bases que tengan al menos un mínimo de 10 observaciones por variable (excluyendo la dependiente) Sin embargo vamos a comenzar con analizar los datos de la base, para visualizar su correlación. Analizaremos los gráficos de dispersión y compararemos los r para medir el nivel de alineamiento de los valores con la recta. Para de esta manera, seleccionar los mejores para la regresión lineal múltiple.

El diagrama de dispersión nos muestra como es el tipo de asociacion si lineal, curvilinea, positiva o negativa e incluso si no hay asociación entre las variables.

El coeficiente de correlación de Pearson (r) nos indica que:

“si es -1, indica una correlación lineal perfectamente negativa entre dos variables”

“si es 0, indica que no hay correlación lineal entre dos variables”

“si es 1, indica una correlación lineal perfectamente positiva entre dos variables”

dispersion

Utilizaremos estos dos recursos para explicar el nivel de asociación lineal de las variables predictoras con la variable H.

plot(detroit$H,detroit$FTP,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Cantidad de policias full time")

correl1 <- cor(detroit$FTP,detroit$H)
correl1
## [1] 0.964061
regre1 <- lm(H~FTP,data=detroit)
summary(regre1)
## 
## Call:
## lm(formula = H ~ FTP, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.649 -2.244 -1.258  3.559  8.254 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -77.63027    8.63094  -8.994 2.11e-06 ***
## FTP           0.33745    0.02804  12.035 1.13e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.547 on 11 degrees of freedom
## Multiple R-squared:  0.9294, Adjusted R-squared:  0.923 
## F-statistic: 144.8 on 1 and 11 DF,  p-value: 1.129e-07

El gráfico de dispersión entre las variables (tasa de homocidio y cantidad de policias fulltime) nos permite visualizar una curva que luego mantiene su tendencia. El (r) nos indica que hay un nivel de asocicación lineal del 0.96; lo cual nos indica una correlación positiva fuerte. El p-value de 1.13e-07 del modelo lineal nos dice que es una relación estadisticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.93

plot(detroit$H,detroit$UNEMP,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Porcentaje de desempleados")

correl2 <- cor(detroit$UNEMP,detroit$H)
correl2
## [1] 0.2101417
regre2 <- lm(H~UNEMP,data=detroit)
summary(regre2)
## 
## Call:
## lm(formula = H ~ UNEMP, data = detroit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.128 -14.059  -1.297  10.354  26.462 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   16.673     12.735   1.309    0.217
## UNEMP          1.460      2.047   0.713    0.491
## 
## Residual standard error: 16.73 on 11 degrees of freedom
## Multiple R-squared:  0.04416,    Adjusted R-squared:  -0.04274 
## F-statistic: 0.5082 on 1 and 11 DF,  p-value: 0.4908

El gráfico de dispersión entre las variables (tasa de homocidio y porcentaje de desempleados) nos permite visualizar una curva que luego mantiene su tendencia. El r nos indica que hay un nivel de asocicación lineal del 0.21; lo cual nos indica una correlación positiva pero mucho mas débil que la variable anterior. El p-value de 0.491 del modelo lineal nos dice que no es una relación estadisticamente significativa y que no podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.04

plot(detroit$H,detroit$M,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Trabajadores en la Industria Manufacturera")

correl3 <- cor(detroit$M,detroit$H)
correl3
## [1] 0.5464227
regre3 <- lm(H~M,data=detroit)
summary(regre3)
## 
## Call:
## lm(formula = H ~ M, data = detroit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.689  -7.559  -3.890   9.953  22.507 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -74.87015   46.38238  -1.614   0.1348  
## M             0.17971    0.08305   2.164   0.0533 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.33 on 11 degrees of freedom
## Multiple R-squared:  0.2986, Adjusted R-squared:  0.2348 
## F-statistic: 4.682 on 1 and 11 DF,  p-value: 0.05334

El gráfico de dispersión entre las variables (tasa de homocidio y trabajadores en la industria manufacturera) nos permite visualizar una relación curvilinea. El r nos indica que hay un nivel de asocicación lineal del 0.54; lo cual nos indica una correlación positiva más débil que la primer variable. El p-value de 0.053 del modelo lineal nos dice que no es una relación estadisticamente significativa y que no podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.29

plot(detroit$H,detroit$LIC,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Cantidad de licencias de portación de armas")

correl4 <- cor(detroit$LIC,detroit$H)
correl4
## [1] 0.7262896
regre4 <- lm(H~LIC,data=detroit)
summary(regre4)
## 
## Call:
## lm(formula = H ~ LIC, data = detroit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.425  -4.930  -3.196   2.583  20.732 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.91138    6.62752   0.741  0.47418   
## LIC          0.03761    0.01073   3.504  0.00493 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.76 on 11 degrees of freedom
## Multiple R-squared:  0.5275, Adjusted R-squared:  0.4845 
## F-statistic: 12.28 on 1 and 11 DF,  p-value: 0.004933

El gráfico de dispersión entre las variables (tasa de homocidio y cantidad de licencias de portación de armas) nos permite visualizar una relación curvilinea. El r nos indica que hay un nivel de asocicación lineal del 0.72; lo cual nos indica una correlación positiva más débil que la primer variable. El p-value de 0.005 del modelo lineal nos dice que es una relación estadisticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.52

plot(detroit$H,detroit$GR,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Cantidad de armas de fuego registradas")

correl5 <- cor(detroit$GR,detroit$H)
correl5
## [1] 0.8162873
regre5 <- lm(H~GR,data=detroit)
summary(regre5)
## 
## Call:
## lm(formula = H ~ GR, data = detroit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.914  -2.901  -2.154   1.398  22.007 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.662063   5.708193   0.291 0.776337    
## GR          0.043003   0.009175   4.687 0.000664 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.886 on 11 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.636 
## F-statistic: 21.97 on 1 and 11 DF,  p-value: 0.0006642

El gráfico de dispersión entre las variables (tasa de homocidio y cantidad de armas de fuego registradas) nos permite visualizar una relación curvilinea. El r nos indica que hay un nivel de asocicación lineal del 0.81; lo cual nos indica una correlación positiva más débil que la primer variable. El p-value de 0.00067 del modelo lineal nos dice que es una relación estadisticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.66

plot(detroit$H,detroit$CLEAR,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Porcentaje de homicidios resueltos")

correl6 <- cor(detroit$CLEAR,detroit$H)
correl6
## [1] -0.9684603
regre6 <- lm(H~CLEAR,data=detroit)
summary(regre6)
## 
## Call:
## lm(formula = H ~ CLEAR, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.385 -1.636 -1.059  2.804  8.737 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 127.22163    8.00767   15.89 6.21e-09 ***
## CLEAR        -1.25352    0.09724  -12.89 5.55e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.264 on 11 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9323 
## F-statistic: 166.2 on 1 and 11 DF,  p-value: 5.553e-08

El gráfico de dispersión entre las variables (tasa de homocidio y porcentaje de homicidios resueltos) nos permite visualizar una dispersión negativa. El r nos indica que hay un nivel de asocicación lineal del -0.97; lo cual nos indica una correlación negativa fuerte. El p-value de 5.55e-08 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.94

plot(detroit$H,detroit$W,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Cantidad de hombres blancos")

correl7 <- cor(detroit$W,detroit$H)
correl7
## [1] -0.9469213
regre7 <- lm(H~W,data=detroit)
summary(regre7)
## 
## Call:
## lm(formula = H ~ W, data = detroit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7133 -4.7647 -0.5768  4.3803  9.1366 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.355e+02  1.140e+01   11.88 1.28e-07 ***
## W           -2.436e-04  2.493e-05   -9.77 9.33e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.502 on 11 degrees of freedom
## Multiple R-squared:  0.8967, Adjusted R-squared:  0.8873 
## F-statistic: 95.44 on 1 and 11 DF,  p-value: 9.328e-07

El gráfico de dispersión entre las variables (tasa de homocidio y cantidad de hombres blancos) nos permite visualizar una dispersión negativa. El r nos indica que hay un nivel de asocicación lineal del -0.95; lo cual nos indica una correlación negativa fuerte pero un poco más débil que la variable anterior. El p-value de 9.33e-07 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.90

plot(detroit$H,detroit$NMAN,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Trabajadores fuera de la industria manufacturera")

correl8 <- cor(detroit$NMAN,detroit$H)
correl8
## [1] 0.9559347
regre8 <- lm(H~NMAN,data=detroit)
summary(regre8)
## 
## Call:
## lm(formula = H ~ NMAN, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.223 -2.888 -1.341  3.425  7.684 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -86.2538    10.4073  -8.288 4.66e-06 ***
## NMAN          0.1653     0.0153  10.799 3.41e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.024 on 11 degrees of freedom
## Multiple R-squared:  0.9138, Adjusted R-squared:  0.906 
## F-statistic: 116.6 on 1 and 11 DF,  p-value: 3.411e-07

El gráfico de dispersión entre las variables (tasa de homocidio y Trabajadores fuera de la industria manufacturera) nos permite visualizar una dispersión positiva. El r nos indica que hay un nivel de asocicación lineal del 0.95; lo cual nos indica una correlación positiva fuerte pero un poco más débil que la primer variable. El p-value de 3.41e-07 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.91

plot(detroit$H,detroit$G,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Trabajadores del gobierno")

correl9 <- cor(detroit$G,detroit$H)
correl9
## [1] 0.9580545
regre9 <- lm(H~G,data=detroit)
summary(regre9)
## 
## Call:
## lm(formula = H ~ G, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.900 -3.857 -1.179  3.360  8.371 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -53.61314    7.23083  -7.415 1.34e-05 ***
## G             0.42386    0.03823  11.087 2.61e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.905 on 11 degrees of freedom
## Multiple R-squared:  0.9179, Adjusted R-squared:  0.9104 
## F-statistic: 122.9 on 1 and 11 DF,  p-value: 2.611e-07

El gráfico de dispersión entre las variables (tasa de homocidio y trabajadores del gobierno) nos permite visualizar una dispersión positiva. El r nos indica que hay un nivel de asocicación lineal del 0.96; lo cual nos indica una correlación positiva fuerte pero un poco más débil que la primer variable. El p-value de 2.61e-07 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.92

plot(detroit$H,detroit$HE,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Ingreso promedio por hora")

correl10 <- cor(detroit$HE,detroit$H)
correl10
## [1] 0.9134063
regre10 <- lm(H~HE,data=detroit)
summary(regre10)
## 
## Call:
## lm(formula = H ~ HE, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.672 -4.505 -1.459  1.682 18.971 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -36.001      8.438  -4.267  0.00133 ** 
## HE            15.484      2.081   7.442 1.29e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.966 on 11 degrees of freedom
## Multiple R-squared:  0.8343, Adjusted R-squared:  0.8192 
## F-statistic: 55.39 on 1 and 11 DF,  p-value: 1.289e-05

El gráfico de dispersión entre las variables (tasa de homocidio y ingreso promedio por hora) nos permite visualizar una dispersión positiva. Lo llamativo aquí es que hay un dato muy distinto a los demás que puede deberse a un error de carga de datos en la base, o a un caso atípico que puede estar perjudicando los calculos y el modelo que pueda construirse. El r nos indica que hay un nivel de asocicación lineal del 0.91; lo cual nos indica una correlación positiva fuerte pero un poco más débil que la primer variable. El p-value de 1.29e-05 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.83

plot(detroit$H,detroit$WE,main = "Gráfico de dispersión",
     pch = 16,
     xlab = "Tasa de Homicidio", 
     ylab = "Ingreso promedio semanal")

correl11 <- cor(detroit$WE,detroit$H)
correl11
## [1] 0.8881526
regre11 <- lm(H~WE,data=detroit)
summary(regre11)
## 
## Call:
## lm(formula = H ~ WE, data = detroit)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.709 -6.182 -1.877  3.361 15.987 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -33.05875    9.33585  -3.541  0.00462 ** 
## WE            0.34233    0.05341   6.410 5.01e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.865 on 11 degrees of freedom
## Multiple R-squared:  0.7888, Adjusted R-squared:  0.7696 
## F-statistic: 41.09 on 1 and 11 DF,  p-value: 5.011e-05

El gráfico de dispersión entre las variables (tasa de homocidio y ingreso promedio semanal) nos permite visualizar una dispersión positiva. Lo llamativo aquí como en el caso anterior, es que hay un dato muy distinto a los demás que puede deberse a un error de carga de datos en la base, o a un caso atípico que puede estar perjudicando los calculos y el modelo que pueda construirse. El r nos indica que hay un nivel de asocicación lineal del 0.88; lo cual nos indica una correlación positiva fuerte pero un poco más débil que la primer variable. El p-value de 5.01e-05 del modelo lineal nos dice que es una relación estadísticamente significativa y que podemos considerar esta variable como predictora. El coeficiente de determinación nos dice que la proporción de varianza de la Tasa de Homocidios que explica ésta relación es del 0.79

En relación a estos resultados, las variables que parecen tener mejor nivel de alineamiento con la recta y pueden ser utilizadas como variables predictoras son:

Cantidad de policias full time (FTP): r= 0.96 p-value=0,000000113 coeficiente de determinación= 0.93

Es llamativa esta correlación debido que indica que a mayor cantidad de policias, mayor cantidad de homicidios. Tambien es dudoso el resultado que arroja que se puede deber por errores de carga en los datos u otros problemas o no.

Porcentaje de homicidios resueltos (CLEAR) r= -0.97 p-value= 0,0000000555 coeficiente de determinación=0.94

Estos datos nos indican que a mayor casos resueltos menor tasa de homicidios.

Otro método para analizar la correlación de una manera más ágil sobre el conjunto de la base de datos en relación a la variable dependiente, es la siguiente: Los valores de Zero Order nos muestran el p-value de las correlaciones que vimos anteriormente.

model1 <- lm(H ~ FTP + LIC + GR + CLEAR + W + NMAN + G + HE + WE, data = detroit)

ols_correlations(model1)
##                Correlations                 
## -------------------------------------------
## Variable    Zero Order    Partial     Part  
## -------------------------------------------
## FTP              0.964      0.669     0.044 
## LIC              0.726      0.779     0.061 
## GR               0.816      0.494     0.028 
## CLEAR           -0.968     -0.615    -0.038 
## W               -0.947      0.431     0.023 
## NMAN             0.956      0.489     0.028 
## G                0.958      0.010     0.001 
## HE               0.913      0.286     0.015 
## WE               0.888      0.166     0.008 
## -------------------------------------------
  1. Para responder al punto 2 del ejercicio vamos a realizar una regresión múltiple con la función ols_regression que nos presenta los mismos resultados que lm pero con un formato de informe:
model_dos <- ols_regress(H ~ FTP + LIC + GR + CLEAR + W + NMAN + G + HE + WE, data = detroit)
model1
## 
## Call:
## lm(formula = H ~ FTP + LIC + GR + CLEAR + W + NMAN + G + HE + 
##     WE, data = detroit)
## 
## Coefficients:
## (Intercept)          FTP          LIC           GR        CLEAR            W  
##  -1.052e+02    8.401e-02    1.505e-02    4.475e-03   -3.810e-01    1.271e-04  
##        NMAN            G           HE           WE  
##   6.842e-02    3.731e-03    3.255e+00    4.647e-02

El método de selección automática de predictores Stepwise, es una combinación del método Forward y Backward. Revisa si agregando o sacando una variable puede afectar al modelo puede afectar el R2, debido a que el ingreso o exclusión de una variable puede afectar la importancia de las variables que ya estaban en el modelo. Para aplicar esta función necesitamos calcular la regresion con la función lm.(model1 <- lm(H ~ FTP + LIC + GR + CLEAR + W + NMAN + G + HE + WE, data = detroit) que usamos anteriormente.

ols_step_best_subset(model1)
##            Best Subsets Regression            
## ----------------------------------------------
## Model Index    Predictors
## ----------------------------------------------
##      1         CLEAR                           
##      2         LIC CLEAR                       
##      3         LIC CLEAR HE                    
##      4         FTP LIC CLEAR HE                
##      5         FTP LIC GR CLEAR HE             
##      6         FTP LIC GR CLEAR W HE           
##      7         FTP LIC GR CLEAR W NMAN HE      
##      8         FTP LIC GR CLEAR W NMAN HE WE   
##      9         FTP LIC GR CLEAR W NMAN G HE WE 
## ----------------------------------------------
## 
##                                                    Subsets Regression Summary                                                   
## --------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                          
## Model    R-Square    R-Square    R-Square     C(p)        AIC       SBIC        SBC        MSEP        FPE       HSP       APC  
## --------------------------------------------------------------------------------------------------------------------------------
##   1        0.9379      0.9323      0.9225    68.0254    78.4277    36.8285    80.1226    236.9508    20.9815    1.8184    0.0847 
##   2        0.9895      0.9874      0.9831     6.0272    57.3254    20.4205    59.5852     44.5278     4.1636    0.3759    0.0168 
##   3        0.9934      0.9912      0.9802     3.1649    53.2520    20.3956    56.0767     31.3968     3.0808    0.2945    0.0124 
##   4        0.9959      0.9938      0.9875     2.1331    49.2182    22.9540    52.6079     22.5583     2.3071    0.2380    0.0093 
##   5        0.9963      0.9936      0.9783     3.6147    49.8341    28.1432    53.7888     23.6600     2.5021    0.2853    0.0101 
##   6        0.9968      0.9936      0.9457     4.9996    49.9744    34.4587    54.4940     24.6076     2.6632    0.3462    0.0107 
##   7        0.9974      0.9937      0.9194     6.2670    49.3442    42.3678    54.4288     25.1253     2.7409    0.4242    0.0111 
##   8        0.9976      0.9927     -1.2856     8.0003    50.2373    51.1205    55.8868     30.7660     3.2964    0.6493    0.0133 
##   9        0.9976      0.9903      -1.374    10.0000    52.2359    59.7880    58.4504     46.1440     4.5944    1.2984    0.0185 
## --------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

El modelo que seleccionamos es el modelo 2 con dos predictores. Tambien respetando las recomendaciones de que con una base de pocas observaciones y muchas variables, es recomendable tener un modelo con pocas predictoras. Además se observa en los resultado que el incremento de R2 ajustado es mayor entre la primera y segunda predictora y disminuye en cuanto a la segunda y tercara predictora. Pareciera suficiente trabajar con dos predictoras sugeridas por R.

En consecuencia, vamos a analizar el modelo 2 si todos sus con sus variables predictoras.

model11 <- lm(H ~ CLEAR + LIC, data = detroit)
summary(model11)
## 
## 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-10

Los resultados nos indican que las variables tienen una alta significación y el p-value de 1.27e-10 no indica que es menor que 0,05 y por ende que es una relación estadísticamente significativa.

  1. en el modelo el Adjusted R-squared es igual a 0.9874. Este coeficiente de determinación ajustado es el que calcula la correlación del modelo pero penalizando por la cantidad de predictoras que integren el modelo. En este caso, la diferencia entre el R2 y el R2 ajustado (R-squared: 0.9895, Adjusted R-squared: 0.9874 ) es mínima.

Siendo el R2 la proporción que explicaría el modelo respecto de la variable dependiente.

5)La bondad de ajuste del modelo es Multiple R-squared: 0.9895

Indica la fuerza de la relación lineal entre las variables predictoras y la variable dependiente.Un resultado cercano a 1 indica una relación lineal perfecta.

4 y 6) Los supuestos que deben cumplirse para aplicar una regresión lineal son: que la relación sea lineal que los errores tengan una distribución normal que la varianza sea constante que la muestra sea independiente

Vamos a suponer que los datos de la muestra son independientes, y los demás supuestos los veremos con el gráfico de residuos R Student

ols_plot_resid_stud_fit(model11, print_plot = TRUE)

Si consideramos que en el gráfico se muestran como datos normales los que estan en la franja de 2 y -2 que representa un 95% de la distribución normal, es posible que el caso 2 que nos muestre outlier pueda serlo o entre dentro del 5% del longtail.

7)Si hay multicolinealidad en el modelo elegido veremos una foto distorisiona, es por ello que se realizan unos chequeos con:

*VIF (Variance inflation factor) y TOL (Tolerancia) que indica cuando un predictor tiene una fuerte relación lineal con el resto de los predictores. El valor de VIF comienza en 1 y no tiene límite superior. Un valor de 1 indica que no hay correlación entre una variable predictora dada y cualquier otra variable predictora en el modelo. Y Tolerancia si es menor que 0.10 indica multicolinealidad.

ols_vif_tol(model11)
##   Variables Tolerance     VIF
## 1     CLEAR 0.6921615 1.44475
## 2       LIC 0.6921615 1.44475
chequeo<-ols_vif_tol(model11)
View(chequeo)

Con estos resultados, entendemos que el modelo 2 seleccionado no tiene problemas de multicolinealidad.

  1. puntos atipicos
ols_plot_cooksd_bar(model11, print_plot = TRUE)

La distancia de Cook para analizar si hay valores atipicos o puntos de influencia nos vuelve a marcar el caso 2.

Consigna del Ejercicio 3

El bajo peso al nacer, definido como 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:

ID: Código de identificación LOW: Bajo peso al nacer. (0 = ≥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 uterine (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.

LOWBWT <- read_sav("LOWBWT.sav")
head(LOWBWT)
## # A tibble: 6 x 11
##      ID      LOW   AGE   LWT      RACE   SMOKE   PTL      HT      UI   FTV   BWT
##   <dbl> <dbl+lb> <dbl> <dbl> <dbl+lbl> <dbl+l> <dbl> <dbl+l> <dbl+l> <dbl> <dbl>
## 1     4   1 [Si]    28   120 3 [Otra]   1 [Si]     1  0 [No]  1 [Si]     0   709
## 2    10   1 [Si]    29   130 1 [Blanc~  0 [No]     0  0 [No]  1 [Si]     2  1021
## 3    11   1 [Si]    34   187 2 [Negra]  1 [Si]     0  1 [Si]  0 [No]     0  1135
## 4    13   1 [Si]    25   105 3 [Otra]   0 [No]     1  1 [Si]  0 [No]     0  1330
## 5    15   1 [Si]    25    85 3 [Otra]   0 [No]     0  0 [No]  1 [Si]     0  1474
## 6    16   1 [Si]    27   150 3 [Otra]   0 [No]     0  0 [No]  0 [No]     0  1588
  1. Calcular el riesgo relativo y los odds ratio de la variable dependiente con todas las variables dicotómicas. Analizar los resultados.

En nuestra base de datos contamos con 3 variable dicotómicas:

SMOKE

Calculamos el Riesgo relativo de la variable SMOKE, generando la tabla de contingencia:

tablesmoke <- table ("BAJO PESO"=LOWBWT$LOW,"FUMA"=LOWBWT$SMOKE)
marginsSMOKE <- addmargins (tablesmoke)
marginsSMOKE
##          FUMA
## BAJO PESO   0   1 Sum
##       0    86  44 130
##       1    29  30  59
##       Sum 115  74 189
#Riesgo o probabilidad de ocurrencia

((30/74)/(44/74))/((86/115)/(29/115))
## [1] 0.2299154

Es decir, que fumar aumenta el riesgo relativo de bajo peso al nacer en un 0.2299

y el ODD Ratio de la variable SMOKE, que nos compara la probabilidad de bajo peso al nacer habiendo fumado o no durante el embarazo:

fumar <- glm(LOW~SMOKE,data=LOWBWT,family="binomial")
summary(fumar)
## 
## 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: 4
or_glm(data = LOWBWT, model = fumar, incr = list(SMOKE=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1     SMOKE     2.022        1.082          3.801         1

En este caso, el odds ratio nos da 2.002. El p-value nos muestra que hay una significación de 0.0276 El oddsratio nos dice que fumar en el embarazo aumento la probabilidad de bajo peso en un 2.002. podemos decir que si hay un efecto de la variable predictora en relación con la variable dependiente.

HT Calculamos el Riesgo relativo de la variable HT

tableht <- table ("BAJO PESO"=LOWBWT$LOW,"HIPERTENSION"=LOWBWT$HT)
marginsht <- addmargins (tableht)
marginsht
##          HIPERTENSION
## BAJO PESO   0   1 Sum
##       0   125   5 130
##       1    52   7  59
##       Sum 177  12 189
#Riesgo o probabilidad de ocurrencia

(((7/12)/(5/12))/((52/177)/(125/177)))
## [1] 3.365385

y el odds ratio de hipertensión:

hipert <- glm(LOW~HT,data=LOWBWT,family="binomial")
summary(hipert)
## 
## 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: 4
or_glm(data = LOWBWT, model = hipert, incr = list(HT=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1        HT     3.365        1.028         11.829         1

En este caso, el riesgo relativo me da 3.36, al igual que el oddsratio. con un p-value que tiene significancia para la variable dependiente y el odds ratio nos indica que aumenta la probabilidad de bajo peso en un 3.36 si la persona es hipertensa

UI Calculamos el Riesgo relativo de la variable UI

tableiu <- table ("BAJO PESO"=LOWBWT$LOW,"IRRITABILIDAD UTERINA"=LOWBWT$UI)
marginsiu <- addmargins (tableiu)
marginsiu
##          IRRITABILIDAD UTERINA
## BAJO PESO   0   1 Sum
##       0   116  14 130
##       1    45  14  59
##       Sum 161  28 189
#Riesgo o probabilidad de ocurrencia

(((14/28)/(14/28))/((45/161)/(116/161)))
## [1] 2.577778

y el odds ratio de hipertensión:

irrit <- glm(LOW~UI,data=LOWBWT,family="binomial")
summary(irrit)
## 
## 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: 4
or_glm(data = LOWBWT, model = irrit, incr = list(UI=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1        UI     2.578        1.133          5.881         1

En este caso, el riesgo relativo me da 2.5, al igual que el oddsratio. con un p-value que tiene significancia para la variable dependiente y el odds ratio nos indica que aumenta la probabilidad de bajo peso en un 3.36 si la persona es hipertensa

  1. 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?

ODDS RATIO es el cociente de ocurrecia del evento en función de la aparición o no de la variable dependiente. Su interpretación nos permite medir el efecto de la variable predictora sobre la dependiente. Si es positivo nos habla de que la variable predictora es un factor de riesgo, es decir que su presencia aumenta la probabilidad de que suceda el evento. o negativa que indicaría que la variable predictora es un factor protector, es decir que su presencia disminuye la probabilidad de que suceda el evento. Se puede calcular con R a través de la función glm, el resultado significa:

*un resultado igual a 1, nos indica que no hay efecto

*un resultado mayor que 1, nos indica que aumenta en un x valor la presencia de la varibale predictora de que suceda el evento

*un resultado menor que 1, nos indica que disminuye el riesgo o oportunidad

  1. Calcule los odds ratio de cada una de las variables predictoras con la variable dependiente? Comentar.
LOWBWT = LOWBWT %>% mutate_each_(funs(scale(.) %>% as.vector),
                                 vars=c("AGE","LWT"))
## Warning: `mutate_each_()` was deprecated in dplyr 0.7.0.
## Please use `across()` instead.
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
head(LOWBWT)
## # A tibble: 6 x 11
##      ID     LOW   AGE      LWT     RACE   SMOKE   PTL      HT     UI   FTV   BWT
##   <dbl> <dbl+l> <dbl>    <dbl> <dbl+lb> <dbl+l> <dbl> <dbl+l> <dbl+> <dbl> <dbl>
## 1     4  1 [Si] 0.899 -0.321   3 [Otra]  1 [Si]     1  0 [No] 1 [Si]     0   709
## 2    10  1 [Si] 1.09   0.00606 1 [Blan~  0 [No]     0  0 [No] 1 [Si]     2  1021
## 3    11  1 [Si] 2.03   1.87    2 [Negr~  1 [Si]     0  1 [Si] 0 [No]     0  1135
## 4    13  1 [Si] 0.333 -0.811   3 [Otra]  0 [No]     1  1 [Si] 0 [No]     0  1330
## 5    15  1 [Si] 0.333 -1.47    3 [Otra]  0 [No]     0  0 [No] 1 [Si]     0  1474
## 6    16  1 [Si] 0.710  0.660   3 [Otra]  0 [No]     0  0 [No] 0 [No]     0  1588

Las variables dicotomicas ya fueron calculadas sus odds ratio, nos resta calcular por un lado la variable cualitativa de 3 modalidades, que se trabajará con variables dummy 1 y 2, ya que la combinación de los resultados de éstas es equivalente a los 3 modalidades de respuesta de la variable predictora original. Esto se realiza de la siguiente manera:

raza <- glm(LOW~factor(RACE),data=LOWBWT,family="binomial")
or_glm(data = LOWBWT, model = raza, incr = list(RACE=1))
##       predictor oddsratio ci_low (2.5) ci_high (97.5)          increment
## 1 factor(RACE)2     2.328        0.926          5.775 Indicator variable
## 2 factor(RACE)3     1.889        0.957          3.758 Indicator variable

Los dos odds ratio contiene en sus intervalos al 1, lo cual significa que no es una relación significativa y que podriamos descartar esta variable como predictora.

La variable AGE es una variable cuantitativa que el R traducirá a una variable dicotómica de la siguiente manera:

edad <-glm(LOW~AGE,data=LOWBWT,family="binomial")
or_glm(data = LOWBWT, model = edad, incr = list(AGE=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1       AGE     0.763        0.543          1.049         1

EL Oddsratio de EDAD es de 0.7 y el 1 está contenido en el intervalo, por lo que podriamos suponer que la variable edad no es una variable predictora del bajo peso.

peso <-glm(LOW~LWT,data=LOWBWT,family="binomial")
or_glm(data = LOWBWT, model = peso, incr = list(LWT=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1       LWT     0.651        0.438          0.922         1

El Oddsratio del peso de la madre es de 0.65, por lo que podriamos suponer que ésta variable puede aumentar la probabilidad de bajo peso.

premat <- glm(LOW~PTL,data=LOWBWT,family="binomial")
or_glm(data=LOWBWT, model=premat, incr=list(PTL=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1       PTL      2.23        1.218          4.284         1

Los antecedentes de partos prematuras parecen tener un efecto mayor en la probabilidad de bajo peso al nacer, de 2.23

consulta <- glm(LOW~FTV, data=LOWBWT,family="binomial")
or_glm(data=LOWBWT, model=consulta, incr=list(FTV=1))
##   predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1       FTV     0.874        0.632          1.174         1

En el caso de la cantidad de consultas realizadas durante el primer trimestre, el oddsratio es del 0.8 y el 1 está contenido en el intervalo de confianza por lo que pareciera no funcionar como una variable predictora.

  1. 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.

El método de selección automática se realiza de la misma manera que en regresión lineal, pero en logística se denomina STEPAIC y los criterios para seleccionar el mejor modelo estarán en la DEVIANCE, AIC y sus variantes. Para realizar la selección automática vamos a solicitarle a la regresión que aplique un modelo con todos las variables y luego incorporaremos las variables que nosotros consideramos como importantes para el modelo.

library(MASS)
logi_sin_pred = glm (LOW ~ 1, family = 'binomial', data =LOWBWT)

step_logi = stepAIC(logi_sin_pred, scope = list(upper = ~AGE+factor(PTL)+factor(FTV)+LWT+factor(RACE)+SMOKE+HT+UI), direction = c("both"), trace = T)
## Start:  AIC=236.67
## LOW ~ 1
## 
##                Df Deviance    AIC
## + factor(PTL)   3   218.80 226.80
## + LWT           1   228.69 232.69
## + UI            1   229.60 233.60
## + SMOKE         1   229.81 233.81
## + HT            1   230.65 234.65
## + factor(RACE)  2   229.66 235.66
## + AGE           1   231.91 235.91
## <none>              234.67 236.67
## + factor(FTV)   5   228.49 240.49
## 
## Step:  AIC=226.8
## LOW ~ factor(PTL)
## 
##                Df Deviance    AIC
## + LWT           1   213.53 223.53
## + AGE           1   213.77 223.77
## + HT            1   214.82 224.82
## + UI            1   214.89 224.89
## + SMOKE         1   215.89 225.89
## + factor(RACE)  2   214.61 226.61
## <none>              218.80 226.80
## + factor(FTV)   5   208.85 226.85
## - factor(PTL)   3   234.67 236.67
## 
## Step:  AIC=223.54
## LOW ~ factor(PTL) + LWT
## 
##                Df Deviance    AIC
## + HT            1   206.24 218.24
## + AGE           1   210.10 222.10
## + UI            1   210.69 222.69
## + SMOKE         1   210.93 222.93
## + factor(RACE)  2   208.97 222.97
## <none>              213.53 223.53
## + factor(FTV)   5   204.86 224.86
## - LWT           1   218.80 226.80
## - factor(PTL)   3   228.69 232.69
## 
## Step:  AIC=218.24
## LOW ~ factor(PTL) + LWT + HT
## 
##                Df Deviance    AIC
## + UI            1   202.59 216.59
## + AGE           1   203.34 217.34
## + SMOKE         1   203.80 217.80
## + factor(RACE)  2   202.21 218.21
## <none>              206.24 218.24
## + factor(FTV)   5   198.65 220.65
## - HT            1   213.53 223.53
## - LWT           1   214.82 224.82
## - factor(PTL)   3   221.14 227.14
## 
## Step:  AIC=216.6
## LOW ~ factor(PTL) + LWT + HT + UI
## 
##                Df Deviance    AIC
## + AGE           1   200.10 216.10
## + SMOKE         1   200.27 216.27
## + factor(RACE)  2   198.40 216.40
## <none>              202.59 216.59
## - UI            1   206.24 218.24
## + factor(FTV)   5   195.81 219.81
## - LWT           1   209.90 221.90
## - HT            1   210.69 222.69
## - factor(PTL)   3   216.61 224.61
## 
## Step:  AIC=216.1
## LOW ~ factor(PTL) + LWT + HT + UI + AGE
## 
##                Df Deviance    AIC
## + SMOKE         1   197.91 215.91
## <none>              200.10 216.10
## - AGE           1   202.59 216.59
## + factor(RACE)  2   196.91 216.91
## - UI            1   203.34 217.34
## - LWT           1   205.78 219.78
## + factor(FTV)   5   194.22 220.22
## - HT            1   207.62 221.62
## - factor(PTL)   3   215.48 225.48
## 
## Step:  AIC=215.91
## LOW ~ factor(PTL) + LWT + HT + UI + AGE + SMOKE
## 
##                Df Deviance    AIC
## + factor(RACE)  2   192.45 214.45
## <none>              197.91 215.91
## - SMOKE         1   200.10 216.10
## - AGE           1   200.27 216.27
## - UI            1   201.04 217.04
## - LWT           1   203.37 219.37
## + factor(FTV)   5   193.23 221.23
## - HT            1   205.32 221.32
## - factor(PTL)   3   211.78 223.78
## 
## Step:  AIC=214.45
## LOW ~ factor(PTL) + LWT + HT + UI + AGE + SMOKE + factor(RACE)
## 
##                Df Deviance    AIC
## - AGE           1   193.59 213.59
## <none>              192.45 214.45
## - UI            1   195.67 215.67
## - factor(RACE)  2   197.91 215.91
## - SMOKE         1   196.91 216.91
## - LWT           1   198.05 218.05
## - HT            1   199.64 219.64
## - factor(PTL)   3   203.95 219.95
## + factor(FTV)   5   188.60 220.60
## 
## Step:  AIC=213.59
## LOW ~ factor(PTL) + LWT + HT + UI + SMOKE + factor(RACE)
## 
##                Df Deviance    AIC
## <none>              193.59 213.59
## + AGE           1   192.45 214.45
## - UI            1   197.17 215.17
## - factor(RACE)  2   200.27 216.27
## - SMOKE         1   198.40 216.40
## - factor(PTL)   3   204.22 218.22
## - LWT           1   200.29 218.29
## - HT            1   200.94 218.94
## + factor(FTV)   5   189.58 219.58
predictores <- glm(LOW ~ 1,family = 'binomial', data = LOWBWT)

birthwt.step <- stepAIC(predictores,scope = list(upper = ~AGE+LWT+factor(RACE)+SMOKE+PTL+HT+UI+FTV, lower = ~1),
                        direction = c("both"),trace = T)
## Start:  AIC=236.67
## LOW ~ 1
## 
##                Df Deviance    AIC
## + PTL           1   227.89 231.89
## + LWT           1   228.69 232.69
## + UI            1   229.60 233.60
## + SMOKE         1   229.81 233.81
## + HT            1   230.65 234.65
## + factor(RACE)  2   229.66 235.66
## + AGE           1   231.91 235.91
## <none>              234.67 236.67
## + FTV           1   233.90 237.90
## 
## Step:  AIC=231.89
## LOW ~ PTL
## 
##                Df Deviance    AIC
## + LWT           1   223.41 229.41
## + HT            1   223.58 229.58
## + AGE           1   224.27 230.27
## + factor(RACE)  2   222.53 230.53
## + SMOKE         1   224.78 230.78
## + UI            1   224.89 230.89
## <none>              227.89 231.89
## + FTV           1   227.30 233.30
## - PTL           1   234.67 236.67
## 
## Step:  AIC=229.41
## LOW ~ PTL + LWT
## 
##                Df Deviance    AIC
## + HT            1   215.96 223.96
## + factor(RACE)  2   217.68 227.68
## + SMOKE         1   220.54 228.54
## + AGE           1   221.05 229.05
## + UI            1   221.23 229.23
## <none>              223.41 229.41
## + FTV           1   223.12 231.12
## - LWT           1   227.89 231.89
## - PTL           1   228.69 232.69
## 
## Step:  AIC=223.96
## LOW ~ PTL + LWT + HT
## 
##                Df Deviance    AIC
## + factor(RACE)  2   210.85 222.85
## + UI            1   213.01 223.01
## + SMOKE         1   213.15 223.15
## <none>              215.96 223.96
## + AGE           1   214.01 224.01
## + FTV           1   215.84 225.84
## - PTL           1   221.14 227.14
## - HT            1   223.41 229.41
## - LWT           1   223.58 229.58
## 
## Step:  AIC=222.85
## LOW ~ PTL + LWT + HT + factor(RACE)
## 
##                Df Deviance    AIC
## + SMOKE         1   204.90 218.90
## + UI            1   207.73 221.73
## <none>              210.85 222.85
## + AGE           1   209.81 223.81
## - factor(RACE)  2   215.96 223.96
## + FTV           1   210.82 224.82
## - PTL           1   216.29 226.29
## - HT            1   217.68 227.68
## - LWT           1   218.69 228.69
## 
## Step:  AIC=218.9
## LOW ~ PTL + LWT + HT + factor(RACE) + SMOKE
## 
##                Df Deviance    AIC
## + UI            1   201.99 217.99
## <none>              204.90 218.90
## + AGE           1   204.11 220.11
## - PTL           1   208.25 220.25
## + FTV           1   204.88 220.88
## - SMOKE         1   210.85 222.85
## - factor(RACE)  2   213.15 223.15
## - HT            1   211.55 223.55
## - LWT           1   211.62 223.62
## 
## Step:  AIC=217.99
## LOW ~ PTL + LWT + HT + factor(RACE) + SMOKE + UI
## 
##                Df Deviance    AIC
## <none>              201.99 217.99
## - PTL           1   204.22 218.22
## - UI            1   204.90 218.90
## + AGE           1   201.43 219.43
## + FTV           1   201.93 219.93
## - SMOKE         1   207.73 221.73
## - LWT           1   208.11 222.11
## - factor(RACE)  2   210.31 222.31
## - HT            1   209.46 223.46
View(birthwt.step$anova)

Los resultados de este procedimiento nos indica que el mejor modelo es el último con AIC=217.9 que incluye: LOW ~ PTL + LWT + HT + factor(RACE) + SMOKE + UI

  1. 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?

LOW ~ PTL + LWT + HT + factor(RACE) + SMOKE + UI

modelo_logistico  = glm(LOW ~ factor(PTL) + LWT + HT + factor(RACE)+ SMOKE + UI,family = "binomial", data = LOWBWT)
summary(modelo_logistico)
## 
## Call:
## glm(formula = LOW ~ factor(PTL) + LWT + HT + factor(RACE) + SMOKE + 
##     UI, 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)    -2.1990     0.4084  -5.385 7.26e-08 ***
## factor(PTL)1    1.4579     0.5074   2.873  0.00406 ** 
## factor(PTL)2    0.2738     0.9808   0.279  0.78007    
## factor(PTL)3  -14.7446   882.7435  -0.017  0.98667    
## LWT            -0.5251     0.2177  -2.412  0.01588 *  
## HT              1.8982     0.7175   2.645  0.00816 ** 
## factor(RACE)2   1.2489     0.5352   2.333  0.01962 *  
## factor(RACE)3   0.7967     0.4474   1.781  0.07493 .  
## SMOKE           0.8854     0.4094   2.163  0.03057 *  
## UI              0.8942     0.4696   1.904  0.05691 .  
## ---
## 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: 13

Según nuestro modelo las variables predictoras más importantes son PTL (Antecendentes de embarazos prematuros), HT(Hipertensión arterial), SMOKE (Fumo durante el embarazo)

  1. Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámetros?

Los supuestos de la estimación de los parametros son los mismos que en la regresión lineal pero en regresión logistica dada la naturaleza de las variables que analizan, el calculo de los parametros no utiliza el metodo de los mínimos cuadrados (minima diferencia con la recta) sino el método de maxima verosimilitud. El criterio del método es maximizar la probabilidad conjunta de la muestra obtenida, bajo el supuesto de que la muestra que tenemos es la que tuvo mayor probabilidad de ser seleccionada. La función de verosimilitud permite además por sus propiedades calcular con el método Deviance, que es el logaritmo de la función de verosimilitud multiplicado por dos. La deviance se utiliza en el test de cociente de máxima verosimilitud o test de bondad de ajuste, que permite comparar y evaluar distintos modelos predictivos para saber cual ajusta mejor.

  1. Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.

El metodo que nos indica la bondad de ajuste es AIC mencionado anteriormente.

  1. Indicar porcentaje de casos bien predichos por el modelo.

Utilizamos la curva ROC para evaluar el modelo segun el test de sensibilidad y de especificidad

Calculo de probabilidades se realiza a través de la siguiente funcion.

predic_modelo <- plogis(predict(modelo_logistico, LOWBWT))
head(predic_modelo)
##         1         2         3         4         5         6 
## 0.8811894 0.2128321 0.7008686 0.9152869 0.5650202 0.1481861

Para medir los errores de calculo, calculamos primero el punto optimo de corte

optCutOff_modelo <- optimalCutoff(LOWBWT$LOW, predic_modelo) 
View(optCutOff_modelo)

El punto de corte óptimo es 0.43. Ahora vamos a medir los errores de clasificacion del modelo a través de:

misClassError(LOWBWT$LOW, predic_modelo, threshold = optCutOff_modelo)
## [1] 0.2328

El modelo tiene un error de clasificación de 0.23 Ahora lo visualimos con la curva de ROC, que grafica la sensibilidad y especificidad del modelo

plotROC(LOWBWT$LOW, predic_modelo)

En término de clasficación de 0.5 test muy malo 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 perfecto

El resultado de la curva de ROC nos indica que nuestro modelo es 0.77, con lo cual podemos considerar que es un modelo bueno estadisticamente hablando.