Entrega Final
Jonathan Rougier
07/10/2021
INTRODUCCION
A continuación se comparten los ejercicios de evaluación para la materia de Estadistica 2. El exámen final consiste en tres ejercicios que abarcan los temas abordados durante la cursada. Los temas y ejercicios son los siguientes:
ANOVA - Ejercicio pasteles.
Regresion lineal multiple - Ejercicio Detroit.
Regresion logistica- Ejercicio Bajo peso al nacer.
Para llevar a cabo el examen se utilizarán las siguientes librerías.
library(olsrr) #Para modelos de regresion
library(haven) #Para cargar archivos con formatos foraneos
library(dplyr) #Para manipular datos
library(plyr) #Para manipular datos
library(ggpubr) #Para graficos
library(rstatix) #Para herramientas estadisticas
library(tidyverse) #Para organizar datos
library(gridExtra) #Para crear grids y tablas
library(ggplot2) #Para graficos
library(goftest) #Para test estadisticos como V de Cramer
library(tibble) #Para dataframes
library(nortest)
library(oddsratio) #Para calculo de oddsratio
library(survey) #Para trabajar con encuestas
library(epiR) #Para correr oddsratios
library(MASS) #Para herramientas estadisticas
library(InformationValue) #Para diagnosticos de probabilitdad
library(reshape2) #Para manipular los datos
library(lmtest) #Para testeo de supuestos de modelos lineales
library(DescTools) #Para testeo de Hosmer Lemeshow
library(corrplot) #Para visualizar la matriz de correlacionEJERCICIO 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
Los resultados obtenidos son los que se muestran en el archivo terapia.sav
1.
Existen diferencias estadísticamente significativas debidas al tipo de harina utilizada?
#Creamos los vectores correspondientes a los datos obtenidos en la consigna.
#Vecttor harina para uso general
harinagen <- c(0.90, 0.87, 0.90, 0.86,0.89,0.91,0.93,0.88,0.87,0.79,0.82,0.80)
#Vector harina torta
harinator <- c(0.91,0.90,0.80,0.88,0.82,0.83,0.86,0.85,0.80,0.86,0.85,0.85)
#Vector concentracion de azucar
concazu<- c(rep(0,3), rep(50,3), rep(75,3), rep(100,3))
#Data frame con las variables
tortaframe <- data.frame(concazu,harinagen,harinator)
View(tortaframe)#Creamos un data frame con los dos vectores que utilizaremos pare el primer ejercicio
#Vector densidad de cada harina
densidad <- c(harinagen,harinator)
#Vector tipo de harina
tipoharina <- c(rep("general",12), rep( "tortas",12))
#Data frame
tortaframe1<- data.frame(tipoharina,densidad)
View(tortaframe1)Para poder efectuar el análisis de varianza debemos comprobar primero que nuestras variables cumplan los supuestos necesarios para realizar un test de ANOVA.
El test de ANOVA puede realizarse bajo tres supuestos:
Las respuestas para cada factor deben tener una distribucion normal.
Las distribuciones deben tener la misma variancia.
La muestras deben ser independientes.
Siguiendo los supuestos necesarios, debemos comprobar que la distribucion de nuestras observaciones sea normal y que tengan la misma variancia.
Para corroborar el primer supuesto podemos aplicar el test de Shapiro-Wilk, la cual plantea como hipotesis nula que una muestra proviene de una poblacion normalmente distribuida. Rechazaremos la hipotesis nula si el valor de P<0.05.
Debajo observamos como a partir del test de Shapiro aceptamos la hipótesis nula con un P valor mayor a 0.05, esto quiere decir que la muestras se distribuyen normalmente. Nuestra muestra cumple el supuesto de normalidad.
#Test de shapiro
tortaframe1 %>% group_by(factor(tipoharina)) %>%
shapiro_test(densidad)## # A tibble: 2 x 4
## `factor(tipoharina)` variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 general densidad 0.924 0.319
## 2 tortas densidad 0.950 0.643El segundo supuesto a corroborar es el de igualdad de las variancias o homoecedasticidas. Es entonces que debemos testear que la varianza de los errores se mantenga igual a lo largo de las muestras. Para ello usamos el test de Levene cuya hipotesis nula nos dice que las varianzas de los grupos son iguales entre si.
En nuestro caso observamos que se cumple la condicón de homoecedasticidad. Se acepta la hipótesis nula que nos dice que las varianzas son iguales debido a que P>0.05.
#Test de levene para homoecedasticidad de la densidad segun tipo de harina.
tortaframe1 %>% levene_test(densidad ~ as.factor(tipoharina))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 22 0.528 0.475#Forma grafica de corroborar la normalidad
ggqqplot(tortaframe1, "densidad", facet.by="tipoharina")
Una vez definida la normalidad de la distribución y la igualdad de las varianzas, realizamos el test de ANOVA.
El test de ANOVA o analisis de la varianza unidireccional nos permite comparar dos medias de grupos independientes utilizando la distribucion F de Snedecor. El principio de este metodo utiliza una variable dependiente y una independiente, en donde la hipotesis nula nos dice que las dos medias son iguales.
#Analisis de varianza de tipo de harina sobre densidad
anova1<- aov(densidad ~ as.factor(tipoharina), tortaframe1)
summary(anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tipoharina) 1 0.00184 0.001837 1.16 0.293
## Residuals 22 0.03486 0.001584Con una Pr(>F) mayor a 0.05 aceptamos la hipótesis nula. Esto significa que debido a que F= 0.293, el tipo de harina utilizada no explicaría la diferencías en las densidades de las tortas.
2.
Existen diferencias estadísticamente significativas debidas al nivel de endulzamiento utilizado?
En el caso que querramos determinar si la concentración de azucar afecta a la densidad debemos repetir el proceso. Primero modificaremos la matriz para luego comprobar la normalidad y homoecedasticidad de la muestras. Segundo procederemos a realizar el análisis de varianza.
#Usamos la funcion melt para cambiar el formato de la tabla y poder analizar los resultados
tablaendu<- melt(tortaframe, id.vars = "concazu" ) %>% as.data.frame()
colnames(tablaendu)[3] <- "densidad"
colnames(tablaendu)[2] <- "tipoharina"
tablaendu## concazu tipoharina densidad
## 1 0 harinagen 0.90
## 2 0 harinagen 0.87
## 3 0 harinagen 0.90
## 4 50 harinagen 0.86
## 5 50 harinagen 0.89
## 6 50 harinagen 0.91
## 7 75 harinagen 0.93
## 8 75 harinagen 0.88
## 9 75 harinagen 0.87
## 10 100 harinagen 0.79
## 11 100 harinagen 0.82
## 12 100 harinagen 0.80
## 13 0 harinator 0.91
## 14 0 harinator 0.90
## 15 0 harinator 0.80
## 16 50 harinator 0.88
## 17 50 harinator 0.82
## 18 50 harinator 0.83
## 19 75 harinator 0.86
## 20 75 harinator 0.85
## 21 75 harinator 0.80
## 22 100 harinator 0.86
## 23 100 harinator 0.85
## 24 100 harinator 0.85#Test de Shapiro para normalidad de las variables endulzantes y densidad
tablaendu %>% group_by(as.factor(concazu)) %>% shapiro_test(densidad)## # A tibble: 4 x 4
## `as.factor(concazu)` 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#Visualizacion de normalidad
ggqqplot(tablaendu, "densidad", facet.by = "concazu" )
Observamos que se cumple la condicion de normalidad para casi todos los factores P>0.05. Sin embargo, el factor de 0 endulzante no cumple la condicion de normalidad (P=0.0122) . Es por ello que deberemos ser cuidadosos al momento de analizar los datos ya que no deberiamos aplicar pruebas parametricas en condiciones de no normalidad.
#Test de levene para homoecedasticidad de densidad segun endulzante.
tablaendu %>% levene_test(densidad ~ as.factor(concazu))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 20 0.0530 0.983A partir del Test de Levene observamos que nuestra P es igual a 0.983. Se acepta la hipótesis nula de normalidad de las varianzas.
Corroborados los supuestos, podemos correr nuestro analisis de ANOVA.
#Analisis de varianza
anovaconcazu<- aov(densidad ~ as.factor(concazu), tablaendu)
summary(anovaconcazu)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(concazu) 3 0.008713 0.002904 2.076 0.136
## Residuals 20 0.027983 0.001399Al igual que con el tipo de harina, se acepta la hipotesis nula. No hay diferencias significativas en la densidad según concentración de azucar (F=0.136 > 0.05).
3.
Existe interacción estadísticamente significativa entre los factores considerados?
Si quisieramos ver de qué manera interactúan la concentración de azucar y el tipo de harina deberíamos realizar un gráfico de interacción. En estos casos que introducimos dos variables independientes interactuando entre si estariamos hablando de un ANOVA de dos direcciones o bidireccional.
Si bien el principio y supuestos del ANOVA bidireccional es el mismo que el del ANOVA unidireccional, nos encontramos con dos variables de tipo factor o categoricas y una variable dependiente de tipo cuantitativa.
Es necesario entender el modo en que las variables independientes interactuan entre si. Cuando las variables interactuan nuestro ANOVA se calculara de manera distinta a cuando no interactuan entre ellas.
Mediante el grafico de interaccion podemos identificar el modo en que las variables interactuan. A continuacion se muestra el grafico de interaccion.
#Grafico de interaccion de medias
with(tablaendu, interaction.plot(concazu,tipoharina, densidad, fun = mean , main="Grafico de interaccion de medias", ylab = "Media de densidad", xlab="Concentracion de azucar"))
A partir del mismo vemos que existe una interacción entre el tipo de harina y la concentración de azucar. Lo cual nos da el indicio de como debera ser ejecutado nuestro analisis de ANOVA.
###ANOVA de dos factores con interaccion. Observando el grafico, estos son los factores que debemos mirar
anovados2 <- aov(densidad ~ as.factor(tipoharina)+as.factor(concazu), tablaendu)
summary(anovados2)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tipoharina) 1 0.001838 0.001837 1.335 0.262
## as.factor(concazu) 3 0.008713 0.002904 2.110 0.133
## Residuals 19 0.026146 0.001376Los coeficientes de F de Snedecor nos arrojan resultados por encima de 0.05. Entonces podemos decir que la relacion no es significativa. Aceptamos la hipotesis nula.
### ANOVA de dos factores sin interaccion
anovados<- aov(densidad ~ as.factor(tipoharina)+as.factor(concazu), tablaendu)
summary(anovados)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tipoharina) 1 0.001838 0.001837 1.335 0.262
## as.factor(concazu) 3 0.008713 0.002904 2.110 0.133
## Residuals 19 0.026146 0.001376En caso de que no hubiera interaccion el ANOVA de dos factores sin interaccion da un resultado similar, aunque no seria el test correcto para nuestro problema.
4.
Realizar los gráficos de medias e interpretar
A continuacion se grafican los procedimientos expuestos en los puntos anteriores. Si bien las visualizaciones ayudan a entender tendencias de manera clara, a veces pueden llevarnos a conclusiones equivocadas. Por ello, se interpretara los graficos pero se concluira a partir de los Test de ANOVA hechos en los ejercicios anteriores.
#Graficos de media para densidad segun tipo de harina
ggplot(tortaframe1, aes(x=densidad, y=tipoharina, fill=densidad)) + geom_boxplot() +
ggtitle("Densidad segun tipo de harina")+
stat_summary(fun = "mean", geom = "point", shape=8, size= 2, color="red") +
theme_cleveland()
Se observa que la media de densidad para la harina general y la harina para torta es muy similar, lo cual confirma el resultado del test que dice que el tipo de harina no explicaria una diferencia en la densidad del producto.
#Grafcos de media para densidad segun nivel de endulzamiento
ggplot(tablaendu, aes(x=densidad, y=as.factor(concazu), fill=densidad)) + geom_boxplot() +
stat_summary(fun = "mean", geom = "point", shape=8, size= 2, color="red") +
theme_cleveland()
En el caso del impacto que tiene el nivel de endulzamiento sobre la densidad de las tortas, se observan leves tendencias hacia una media cercana a una densidad del 0.90 en aquellas endulzadas con 0, 50 y 75. Uno podria deducir que hay cierta influencia del endulzante sobre la densidad, pero sin embargo el test de ANOVA no mostraba diferencias significativas.
#Graficos de media para densidad segun nivel de endulzamiento y tipo de harina
ggplot(tablaendu, aes(x=densidad, y= interaction(as.factor(concazu), tipoharina), fill=densidad)) + geom_boxplot() +
ggtitle("Densidad segun nivel de endulzamiento y tipo de harina")+
stat_summary(fun = "mean", geom = "point", shape=8, size= 2, color="red") +
theme_cleveland()
Cuando observamos las medias de densidad segun el tipo de harina y endulzante, la combinacion de cualquier tipo de endulzante con la harina de tortas tiende a arrojar una densidad con menos dispersion. En cambio, la harina general es suceptible al endulzante produciendo asi densidades con mayor dispersion. Retomando el test de ANOVA y el grafico de interaccion, vemos que el efecto de la azucar 100 con la harina general arroja densidades mas bajas.
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 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 (archivo p301.sav):
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.
#Cargamos la base
detroit<- read_sav("p301.sav")1.
Cuál es el nivel de asociación lineal de las variables predictoras con la variable H? Comentar.
Para entender cual es el nivel de asociacion lineal de las variables predictoras con la variable explicada debemos observar el coeficiente de correlación que índica la fuerza con que dos variables se asocian.
Este coeficiente se expresa entre 1 y -1, siendo que 1 es la máxima asociación y que -1 es la máxima asociación inversa. El valor 0 o cercanos a 0 expresa la ausencia de asociación.
Para dar cuenta de la asociación correremos una correlacion de cada una de las variables independientes con la variable dependiente.
#Creamos el modelo de regresion lineal para usarlos en la correlacion
modelodetroit<- lm( H ~ FTP + UNEMP + M + LIC + GR + CLEAR + W + NMAN + G + HE + W, data = detroit )
#Corremos la correlacion de las variables independientes con la dependiente
ols_correlations(modelodetroit)## Correlations
## -------------------------------------------
## Variable Zero Order Partial Part
## -------------------------------------------
## FTP 0.964 0.392 0.014
## UNEMP 0.210 0.608 0.025
## M 0.546 0.205 0.007
## LIC 0.726 0.898 0.067
## GR 0.816 0.178 0.006
## CLEAR -0.968 -0.803 -0.044
## W -0.947 -0.486 -0.018
## NMAN 0.956 -0.154 -0.005
## G 0.958 -0.664 -0.029
## HE 0.913 0.725 0.034
## -------------------------------------------A partir del resultado de las correlacion de las variables independientes con la dependiente podemos decir que:
Hay una fuerte correlacion positiva en las variables FTP (Cantidad de policias full time cada 100.00 habitantes) con un R 0.964 , NMAN( Cantidad de trabajadores fuera de la industria manufacturera) con un R de 0.956, G (Cantidad de trabajadores del gobierno) con un R de 0.958 y HE (Ingreso promedio por hora) con un R de 0.913.
Hay una fuerte correlacion negativas en las variables CLEAR (Porcentaje de homicidios resueltos con arresto del responsable) con un R de -0.968 y W (Cantidad de hombres blancos en la poblacion) con un R de -0.947.
Las variables UNEMP (Porcentaje de desempleados) con un R de 0.210, M (Cantidad de trabajadores en la industria manufacturera en miles) con un R de 0.546, LIC (Cantidad de licencias de portacion de armas cada 100.000 habitantes) con un R de 0.726 y GR (Cantidad de armas de fuego registradas cada 100.000 habitantes) con un R de 0.816, presentan una baja o incluso casi nula correlacion.
#Visualizacion matriz de correlacion
matrizcorreladetroit<- cor(detroit) #correlacion de todas las variables
corrplot(matrizcorreladetroit, method = "color" ) #grafico
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.)
El mayor inconveniente de trabajar con grandes set de variables que podrian actuar como independientes es definir cuales son las que mejor se ajustan para nuestro modelo. Los metodos de seleccion automatica tienen como fin reducir el subset de predictores a solo aquellos que mejor predicen la variable dependiente. Los metodos que se pueden utilizar son:
Forward
Backward
Stepwise
Para seleccionar los mejores predictores entre las variables independientes de nuestro set se utilizara el metodo de stepwise. Este metodo secuencial de seleccion es una combinacion entre los metodos de seleccion forward y backwards.
El metodo forward incorpora variables independientes a medida que estas ingresan en la secuencia. El primer paso del modelo consiste en seleccionar la variable independiente que presente el valor mas alto en la suma de cuadrados de la regresion. Una vez incorporada esta primera variable, el metodo continua incorporando las variables que mas aportan a explicar la dependiente.
El metodo backward procede con una logica similar pero eliminando las variables que arrojen el valor mas pequeno de SSR en el modelo. Este metodo elimina variables del modelo hasta dejar solo las variables significativas.
Finalmente, el metodo stepwise combina los elementos anteriores pero incorporando a la vez criterios tales como el R2 ajustado, el Cp de Mallows, los indicadores de Akaike y criterios Bayesianos. Lo cual no solo automatiza el proceso sino que nos da herramientas para tomar la decision de cual sera el mejor de los modelos segun nuestros criterio.
Un bueno modelo sera aquel que:
Su R2 ajustado explique lo mas posible la variacion de la variable/s indepente/s sobre la dependiente.
El Cp de Mallows sea lo menor o lo mas cercano a la cantidad de parametros que tiene el modelo.
El indicador de Akaike sea el mas bajo.
El indicador de BIC sea el mas bajo de los posibles modelo (al igual que Akaike).
Teniendo en cuenta los metodos y criterios de seleccion podemos correr el analisis de stepwise para nuestro set de variables y asi obtener los modelos que mejor se ajusten.
#Los mejores modelos posibles
bestdetroit<- ols_step_best_subset(modelodetroit, pent=0.10, prem=0.15 )#Observamos el resultado del stepwise
bestdetroit## Best Subsets Regression
## ---------------------------------------------------
## Model Index Predictors
## ---------------------------------------------------
## 1 CLEAR
## 2 LIC CLEAR
## 3 LIC CLEAR HE
## 4 FTP LIC CLEAR HE
## 5 FTP UNEMP LIC CLEAR HE
## 6 UNEMP LIC CLEAR W G HE
## 7 FTP UNEMP LIC CLEAR W G HE
## 8 FTP UNEMP LIC GR CLEAR W G HE
## 9 FTP UNEMP M LIC GR CLEAR W G HE
## 10 FTP UNEMP M LIC GR CLEAR W NMAN G HE
## ---------------------------------------------------
##
## 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 107.7086 78.4277 36.4016 80.1226 236.9508 20.9815 1.8184 0.0847
## 2 0.9895 0.9874 0.9831 12.7387 57.3254 18.1515 59.5852 44.5278 4.1636 0.3759 0.0168
## 3 0.9934 0.9912 0.9802 7.3715 53.2520 16.7608 56.0767 31.3968 3.0808 0.2945 0.0124
## 4 0.9959 0.9938 0.9875 4.7777 49.2182 18.1386 52.6079 22.5583 2.3071 0.2380 0.0093
## 5 0.9975 0.9956 0.9813 3.7934 44.9261 22.7160 48.8808 16.2200 1.7153 0.1956 0.0069
## 6 0.9985 0.9971 0.9854 3.7495 39.7003 27.2038 44.2199 11.1645 1.2083 0.1571 0.0049
## 7 0.9988 0.9972 0.9863 5.1685 38.6146 31.7429 43.6991 11.0069 1.2007 0.1858 0.0048
## 8 0.9989 0.9967 0.9745 7.0883 40.1245 42.6806 45.7740 14.1329 1.5142 0.2983 0.0061
## 9 0.9989 0.9956 -3.8595 9.0486 41.8747 54.7430 48.0892 20.7958 2.0706 0.5852 0.0084
## 10 0.9989 0.9936 -4.2198 11.0000 43.5628 67.1704 50.3422 40.6055 3.1641 1.7139 0.0128
## ---------------------------------------------------------------------------------------------------------------------------------
## 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 CriteriaA partir del metodo de stepwise obtenemos 10 modelos que incluyen las variables independientes con las que estamos trabajando. Estos modelos contemplan una probabilidad de entrada de 0.10 y una probabilidad de salida del modelo de 0.15. Es asi como observamos que:
El modelo 1 con un predictor (CLEAR) explica el 94% de la variacion de la variable independiente sobre la dependiente. A pesar de su alto valor explicativo, la posibilidad de sesgos es bastante alta y se observa en los altos valores del indicador de AIC.
El modelo 2 con dos predictores (CLEAR y LIC) explica casi el 99% de la variacion de las variables independientes sobre la dependiente. Este modelo presenta mejores indicadores en AIC y Cp de Mallows.
El resto de los modelos con 3 variables o mas explican el 99% o mas de la variacion de las variables independientes sobre la dependiente. A pesar de su alto valor explicativo y sus mejores indicadores de Cp de Mallows y AIC, este modelo incorpora varios parametros pero sin ser significativos contra el modelo 2 en terminos de coeficiente de determinacion.
Entonces, para seleccionar el mejor modelo para nuestro analisis no vasta con observar el R2 (coeficiente de determinacion), sino que debemos tener en cuenta todos los indicadores.
Tomando en consideracion que, a partir de la incorporacion de un tercer parametro, el R2 ajustado de nuestro modelo no mejora sustancialmente podemos prescindir de observar los estimadores de AIC y Cp de Mallows. Asi la decision que deberiamos tomar es la de quedarnos con el modelo 2 que incluye las variables CLEAR y LIC.
3.
Qué información da el coeficiente de determinación?
Tomando en consideracion las variables disponibles en nuestro data set, la variable que mas explica la tasa de homicidios cada 100.000 habitantes es el porcentaje de homicidios resueltos con arresto de los responsables en donde el 94% de la variabilidad de la variable dependiente es explicada solo por esta variable.
Al introducir la variable “Cantidad de licencias de portacion de armas cada 100.000 habitantes” podemos mejorar las explicacion en un 5%. Sin embargo, el resto de las variables otorgan muy poco explicacion al modelo predictivo. Entonces, lo que el coeficiente de determinacion nos dice es cuanto varia “Y” en funcion de las “Xs” que introduzcamos en nuestro modelo. Este coeficiente lo podemos identificar como R2 y es la proporcion de la varianza de la variable dependiente sobre las dependientes. Siendo R2 un estimador sesgado, el R2 adjusted corrige los errores que produce la incorporacion de multiples variables explicativas.
Tomando esto en cuenta, a continuacion creamos un nuevo modelo con las variables que mejor explican la variacion de “Y” segun “X”.
#Creamos el nuevo modelo de regresion a partir de los resultados del stepwise
modelodetroitimproved<- lm( H ~ CLEAR + LIC, data = detroit )
modelodetroitimproved##
## Call:
## lm(formula = H ~ CLEAR + LIC, data = detroit)
##
## Coefficients:
## (Intercept) CLEAR LIC
## 103.65554 -1.05747 0.014144.
Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámeros?
Los supuestos necesarios para definir la prueba inferencial son los de:
Linealidad.
Variancia constante u homoecedasticidad.
Distribucion normal de los errores/residuos.
Independencia de los residuos.
Lo que esperariamos es que al graficar los puntos de la variable independiente con respecto a la dependiente los mismos presenten un comportamiento lineal ya sea positivo o negativo. A su vez,uno esperaria que los residuos tengan una variancia constante mientras que se distribuyen normalmente. Finalmente, buscaremos corroborar que los residuos constituyan una variable aleatoria o en otras palabras, que sean independientes.
Para corroborar estos supuestos disponemos de herramientras graficas como matematicas.
#Testeo de normalidad
ols_plot_resid_qq(modelodetroitimproved)
A partir del grafico de los quantiles teoricos vs los quantiles de la muestra podemos trazar una linea recta y agrupar a las observaciones en torno a ella. En nuestro caso diriamos que se cumple el supuesto de normalidad. En caso de querer corroborarlo con un modelo matematico utilizariamos la siguiente funcion.
#Test de Shapiro-Wilk para normalidad
shapiro.test(residuals(modelodetroitimproved))##
## Shapiro-Wilk normality test
##
## data: residuals(modelodetroitimproved)
## W = 0.90239, p-value = 0.1442La hipotesis nula del test de Shapiro-Wilk dice que la distribucion de los residuos es normal. Al ser el P valor mayor a 0.05, no rechazamos la hipotesis nula. Entonces diriamos que la distribucion de los valores de los residuos de nuestro modelo es normal.
#Testeo de homoecedasticidad visual
plot(fitted(modelodetroitimproved), residuals(modelodetroitimproved), xlab = "Fitted", ylab="Residuals") + abline(h = 0) + plot(fitted(modelodetroitimproved), abs(residuals(modelodetroitimproved)), xlab = "Fitted", ylab = "|Residuals|")
## integer(0)#Testeo de homoecedasticidad
#Calculo
ols_test_breusch_pagan(modelodetroitimproved)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## -----------------------------
## Response : H
## Variables: fitted values of H
##
## Test Summary
## ----------------------------
## DF = 1
## Chi2 = 0.6554146
## Prob > Chi2 = 0.4181835El test de Breusch-Pagan es utilizado para dar cuenta de la heteroecedasticidad de un modelo de regresion lineal. En nuestro modelo este test arroja un Chi cuadrado de 0.65, lo cual permite aceptar la hipotesis nula y decir que hay homoecedasticidad.
El resto de los supuestos son corroborados en los puntos 6 y 7.
5.
Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
Para evaluar la bondad del ajuste del modelo debemos tomar en consideracion los componentes de la ecuacion lineal. En otras palabras esto significa observar la variancia de los betas de la ecuacion y las variancias de los errores.
La primera medida de bondad de ajuste que podemos tomar es el R2 el cual mide que proporcion de la variacion de la variable dependiente es explicada por la independiente. La segunda medida de bondad de ajuste seria el error standard de la regresion.
“It is useful to have some measure of how well the model fits the data. One common choice is R2 , the so-called coefficient of determination or percentage of variance explained(…)An alternative measure of fit is This quantity is directly related to the standard errors of estimates of ! and predictions. The advantage is that is measured in the units of the response and so may be directly interpreted in the context of the particular dataset. This may also be a disadvantage in that one must understand whether the practical significance of this measure whereas R2 , being unitless, is easy to understand.”(Faraway.J 2009)
Tomando en cuenta estos criterios, debajo observamos que las diferencias en el R2 y el error standard del modelo completo y el modelo mejorado no son signifcativas. Mientras el primer modelo tiene un R2=0.998 y un error standar = 1.309, el modelo mejorado tiene un R2=0.985 y un error estandard= 1.839. Si bien el R2 del modelo mejorado es algunos decimales menos y su error es apenas un poco mayor, dichas diferencias no son importantes.
Uno podria pensar entonces que tanto el primer modelo como el segundo podrian ser buenos. Sin embargo, debido al principio de simplicidad el modelo mejorado es superior al modelo original reduciendo la explicacion a solo dos parametros.
#Primer modelo
summary(modelodetroit)##
## Call:
## lm(formula = H ~ FTP + UNEMP + M + LIC + GR + CLEAR + W + NMAN +
## G + HE + W, data = detroit)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.12195 0.08714 0.55769 -1.29578 1.04816 -0.22445 0.13999 0.12126
## 9 10 11 12 13
## -0.42310 -0.14609 0.15417 0.01634 0.08662
## 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.044e+02 1.405e+02 0.743 0.535
## FTP 3.115e-02 5.166e-02 0.603 0.608
## UNEMP 1.137e+00 1.051e+00 1.082 0.392
## M 1.666e-02 5.632e-02 0.296 0.795
## LIC 1.619e-02 5.610e-03 2.886 0.102
## GR 1.093e-03 4.262e-03 0.256 0.822
## CLEAR -3.643e-01 1.910e-01 -1.908 0.197
## W -1.361e-04 1.733e-04 -0.785 0.514
## NMAN -2.531e-02 1.149e-01 -0.220 0.846
## G -1.300e-01 1.035e-01 -1.257 0.336
## HE 4.730e+00 3.176e+00 1.489 0.275
##
## Residual standard error: 1.309 on 2 degrees of freedom
## Multiple R-squared: 0.9989, Adjusted R-squared: 0.9936
## F-statistic: 187.8 on 10 and 2 DF, p-value: 0.005308#Modelo mejorado
summary(modelodetroitimproved)##
## 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-106.
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.
Los valores atipicos de un modelo de regresion pueden ser observados tanto en las variables independientes como en la dependiente. Estos valores atipicos se conocen como outliers de los residuos en las variables dependiente y como puntos de influencia en la variable independientes. El problema con estos residuos reside en que su influencia sobre las medias puede cambiar los resultados de los estimadores de un modelo. Es por ello que es importante poder identificarlos y someterlos a pruebas que den cuenta de dicha influencia.
#Grafico para residuos studentizados
ols_plot_resid_stud_fit(modelodetroitimproved)
En el grafico superior podemos observar como uno de los casos se encuentra por fuera de los dos devios standard con respecto a la media. Entonces, el caso 2 que estaria perturbando nuestro modelo.
7.
Analizar la colinealidad de las variables predictoras presentes en la ecuación.
Si comparamos nuestro primer modelo vs. el modelo mejorado podemos observar con claridad las diferencias en la colinealidad. El problema de multicolinealidad se plantea de manera tal que las variables predictoras se correlacionan fuertemente entre ellas a punto tal que son casi perfectas y por ende, explican lo mismo.
Siguiendo con nuestro primer modelo y nuestro modelo mejorado observamos que:
#Testeo de multicolinealidad primer modelo
ols_vif_tol(modelodetroit)## Variables Tolerance VIF
## 1 FTP 0.024424080 40.94320
## 2 UNEMP 0.023224490 43.05800
## 3 M 0.018141316 55.12279
## 4 LIC 0.045316150 22.06719
## 5 GR 0.081268646 12.30487
## 6 CLEAR 0.024436084 40.92309
## 7 W 0.001171905 853.31131
## 8 NMAN 0.001205293 829.67373
## 9 G 0.009728401 102.79181
## 10 HE 0.015155411 65.98304#Testeo de multicolinealidad modelo mejorado
ols_vif_tol(modelodetroitimproved)## Variables Tolerance VIF
## 1 CLEAR 0.6921615 1.44475
## 2 LIC 0.6921615 1.44475Tanto en el modelo original como en el modelo mejorado, no hay multicolineadidad entre las variables predictoras. Debido a que los valores de la tolerancia son menores a 10 y/o el VIF es mayor a 0,10 podemos decir que las variables introducidas en los modelos no se correlacionan entre si.
8.
Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso
Cuando se trata de analizar la presencia de observaciones atipicas o casos influyentes un indicador posible a utilizar es la distancia de Cook. Este indicador provee informacion tanto en X como en Y, lo cual significa que nos da aquellos individuos que demuestran alguna anomalia en la variable dependiente tanto como en la independiente.
En cambio, los Dffits se utilizan para identificar solo los puntos de influencia. Los Dfbetas comparan el estimador del parametro incluyendo el caso y excluyendo el caso influyente. El dffits hacen lo mismo que el dfbetas pero con el predictor. El caso que se saca es el mas extremo.
#Grafico de Cook para saber en que casos hay un outlier en la variable tanto dependiente como independiente
ols_plot_cooksd_bar(modelodetroitimproved)
Al igual que lo observamos en el grafico de residuos, se detecta que el caso 2 se halla por fuera del umbral aceptable de Cook=0.308. El valor del caso 2 incluso duplica y casi triplica el tamano de los valores de las observaciones 8,11 y 12. Esto puede corroborarse mediante el calculo de Cook para los puntos de influencia. Debajo se observa que la distancia de Cook para el caso 2 es de 0.34.
#Calculo de los puntos de influencia. Mayor el numero, mayor la probabilidad de ser un punto de influencia
cooks.distance(modelodetroitimproved)## 1 2 3 4 5 6
## 0.042696165 0.349942740 0.105916409 0.009214249 0.034400845 0.058659257
## 7 8 9 10 11 12
## 0.063022059 0.126814751 0.005132083 0.007362778 0.157195255 0.140192490
## 13
## 0.076944289
## attr(,"label")
## [1] "H"
## attr(,"format.spss")
## [1] "F9.2"
## attr(,"display_width")
## [1] 9A partir del analisis de Cook podemos observar que el caso dos es un outlier tanto para las variables independientes como la dependiente. Sin embargo el comportamiento de este caso puede ser distinto segun el parametro que incluyamos en el modelo. En caso de que necesitemos analizar la influencia de los casos sobre las variables predictoras y predichas por separado, debemos analizar los Dffits y los Dfbetas.
#Visualizacion de los dfitts. Es la influencia de un caso sobre la variable dependiente
ols_plot_dffits(modelodetroitimproved)
Si observamos mas detenidamente los casos atipicos sobre la variable dependientes, vemos que el Dffit muestra al caso 2 por encima del umbral. Esto nos da la idea de que si bien es un outlier dentro del modelo, el caso podria no comportarse de la misma manera con respecto a las variables independientes. Es por ello que resulta pertinente analizar los outliers sobre las variables independientes.
#Visualizacion de los dfbetas. Es la influencia de un caso sobre la variable independiente
ols_plot_dfbetas(modelodetroitimproved)
Este caso presenta un comportamiento similar en la variable predictora LIC (“Cantidad de licencias de portacion de armas cada 100.000 habitantes”) mientras que en el resto se halla por dentro de los umbrales aceptables tanto en esa variable predictora como en el resto. Por lo tanto, resultaria pertinente excluir el caso 2 del modelo para ajustarlo y ver asi si existen o no diferencias.
#Excluimos caso dos de la base
detroitclean<- detroit[-2,]
#Corremos nuevamente el modelo mejorado ahora con el ajuste excluyendo el caso 2
modelodetroitimproved2<- lm(H ~ CLEAR + LIC, data = detroitclean)
summary(modelodetroitimproved2)##
## Call:
## lm(formula = H ~ CLEAR + LIC, data = detroitclean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7013 -0.6165 0.5313 0.8390 1.5608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 105.116272 3.905333 26.916 6.53e-10 ***
## CLEAR -1.061932 0.040445 -26.256 8.14e-10 ***
## LIC 0.012698 0.001711 7.421 4.01e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.474 on 9 degrees of freedom
## Multiple R-squared: 0.9933, Adjusted R-squared: 0.9919
## F-statistic: 671.1 on 2 and 9 DF, p-value: 1.607e-10En comparacion con el modelo mejorado con el caso 2, el modelo actual presenta una leve mejoria con respecto al R2 y el error residual. Estos ajustes podrian continuarse pero uno correria el riesgo de sobreajustar el modelo transformandolo en util solo para predicciones en nuestra muestra pero no en la poblacion.
EJERCIO 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.
#Levantamos la base
bajopeso<- read_sav("LOWBWT.sav")1.
Calcular el riesgo relativo y los odds ratio de la variable dependiente con todas las variables dicotómicas. Analizar los resultados.
El riesgo relativo (RR) y el oddsratio corresponden a medidas que han sido utilizadas tradicionalmente en epidemiología pero luego ampliadas a otras ciencias incluyendo las ciencias sociales. El concepto de riesgo relativo refiere a la probabilidad de que un grupo en particular pueda sufrir algún riesgo en contraposición de otro grupo.
En nuestro caso y a partir de la información proveniente de nuestro set de datos calcularemos el riesgo relativo (RR) de tener un hijo con bajo peso al nacer según estén presentes o no los distintos factores de riesgo considerados. Para ello vamos a tener en cuenta que el RR de tener un hijo con bajo peso al nacer para el factor T se define como:
RR = (P (BAJOPESO |T+ )/ P (BAJOPESO |T−))
donde BAJOPESO significa que el sujeto padece la enfermedad, T+ significa que el posible factor de riesgo está presente y T− que está ausente.
A continuación calculamos los riesgos relativos para todas las variables dependientes:
# Creamos las tablas de doble entrada para cada una de las variables dicotomicas y asi luego calcular los riesgos relativos y odds ratio.#CREAMOS TABLAS CON FRENCUENCIAS PARA LAS VARIABLES BAJO PESO X FUMO DURENTE EL EMBARAZO
towaylowandsmoke<- table(bajopeso$SMOKE, bajopeso$LOW )
towaylowandsmoke##
## Peso aceptable Bajo peso
## No fumo 86 29
## Fumo 44 30#TABLA PORCENTUALIZADA
towaylowandsmoke2<- prop.table(towaylowandsmoke,2)
towaylowandsmoke2##
## Peso aceptable Bajo peso
## No fumo 0.6615385 0.4915254
## Fumo 0.3384615 0.5084746El 51% de los niños con bajo peso al nacer son de madres fumadoras, también 33% de los niños que nacieron con un peso aceptable son de madres fumadoras. Para calcular el riesgo relativo, debemos determinar cual es la probabilidad de tener hijos con bajo peso al nacer si la madre fuma y la probabilidad de tener hijos si la madre no fuma.
#riesgo relativo
#Invertimos las probabilidades para calcularlas por fila
towaylowandsmoke2<- prop.table(towaylowandsmoke,1)
towaylowandsmoke2##
## Peso aceptable Bajo peso
## No fumo 0.7478261 0.2521739
## Fumo 0.5945946 0.4054054#Calculamos la probabilidad de tener un hijo con bajo peso al nacer cuando la madre fuma sobre la probabilidad de tener un hijo con bajo peso al nacer cuando la madre no fuma.
rrbxf<- towaylowandsmoke2[2,2]/towaylowandsmoke2[1,2]
rrbxf## [1] 1.607642Entonces viendo el riesgo relativo podemos decir que el riesgo de tener un hijo con bajo peso al nacer es 1.61 mas alta cuando la madre es fumadora que entre las no fumadoras.
El oddsratio en cambio mide la asociacion entre dos eventos A y B, medida que el riesgo relativo no arroja ya que no da cuenta de asociacion. A modo de ejemplo, si volvemos a la relacion entre bajo peso al nacer y fumadoras, si consideramos que el suceso A es bajo peso, entones el contrario Ac es peso aceptable. Asimismo, si el suceso B es fumar, el Bc seria no fumar. Entonces, podriamos calcular el Oddsratio mediante.
OR (BAJOPESO, T) = (P(BAJOPESO | T+)/P(PESOACEPTABLE | T+))/(P(BAJOPESO | T−)/P(PESOACEPTABLE | T−))
#Calculo para el oddsratio
ORbxf<- (towaylowandsmoke2[2,2]/towaylowandsmoke2[2,1])/(towaylowandsmoke2[1,2]/towaylowandsmoke2[1,1])
ORbxf## [1] 2.021944Viendo el odds ratio podemos decir que hay asociacion entre las variables. El efecto de la variable SMOKE sobre el bajo peso al nacer influye de manera tal que se duplica (2.02) la probabilidad de tener un hijo con bajo peso al nacer cuando la mujer fuma.
#CREAMOS LAS TABLAS PARA LAS FRECUENCIAS DE LAS VARIABLES bajo peso x hipertension
towaylowandht<- table(bajopeso$HT, bajopeso$LOW)
towaylowandht##
## Peso aceptable Bajo peso
## No 125 52
## Yes 5 7#TABLA PORCENTUALIZADA
towaylowandht2<- prop.table(towaylowandht,2)
towaylowandht2##
## Peso aceptable Bajo peso
## No 0.96153846 0.88135593
## Yes 0.03846154 0.11864407El 11% de los ninos con bajo peso al nacer son de madres hipertensas, tambien un 3% de los ninos con peso aceptable al nacer son de madres hipertensas.
#Riesgo relativo
towaylowandht2<- prop.table(towaylowandht,1)
rrbxht<- towaylowandht2[2,2]/towaylowandht2[1,2]
rrbxht## [1] 1.985577El riesgo relativo de tener un hijo con bajo peso al nacer cuando la madre tiene hipertension es 1.98 mas alta que cuando la madre no tiene hipertension.
#Odds ratio
ORbxht<- (towaylowandht2[2,2]/towaylowandht2[2,1])/(towaylowandht2[1,2]/towaylowandht2[1,1])
ORbxht## [1] 3.365385Hay una asociacion entre la hipertension y el bajo peso al nacer. De modo tal que la presencia de hipertension triplica (3.37) la probabilidad de presencia de bajo peso en el nino al nacer.
#BAJO PESO X IRRITBALIDAD UTERINA
towaylowandui<- table(bajopeso$UI, bajopeso$LOW )
towaylowandui##
## Peso aceptable Bajo peso
## No 116 45
## Yes 14 14#Tabla porcentualizada
towaylowandui2 <- prop.table(towaylowandui,2)
towaylowandui2##
## Peso aceptable Bajo peso
## No 0.8923077 0.7627119
## Yes 0.1076923 0.2372881El 24% de los ninos con bajo peso al nacer son de madres que sufren irritabilidad urinaria, tambien 11% de los ninos con peso aceptable al nacer son de madres que no sufren irritabilidad urinaria.
#Riesgo relativo
towaylowandui2 <- prop.table(towaylowandui,1)
rrbxui<- towaylowandui2[2,2]/towaylowandui2[1,2]
rrbxui## [1] 1.788889El riesgo de tener hijos con bajo peso al nacer es casi dos veces (1.8) mayor para las madres con problemas de irritabilidad urinaria.
#Odds ratio
ORbxui<- (towaylowandui2[2,2]/towaylowandui2[2,1])/(towaylowandui2[1,2]/towaylowandui2[1,1])
ORbxui## [1] 2.577778Existe entonces una asociacion entre la variable irritabilidad urinaria y el bajo peso al nacer. Viendo el odds ratio podemos decir que el efecto de la variable UI sobre el peso del bebe influye de manera tal que se duplica (2.58) cuando la madre tiene irritabilidad urinaria.
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?
El odds ratio es una medida de asociacion que es ampliamente usada en epidemiologia y tiene por objeto determinar la probabilidad de ocurrencia de un fenomeno presente en la variables independientes. Esto quiere decir la probabilidad de que x=1 o x=0.
Este principio es utilizado en los modelos de regresion logistica ya que permite determinar la probabilidad de ocurrencia de la variable “Y” en tanto “X”. Lo que nos dice el Odds ratio es en cuanto la variable independiente influye sobre la variable dependiente.
A partir de la noción de Odd podemos pensar la forma de mostrar una relación algebraica que nos permita indicar la probabilidad de respuesta -afirmativa- que tiene un entrevistado escogido al azar. En su lógica el Odds ratio se expresa de la siguiente manera:
Podemos entonces definir el logit, simplemente como el logaritmo del Odd:
El modelo de regresion logistica parte del mismo razonamiento y se expresa como el logaritmo natural de la probabilidad de ocurrencia de un fenomeno sobre la probabilidad de ocurrencia menos uno. Entonces, la información que suministra el odds es sobre la fuerza de la probabilidad ocurrencia o no de un evento en el cálculo de la regresión logística.
A continuación se calculan a modo de ejemplo las variables del punto anterior pero mediante el método de regresión logística simple:
3.
Calcule los odds ratio de cada una de las variables predictoras con la variable dependiente? Comentar.
En nuestro set de variable nos encontramos que algunas de las predictoras no son dictomicas, esto implica un tratamiento sobre las mismas para poder calcular los oddsratio. Esta medida de asociacion se puede calcular solo a partir de la ocurrencia o no de los fenomenos. Por lo tanto, es necesario transformar las variables de numericas, ordinales y nominales con mas de dos alternativas, en variables dicotomicas.
En el caso de nuestro ejemplo, las variables que se dicotomizaran son:
EDAD RAZA CANTIDAD DE EMBARAZOS PREMATUROS CONSULTAS MEDICAS
Para dicotomizar estas variables observaremos la distribucion de las mismas para observar tendencias. Asimismo, otras variables como las categoricas, se dicotomizaran a partir de una aproximacion teorica.
#Observamos la distribucion para dicotomizar las variables cuantitativas. Corregir y crear variables dummy como simulacion
#Histograma edad
hist.scott(bajopeso$AGE)
#Histograma raza
bajopeso$RACE2<- as.numeric(bajopeso$RACE)
hist(bajopeso$RACE2)
#Histograma embarazon prematuros previos
hist(bajopeso$PTL)
#Histograma conusltas medicas
hist(bajopeso$FTV)
#Histograma peso al nacer
hist(bajopeso$BWT)
A partir de la distribucion de frecuencias, las variables se dicotomizaran de la siguiente manera.
EDAD: Mayor de 25 anos y Menor de 25 anos RAZA: Blaco y otras razas. EMBARAZOS PREMATUROS: Sin embarazos prematuros y con embarazos prematuros. CONSULTAS MEDICAS: Sin consultas medicas y con consultas medias.
#Creamos nuevas variables dicotomizandos tentiendo en cuenta las distribuciones observadas en los histogramas
bajopeso2<- bajopeso %>% mutate(AGE2 = case_when(AGE < 25 ~ "Menor de 25",
AGE>=25 ~ "Mayor de 25"),
RACE2 = case_when(as.numeric(RACE) == 1 ~ "Blanco",
as.numeric(RACE) >= 2 ~ "Otra raza"),
PTL2= case_when(PTL ==0 ~ "Sin embarazo prematuro",
PTL >= 1 ~ "Con embarazo prematuro"),
FTV2= case_when(FTV == 0 ~ "Sin consultas",
FTV>= 1~ "Con consultas"))#Creamos las tablas de doble entrada para la edad recodificadas y el bajo peso al nacer y corremos el analisis del oddsratio.
#Tabla de doble entradada
towaylowandage<- table(bajopeso2$AGE2, bajopeso2$LOW)
#Tabla con proporciones por fila
towaylowandage2<- prop.table(towaylowandage,1)
#Calculo para el oddsratio
ORbxa<- (towaylowandage2[1,2]/towaylowandage2[1,1])/(towaylowandage2[2,2]/towaylowandage2[2,1])
ORbxa## [1] 0.76Para las variables edad de la madre y el bajo peso al nacer no hay asociacion. El oddsratio que arroja es menor a 1 (0.76), entonces se puede decir que no hay asociacion entre las variables.
#Creamos las tablas de doble entrada para la raza recodificadas y el bajo peso al nacer
towaylowandrace<- table(bajopeso2$RACE2, bajopeso2$LOW)
#Tabla con proporciones por fila
towaylowanrace2<- prop.table(towaylowandrace,1)
#Calculo para el oddsratio
ORbxr<- (towaylowanrace2[1,2]/towaylowanrace2[1,1])/(towaylowanrace2[2,2]/towaylowanrace2[2,1])
ORbxr## [1] 0.4988584Para las variables raza de la madre y el bajo peso al nacer no hay asociacion. El oddsratio que arroja es de 0.50.
#Creamos las tablas de doble entrada para embarazo prematuro recodificadas y el bajo peso al nacer
towaylowandprematuros<- table(bajopeso2$PTL2, bajopeso2$LOW)
#Tabla con proporciones por fila
towaylowandprematuros2<- prop.table(towaylowandprematuros,1)
#Calculo para el oddsratio
ORbxpr<- (towaylowandprematuros2[1,2]/towaylowandprematuros2[1,1])/(towaylowandprematuros2[2,2]/towaylowandprematuros2[2,1])
ORbxpr## [1] 4.317073Para las variables embarazo prematuros previos y hijos con bajo peso al nacer se observa una fuerte asociacion. Las madres que tuvieron un embarazo prematuro anterior tienen un 4.4 de probabilidades de tener un hijo con bajo peso.
#Creamos las tablas de doble entrada para consultas al medico recodificadas y el bajo peso al nacer
towaylowandmedical<- table(bajopeso2$FTV2, bajopeso2$LOW)
#Tabla con proporciones por fila
towaylowandmedical2<- prop.table(towaylowandmedical,1)
#Calculo para el oddsratio
ORbxcons<- (towaylowandmedical2[1,2]/towaylowandmedical2[1,1])/(towaylowandmedical2[2,2]/towaylowandmedical2[2,1])
ORbxcons## [1] 0.6195286Para las variables consultas previas al medico y el bajo peso al nacer no hay asociacion. El oddsratio que arroja es de 0.62.
Podriamos calcular el Odds ratio sin necesidad de recategorizar las variables en dicotomicas, para ello debemos utilizar UN MODELO DE REGRESION LOGISTICA. Se debe asignar el tipo de variable con la que queremos que el modelo tome nuestras variables y luego utilizar la funcion de glm para determinar la asociacion.
# Regresion logistica bajo peso al nacer segun edad.
lowxage<- glm(LOW ~ AGE, family = "binomial", data = bajopeso)
summary(lowxage)##
## Call:
## glm(formula = LOW ~ AGE, family = "binomial", data = bajopeso)
##
## 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: 4#El incremento usado para AGE es distinto ya que solo usamos incremento ! para variables dicotomicas o cualitativas. El caso de edad usamos un incremento de 15 para que nos compare aquellos grupos de edad que difieran en 15 unidades.
gg<- or_glm(data = bajopeso, model = lowxage, incr = list(AGE=15))
gg## predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1 AGE 0.464 0.178 1.144 15Utilizando la funcion de or_glm con un incremento de 15 para la edad observamos que no hay asociacion significativa entre la variable edad y peso debido a que el intervalo del confianza inferior (0.494) se encuentra por debajo de 1.
# Regresion logistica bajo peso al nacer segun embarazo prematuro previo
lowxptl<- glm(LOW ~ PTL, family = "binomial", data = bajopeso)
gg3<- or_glm(data = bajopeso2, model = lowxptl, incr = list(PTL=1))
gg3## predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1 PTL 2.23 1.218 4.284 1Utilizando la funcion de or_glm con un incremento de 1 para los embarazos prematuros previos observamos que hay asociacion significativa entre la variable embarazos prematuros y bajo peso al nacer debido a que el intervalo del confianza inferior (1.218) se encuentra por encima de 1.
# Regresion logistica bajo peso al nacer segun raza
lowxrace<- glm(LOW ~ factor(RACE), family = "binomial", data = bajopeso)
gg2<- or_glm(data = bajopeso, model = lowxrace, incr = list(RACE=1))
gg2 ## predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1 factor(RACE)Black 2.328 0.926 5.775 Indicator variable
## 2 factor(RACE)Other 1.889 0.957 3.758 Indicator variablePara el caso de la variable RACE, observamos que el intervalo de confianza del limite inferior se halla por debajo de 1 y por lo tanto, podemos decir que no hay una asociacion entre el bajo peso al nacer y la raza de la madre.
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.)
El modelo de regresion logistica multiple nos permitira determinar cual es el impacto de las variables independientes sobre la variable dependiente.
Cuando nuestro set de variables explicativas es tan amplio podemos encontrarnos con algunas dificultades para definir aquellas que seran ingresadas en el modelo. Muchas veces las variables independientes no impactan significativamente en nuestra variable explicada, como asi tambien puede que suceda que entre las variables explicativas se den fenomenos de multicolinealidad. Por ello, el modo de buscar el modelo mas optimo removiendo aquellas variables que no sumen a la explicacion y que presenten multicolinealidad es mediante el metodo de Akaike Information Critera (AIC).
El metodo de AIC funciona de manera similar al R2 ajustado siguiendo la logica de penalizar por agregar mas variables al modelo. Es una medida que funciona siempre y cuando en referencia a otras, ya que un AIC no nos dice nada por si mismo. Solamente podemos determinar cual es el modelo mas optimo de nuestro set de variables disponibles a partir de seleccionar aquel que tenga el valor de AIC mas bajo.
El metodo StepAIC permite tres tipos de tecnicas de seleccion, las cuales son:
Backward
Forward
Both
La explicacion de estas tecnicas de seleccion se halla en el punto 2 de regresion lineal multiple siendo que la lógica para la selección es similiar. Para nuestro objetivo, vamos a tomar la tecnica de seleccion “both” que combina el metodo de seleccion “backward” y “forward”.
#Primero corremos un modelo binomial sin predictores. Por eso se agrega ~1 en lugar de los predictores.
logistica <- glm(LOW~1, family = "binomial", data = bajopeso )
#Luego cargamos el modelo a evaluar. EL objeto creado en la linea anterior sera usado como modelo inicial asignando como modelo maximo el set de predictores y como minimo 1. Asimismo, la probabilidad de entrada de 0.10 y de salida 0.15.
modelostep<- stepAIC(logistica, scope= list( upper= ~ AGE + LWT + factor(RACE) + factor(SMOKE) + PTL + factor(HT) + factor(UI) + FTV, lower=~1), direction = c("both"), trace = T, pent=0.10, prem=0.15)## Start: AIC=236.67
## LOW ~ 1
##
## Df Deviance AIC
## + PTL 1 227.89 231.89
## + LWT 1 228.69 232.69
## + factor(UI) 1 229.60 233.60
## + factor(SMOKE) 1 229.81 233.81
## + factor(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
## + factor(HT) 1 223.58 229.58
## + AGE 1 224.27 230.27
## + factor(RACE) 2 222.53 230.53
## + factor(SMOKE) 1 224.78 230.78
## + factor(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
## + factor(HT) 1 215.96 223.96
## + factor(RACE) 2 217.68 227.68
## + factor(SMOKE) 1 220.54 228.54
## + AGE 1 221.05 229.05
## + factor(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 + factor(HT)
##
## Df Deviance AIC
## + factor(RACE) 2 210.85 222.85
## + factor(UI) 1 213.01 223.01
## + factor(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
## - factor(HT) 1 223.41 229.41
## - LWT 1 223.58 229.58
##
## Step: AIC=222.85
## LOW ~ PTL + LWT + factor(HT) + factor(RACE)
##
## Df Deviance AIC
## + factor(SMOKE) 1 204.90 218.90
## + factor(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
## - factor(HT) 1 217.68 227.68
## - LWT 1 218.69 228.69
##
## Step: AIC=218.9
## LOW ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE)
##
## Df Deviance AIC
## + factor(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
## - factor(SMOKE) 1 210.85 222.85
## - factor(RACE) 2 213.15 223.15
## - factor(HT) 1 211.55 223.55
## - LWT 1 211.62 223.62
##
## Step: AIC=217.99
## LOW ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE) +
## factor(UI)
##
## Df Deviance AIC
## <none> 201.99 217.99
## - PTL 1 204.22 218.22
## - factor(UI) 1 204.90 218.90
## + AGE 1 201.43 219.43
## + FTV 1 201.93 219.93
## - factor(SMOKE) 1 207.73 221.73
## - LWT 1 208.11 222.11
## - factor(RACE) 2 210.31 222.31
## - factor(HT) 1 209.46 223.46A partir del analisis de stepAIC podemos observar que el modelo que arroja el valor de AIC mas bajo es
LOW ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE) + factor(UI)
Es en este punto donde vale la pena aclarar que el metodo de regresion logistica utiliza el metodo de maxima verosimilitud para ajustar nuestro modelo a un modelo ideal. A diferencia del metodo de regresion lineal, el metodo de regresion logistica no puede asumir el metodo de los minimos cuadrados para trazar una linea perfecta de asociacion entre las variables. Por lo tanto, se transforma esta ultima en una funcion logaritmica como la que se muestra a continuacion:
Es a partir de este metodo que se introduce, junto al concepto de AIC, el concepto de “Deviance” como un valor que mide la bondad de ajuste de una regresion logistica. Se entiende como la distancia que tenemos al ajuste perfecto del modelo.
Se observa entonces que este modelo seleccionado no solo presenta el menor valor de AIC, sino que tambien las menores deviancias para cada variable explicativa.
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?
Para determinar cuales son los principales factores de riesgo de bajo peso al nacer y la magnitud de su efecto debemos tomar el modelo obtenido y observar su significancia sobre la variable dependiente.
#El modelo obtenido es el de menor aic (AIC=217.09)
modeloreglog2<- glm(LOW ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE) +
factor(UI), family = "binomial", data = bajopeso)
summary(modeloreglog2)##
## Call:
## glm(formula = LOW ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE) +
## factor(UI), family = "binomial", data = bajopeso)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9049 -0.8124 -0.5241 0.9483 2.1812
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.086550 0.951760 -0.091 0.92754
## PTL 0.503215 0.341231 1.475 0.14029
## LWT -0.015905 0.006855 -2.320 0.02033 *
## factor(HT)Yes 1.855042 0.695118 2.669 0.00762 **
## factor(RACE)Black 1.325719 0.522243 2.539 0.01113 *
## factor(RACE)Other 0.897078 0.433881 2.068 0.03868 *
## factor(SMOKE)Fumo 0.938727 0.398717 2.354 0.01855 *
## factor(UI)Yes 0.785698 0.456441 1.721 0.08519 .
## ---
## 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: 201.99 on 181 degrees of freedom
## AIC: 217.99
##
## Number of Fisher Scoring iterations: 4Vemos entonces que en nuestro modelo solo las variable UI (irritabilidad urinaria) y PTL (embarazos prematuros previos) no son significativos. Si bien estas variables podrian ser retiradas del modelo para asi optimizarlo, observaremos a continuacion que el valor de AIC seria mayor.
#Modelo de prueba sin las variables significativas estadisticamente
modeloreglog3<- glm(LOW ~ LWT + factor(HT) + factor(RACE) + factor(SMOKE), family = "binomial", data = bajopeso)
summary(modeloreglog3)##
## Call:
## glm(formula = LOW ~ LWT + factor(HT) + factor(RACE) + factor(SMOKE),
## family = "binomial", data = bajopeso)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7751 -0.8747 -0.5712 0.9634 2.1131
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.352049 0.924438 0.381 0.70333
## LWT -0.017907 0.006799 -2.634 0.00844 **
## factor(HT)Yes 1.749163 0.690820 2.532 0.01134 *
## factor(RACE)Black 1.287662 0.521648 2.468 0.01357 *
## factor(RACE)Other 0.943645 0.423382 2.229 0.02583 *
## factor(SMOKE)Fumo 1.071566 0.387517 2.765 0.00569 **
## ---
## 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: 208.25 on 183 degrees of freedom
## AIC: 220.25
##
## Number of Fisher Scoring iterations: 4A modo de conclusion uno podria decir que detras de las variables de nuestro modelo hay otras variables con las cuales no contamos que podrian explicar mejor el bajo peso al nacer del nino. Se podria decir incluso que en realidad la variable raza (RACE) esta funcionando como proxy a una variable que daria cuenta de las condiciones de vida de la madre. Una madre con problemas alimenticios producto de la marginalidad, no satisfaceria las necesidades nutricionales basicas para el desarrollo del nino en el utero.
En el caso de la magnitud, podemos abordar el problema mediante los oddsratio de cada una de las variables independientes con la variable dependiente. Para nuestro modelo, la maginitud del impacto de los factores de riesgo son los siguientes:
#Oddsratio para nuestro modelo de regresion logistica seleccionado, Para ello debemos introducir los incrementos correspondientes a las variables independientes.
bajopeso<- bajopeso %>% mutate(LOW2 = case_when(LOW == "Bajo peso" ~ 1, LOW == "Peso aceptable" ~ 0)) #recodificamos para asignar valores 1 y 0, presencia y ausencia del evento.
modeloreglog2<- glm(LOW2 ~ PTL + LWT + factor(HT) + factor(RACE) + factor(SMOKE) +
factor(UI), family = "binomial", data = bajopeso)
summary(modeloreglog2)##
## Call:
## glm(formula = LOW2 ~ PTL + LWT + factor(HT) + factor(RACE) +
## factor(SMOKE) + factor(UI), family = "binomial", data = bajopeso)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9049 -0.8124 -0.5241 0.9483 2.1812
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.086550 0.951760 -0.091 0.92754
## PTL 0.503215 0.341231 1.475 0.14029
## LWT -0.015905 0.006855 -2.320 0.02033 *
## factor(HT)Yes 1.855042 0.695118 2.669 0.00762 **
## factor(RACE)Black 1.325719 0.522243 2.539 0.01113 *
## factor(RACE)Other 0.897078 0.433881 2.068 0.03868 *
## factor(SMOKE)Fumo 0.938727 0.398717 2.354 0.01855 *
## factor(UI)Yes 0.785698 0.456441 1.721 0.08519 .
## ---
## 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: 201.99 on 181 degrees of freedom
## AIC: 217.99
##
## Number of Fisher Scoring iterations: 4modeloreglog2oddsratio<- or_glm(data = bajopeso, model = modeloreglog2, incr = list(PTL=1, LWT=10, RACE=1, SMOKE=1, UI=1))
modeloreglog2oddsratio## predictor oddsratio ci_low (2.5) ci_high (97.5) increment
## 1 PTL 1.654 0.855 3.305 1
## 2 LWT 0.853 0.739 0.969 10
## 3 factor(HT)Yes 6.392 1.694 27.265 Indicator variable
## 4 factor(RACE)Black 3.765 1.356 10.681 Indicator variable
## 5 factor(RACE)Other 2.452 1.062 5.877 Indicator variable
## 6 factor(SMOKE)Fumo 2.557 1.185 5.709 Indicator variable
## 7 factor(UI)Yes 2.194 0.888 5.388 Indicator variableSiguiendo con lo expuesto arriba, diriamos que los factores de riesgo con mayor magnitud de impacto en la probabilidad de tener un nino con bajo peso al nacer son; la hipertension (6.4), el tabaquismo (2.6) y ser de raza negra (3.8).
6.
Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámetros?
La regresion logistica no parte de los mismos supuestos de los que parten la regresion lineal y los modelos lineales. Estos ultimos se basan en el metodo de los minimos cuadrados tales como la normalidad y la homoecedasticidad, mientras que la regresion logistica trabaja con todo tipo de relaciones, debido a que aplica una transformacion no linear logaritmica a los Oddsratio predichos. Por lo tanto, los supuestos que se aplican a estos modelos difieren levemente de los supuestos lineales.
Entre los principales supuestos necesarios para la regresion logistica necesitamos que:
- Contemos con variables dicotomicas. Si bien las variables independientes pueden ser o no dicotomicas, la variable explicada debe ser binaria.
*Debido a que la regresion logistica asume que P(Y=1) es la probabilidad de que un evento ocurra, es necesario que la variable explicada se encuentre codificada acorde al resultado deseado.
*El modelo debe ser ajustado correctamente, se debe evitar el overfitting u el underfitting. Solo las variables que optimicen el modelo deben ser incluidas.
*Se requiere que las observaciones sean independientes y que las variables independientes sean independientes entre si.
*De la mano del punto anterior, se requiere que las variables que presenten multicolinealidad sean removidas o, al menos, reducir la multicolinealidad lo mas posible.
*Si bien la linealidad en sentido de la regresion lineal no aplica la regresion logistica, se busca que la linealidad se asuma entre el logaritmo de los odds y las variables independientes. De otro modo, el test subestima la fuerza de la relacion rechazandola.
*Por ultimo, es necesario trabajar con muestras grandes. Debido a que el metodo de maxima verosimilitud es menos “poderoso” que el metodo de los minimos cuadrados, la regresion logistica requiere al menos un minimo de 10 observaciones en las variables independientes.
7.
Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
La bondad del ajuste de un modelo de regresion logistica sirve para comprobar que tan bueno es el ajuste de los valores predichos por el modelo a los valores observados.
Para ello, buscaremos contrastar la hipotesis nula con respecto a la hipotesis alternativa. Nuestra hipotesis nula dice que “los valores predichos por el modelo se ajustan a los valores observados” y nuestro objetivo es no rechazarla.
Los principales modos de analizar la bondad del ajuste del modelo obtenido son mediantes distintos metodos, los cuales se pueden dividir en los siguientes grupos:
*Test basados en patrones de las covariables (Devianza D y Estadistico Chi cuadrado de Pearsons X2).
*Test basados en probabilidades estimadas (Estadistico de Hosmer-Lemeshow C9 y Estadistico de Hosmer-Lemeshow H9).
*Test Score (Test Score de Brown B, Estadistico Score de Stukel Sst y Estadistico de Tsiatis T).
*Test basados en residuos suavizados (Estadisticos de le Cessie y Van Houwelingen Tlc).
*Medidas tipo R2.
*Curva de ROC
A continuacion procedemos a analizar la bondad de ajuste con los metodos basados en patrones de la covariabilidad y mediante la curva de ROC.
Metodos basados en patrones de la covariabilidad
#Test de bondad de ajuste para la devianza.
sum(residuals(modeloreglog2, type = "deviance")^2)## [1] 201.98561-pchisq(sum(residuals(modeloreglog2, type = "deviance")^2),181)## [1] 0.1361586Como detallamos antes, nuestro objetivo es no rechazar la hipotesis nula y por ende, nuestro test de bondad de ajuste para la devianza deberia ser mayor a 0.05. Siendo que P= 0.13 entonces diriamos que aceptamos la hipotesis nula. Nuesro modelo global se ajusta bien a los datos.
#Test de bondad de ajuste para X2
sum(residuals(modeloreglog2, type = "pearson")^2)## [1] 183.82581-pchisq(sum(residuals(modeloreglog2, type = "pearson")^2),181)## [1] 0.4274947Al igual que con el test de bondad de ajuste de la deviance, observamos que el test de X2 arroja la misma conclusion. Siendo que X2 > 0.05, no rechazamos la hipotesis nula y diriamos tambien que nuestro modelo se ajusta bien.
Metodo curva de ROC.
#Calculo de probabilidades
predicmodel<- plogis(predict(modeloreglog2, bajopeso))
#calculo del punto de corte optimo
optCutOff_modelo <- optimalCutoff(bajopeso$LOW, predicmodel)
#Error de clasificacion
errorclasmodelo<- misClassError(bajopeso$LOW, predicmodel, threshold = optCutOff_modelo)
#Curva de ROC
plotROC(bajopeso$LOW2, predicmodel)
El metodo de la CURVA ROC (Receiving operating characteristics curve) nos dice que tan bien se ajusta nuestro modelo al momento de predecir la variable dependiente. Para ello trabaja agrupando las observaciones de nuestro modelo en dos vectores que permiten identificar la sensibilidad y la especificidad del mismo.
La sensibilidad representa el porcentaje de los verdaderos positivos de nuestra prueba, mientras que la especifidad da cuenta del porcentaje los falsos positivos.
#Calculo de sensibilidad del modelo
sensitivity(bajopeso$LOW2, predicmodel, threshold = optCutOff_modelo)## [1] 1En la primera se hallan las observaciones que se predijeron correctamente como positivas de todas las observaciones positivas.
#Calculo de especificidad del modelo
specificity(bajopeso$LOW, predicmodel, threshold = optCutOff_modelo)## [1] 0.005291005La segunda es la proporción de observaciones que se predicen incorrectamente como positivas de todas las observaciones negativas.
La curva ROC muestra el equilibrio entre estos dos vectores y funciona como clasificador de nuestro modelo, cuanto mas diagonal la curva menor precision tendra nuestro modelo. El area bajo la curva (AUC) nos dice que tan preciso sera nuestro estimador.
Si observamos los resultados podemos dar cuenta de que el modelo tiene una alta sensibilidad y una baja especificidad. Si bien puede predecir correctamente como positivas todas las positivas, predice incorrectamente positivas las observaciones negativas. Con un AUC de 74% podemos decir que el caracter del predictor del modelo es regular.
8.
Indicar porcentaje de casos bien predichos por el modelo.
Para saber cual es el porcentaje de casos bien predichos por el modelo es necesario realizar una matriz de clasificacion que agrupe los casos de nuestra muestra de manera dicotomica ordenando los casos que nuestro modelo predice y los reales. De esta forma, un buen modelo seria aquel que tenga la mayor cantidad de verdaderos positivos y verdaderos negativos.
A continuacion realizamos una matriz de clasificacion con nuestro modelo. Como criterio de corte se toman los valores obtenidos en el punto anterior (optcutoff_modelo = 0.37) para confeccionar la matriz con nuestros datos.
#Matriz de clasificacion
cmatrix_modelo<- confusionMatrix(bajopeso$LOW2, predicmodel, threshold = optCutOff_modelo)
cmatrix_modelo## 0 1
## 0 1 0
## 1 129 59Una vez obtenida nuestra matriz podemos identificar los casos falsos positivos y verdaderos positivos, asi podremos saber cual es el porcentaje de casos bien predichos con una simple formula.
#Calculo para saber cuantos casos son bien predichos por el modelo
porc_predmodelo<- (cmatrix_modelo[1,1]+cmatrix_modelo[2,2])/ (cmatrix_modelo[1,1]+cmatrix_modelo[1,2]+ cmatrix_modelo[2,1]+cmatrix_modelo[2,2])
porc_predmodelo## [1] 0.3174603Nuestro modelo predice un 32% de los casos de manera correcta. Por lo tanto, la tasa de clasificación es baja y podría ser mejorada. Para ello, nuestro modelo necesitaría ser ajustado con variables que puedan predecir mejor el bajo peso del niño al nacer.
Bibliografía.
Hilbe Joseph “Practical Guide to Logistic Regression”
Anderson David y Burnham Kenneth “Multimodel Inference: Understanding AIC and BIC in Model Selection”
Bliese Paul " Multilevel modeling in R. A brief introduction to R, the multilevel package and the nlme package".
Cabo Tania Iglesias “Metodos de Bondad de Ajuste en Regresion Logıstica”
Faraway Julian “Linear models with R”
Hosmer David y Lemeshow Stanley “Applied logistic regression”.
Yang Hongwei “The Case for Being Automatic: Introducing the Automatic Linear Modeling (LINEAR) Procedure in SPSS Statistics”
Lopez Abelardo “Estudio de AIC y BIC en la seleccion de modelos de vida con datos censurados”
Reche Jose Luis “Regresión logística. Tratamiento computacional con R.”
Rioja Luis, Llorente Alejandro y Ramirez Beatriz “Regresión Logística: Fundamentos y aplicación a la investigación sociológica”
Santana Angelo “Analisis y tratamiento de datos:Aplicaciones con R”
Triola Mario “Estadística. Decimoprimera edición”
Walpole Ronald, Myers Raymond, Myers Sharon, Ye Keying .“Probabilidad y estadística para ingeniería y ciencias. Novena edición”