Apunte_Estadística_II
Juan Nieve
#install.packages("pacman") #con este paquete levantas todes les paquetes que querés juntes y los lee sin tener que poner library() y si no están te los busca y te los instala y después los lee. Lo +.
pacman::p_load(tidyverse,ggpubr,rstatix,plyr,pastecs,summarytools,RCurl,reshape2,pracma,raster,dplyr,gmodels,DescTools, reshape2, Hmisc, olsrr, fastDummies)I. ANOVA (Análisis de la Varianza)
a. Antes de ANOVA: análisis de dos parámetros
La estadística inferencial implica la formación de conclusiones (o inferencias) sobre un parámetro poblacional. Dos actividades cruciales de la estadística inferencial son la estimación de parámetros poblacionales utilizando intervalos de confianza (IC) y las pruebas de hipótesis sobre parámetros poblacionales (PdH).
En este apartado se verán los métodos para distinguir entre una situación que involucra dos muestras independientes y una que involucra dos muestras dependientes. A partir de eso se realizará una PdH formal de una afirmación sobre dos medias poblacionales diferentes (𝝻1 y 𝝻2) y también una estimación del IC para la diferencia entre dos medias poblacionales.
Pueden darse dos situaciones diferentes:
Una comparación entre dos medias en las que las desviaciones estándar de las dos poblaciones son desconocidas y no se supone que sean iguales
Una comparación entre dos medias en las que las desviaciones estándar de las dos poblaciones son desconocidas pero se suponen iguales (no se analizará el caso “poco realista” en el que se conocen las dos desviaciones estándar poblacionales). Este segundo punto es la base para después comprender la comparación entre más de dos medias, a partir del método conocido como ANOVA (Analisys of the Variance)
Previamente, unas escuetas (?) definiciones:
Dos muestras (n1 y n2) son independientes si los valores muestrales de una población no están relacionados .o de alguna manera naturalmente emparejados o combinados. con los valores muestrales de otra población. Por ejemplo, si las dos muestras tienen diferentes tamaños n sin datos faltantes, deben ser independientes.
Dos muestras son dependientes (o constan de pares relacionados) si los vlaores muestrales se corresponden de alguna manera, donde la correspondencia se basa en una relación inherente. Es decir, cada par de valores muestrales consiste en dos medidas del mismo sujeto -como datos de antes y después-, o cada par de valores muestrales consiste en pares coincidentes, como datos de marido y mujer, y donde la coincidencia se basa en alguna relación significativa. Recordar que dependencia no requiere una relación de causalidad.
Completar con todo el análisis para dos muestras de Triola
Introducción a ANOVA
En el apartado anterior se analizó cómo realizar una PdH para una afirmación que involucraba dos medias poblacionales (𝝻1, 𝝻2) con sus respectivas varianzas (𝞂21, 𝞂22). Como todo análisis inferencial parte de una muestra (en este caso, dos muestras n1 y n2), calculamos las medias muestrales x̄1 y x̄2 y sus desvíos (s21 y s22) con el fin de realizar nuestra PdH.
A grandes rasgos, esta prueba no difirió demasiado de una PdH para una sola población con la única salvedad de que se hizo necesario reconocer previamente si nuestros datos muestrales correspondían a dos muestras pareadas o si eran dos muestras independientes; esto condiciona con qué fórmula debemos trabajar y la forma de calcular nuestro estadístico de prueba. Fuera de eso, el procedimiento fue similar:
- Se analizó si, cumplidos ciertos supuestos, mismas podían ser consideradas iguales (H0: 𝝻1 - 𝝻2 = 0) o no (Ha: 𝝻1 - 𝝻2 ≠ 0 | Ha: 𝝻1 - 𝝻2 > 0 | Ha: 𝝻1 - 𝝻2 < 0).
- Se definió el tipo de prueba (de una cola o dos colas) y el nivel de significancia (α) con el que definimos nuestra zona de rechazo.
- Se identificó y calculó el p-value o estadístico de prueba (t en general, aunque de ser conocida la varianza 𝞂2 e las poblaciones puede tambien utilizarse el z score)
- En función de eso se interpretaron los resultados
Ahora bien, ¿qué onda si tenemos más de dos poblaciones?
Uno podría decir que una forma válida para comparar tres o más medias poblacionales poblaciones (𝝻1, 𝝻2, 𝝻3…𝝻k) sería aplicar el mismo método que antes, nada más que calculando la relación entre las medias poblacionales de a pares, algo como: H0: 𝝻1 = 𝝻2 | 𝝻1 = 𝝻3 | 𝝻2 = 𝝻3 . El problema de esto es que estaríamos realizando tres PdH con un nivel de significancia (α) para cada una que aumentaría la probabilidad de que algunas diferencias resulten significativas debido al azar. Por ejemplo, si tomasemos para nuestra PdH un α = 0,05 e hiciesemos las comparaciones de a pares (por ejemplo entre 𝝻1, 𝝻2 y 𝝻3 ) la probabilidad de encontrar al menos una diferencia significativa por azar sería, en los hechos, del 0,143 (0,053 = 0,05*0,05*0,05 =0,1426…) por la combinatoria de los α de cada PdH. Es decir, aumentamos el riesgo de un error de Tipo I.
Por lo tanto, una prueba estadística basada en todos los datos utilizados simultáneamente, es más estable que la prueba o análisis que parcializa los datos y no los examina todos juntos.
Entonces, muchache, lo que debes usar es, como PdH, el famoso Análisis de la Varianza (ANOVA). Aclaración: el ANOVA no constituye un método o procedimiento único sino que según los diseños y datos disponibles existen diversos modelos de análisis de la varianza:
- ANOVA de un factor
- ANOVA de varios factores
- Diseños factoriales
- ANOVA de la varianza de medidas repetidas
Origen y conceptos básicos de ANOVA
La técnica de análisis de varianza (ANOVA) también conocida como análisis factorial y desarrollada por Fisher en 1930, constituye la herramienta básica para el estudio del efecto de uno o más factores (cada uno con dos o más niveles) sobre la media de una variable continua. Es por lo tanto el test estadístico a emplear cuando se desea comparar las medias de dos o más grupos. Esta técnica puede generalizarse también para estudiar los posibles efectos de los factores sobre la varianza de una variable.
La hipótesis nula de la que parten los diferentes tipos de ANOVA es que la media de la variable estudiada es la misma en los diferentes grupos, en contraposición a la hipótesis alternativa de que al menos dos medias difieren de forma significativa. ANOVA permite comparar múltiples medias, pero lo hace mediante el estudio de las varianzas.
El funcionamiento básico de un ANOVA consiste en calcular la media de cada uno de los grupos para a continuación comparar la varianza de estas medias (varianza explicada por la variable grupo, intervarianza) frente a la varianza promedio dentro de los grupos (la no explicada por la variable grupo, intravarianza). Bajo la hipótesis nula de que las observaciones de los distintos grupos proceden todas la misma población (tienen la misma media y varianza), la varianza ponderada entre grupos será la misma que la varianza promedio dentro de los grupos. Conforme las medias de los grupos estén más alejadas las unas de las otras, la varianza entre medias se incrementará y dejará de ser igual a la varianza promedio dentro de los grupos.1
Resumiendo:
El concepto básico de ANOVA es simple: se compara la varianza que hay entre todas las unidades dentro de cada grupo (Variabilidad Intragrupal VI o SSIntra) con la que hay entre el promedio de los grupos (Variabilidad Intergrupal VE o SSEntre)2. Si el primero es mayor (SSIntra>SSEntre) entonces la variación entre los grupos o muestras no representa una variación real y no se rechaza la H0.
Es importante recalcar que la prueba de hipótesis que se realiza a través del ANOVA está orientada a discernir si las medias poblacionales 𝝻k son iguales (H0). Peeeeeeeeeeeero en caso de rechazar H0 (todas las medias son iguales) y aceptar Ha (“al menos dos de las medias no son iguales”), este tipo de test no brinda información sobre cuáles son las diferencias en sí. Por eso se lo denomina Test Global, y para escudriñar sobre las diferencias entre las medias (for example, qué medias no son iguales) será necesario realizar algún tipo de Test Post Hoc (que se verá más adelante).
El estadístico F y la distribución de probabilidad F de Snedecor
El estadístico estudiado en el ANOVA para la PdH sobre dos o más medias muestrales se calcula a partir de la ratio (la división) entre la varianza de las medias de los grupos (SSEntre) y el promedio de la varianza dentro de los grupos (SSSIntra). Así se obtiene el estadístico F,llamado así debido a que sigue una distribución de probabilidad conocida como “F de Fisher-Snedecor”. Al igual que la t de Student es una familia de curvas cuya curva exacta a usar está determinada por dos grados de libertad (a diferencia chi cuadrado, por ejemplo, ya que justamente se están analizando las distribuciones de dos varianzas simultáneamente).
Si se cumple la hipótesis nula, el estadístico F adquiere el valor de 1 ya que la intervarianza (SSEntre) será igual a la intravarianza (SSIntra). Cuanto más difieran las medias de los grupos mayor será la varianza entre medias en comparación al promedio de la varianza dentro de los grupos (o sea, cuando SSEntre>SSIntra), obteniéndose valores de F superiores a 1 y por lo tanto un menor p-value.
En concreto, si s12 es la la varianza de una muestra de tamaño n1 extraída de una población normal de varianza σ12 y s22 es la la varianza de una muestra de tamaño n2 extraída de una población normal de varianza σ22, y ambas muestras son independientes, el cociente F= (s12/σ12) / (s22/σ22) se distribuye como una variable F de Snedecor.
Grados de libertad de la distribución F: determinación del valor crítico y el estadístico de prueba.
El valor crítico F se encuentra suponiendo una prueba de cola derecha, ya que los valores grandes de F corresponden a diferencias significativas entre las medias (es decir, nos acercan a la zona de rechazo). Para encontrar el valor crítico y el estadístico de prueba, se tomará para el denominador k-1 grados de libertad (donde k es el número de muestras a comparar) y para el numerador con n-k grados de libertad (donde n es el tamaño de muestra y k el número de muestras)
Existen diferentes tipos de ANOVA dependiendo de la si se trata de datos independientes (ANOVA entre sujetos), si son pareados (ANOVA de mediciones repetidas), si comparan la variable cuantitativa dependiente contra los niveles de una única variable explicatoria o factor (ANOVA de una vía) o frente a dos factores (ANOVA de dos vías). Este último puede ser a su vez aditivo o de interacción (los factores son independientes o no lo son). Cada uno de estos tipos de ANOVA tiene una serie de requerimientos.
b. ANOVA de un factor
Pos bien, vamos a arrancar viendo cómo se realiza un ANOVA para un factor en primer lugar.
Consideremos que se seleccionan k muestras aleatorias de tamaño n de una misma variable y,correspondiente a k poblaciones, y queremos saber si hay diferencias estadísticamente significativas en la media (Ȳ) de la variable entre los distintos k grupos definidos.
Para eso, dispondríamos nuestros datos relevados en una tabla similar a esta:
Vamos a un ejemplo más concreto:
Ejemplo 1
Se ha observado que el insomnio es uno de los síntomas asociados a determinados trastornos de ansiedad. Para aliviarlo, un gabinete psicológico propone una terapia para aumentar el número de horas dormidas.
Se conforman dos grupos de sujetos: al primer grupo se le aplica la terapia y al segundo no (grupo placebo), y se observa la cantidad de horas dormidas por los sujetos en tres noches consecutivas.
Se relevaron las siguientes variables (que pueden verse descargando el archivo .sav de aquí):
Primer Noche, Segunda Noche y Tercer Noche corresponden a la cantidad de horas dormidas por los sujetos en las tres noches.
Sexo: Asume los valores 0 para las mujeres y 1 para los hombres.
Edad Asume los valores 1 (menores de 20 años), 2 (entre 20 y 25 años) y 3 (más de 25 años)
La variable Terapia 1 hace referencia a si el sujeto ha recibido terapia contra el insomnio, Terapia 2 refiere a si el sujeto ha recibido terapia contra estados de ansiedad generalizada y Terapia 3 indica si el sujeto ha recibido terapia contra algún tipo de fobia. Todas ellas asumen los valores 1, si el sujeto ha recibido la terapia y o en caso contrario.
library(haven) #para leer archivos en formato .sav de SPSS
urlfile_1 = 'https://github.com/oblitterator/estadistica_II_UNTREF/blob/main/Insomnio_2.sav?raw=true'
Insomnio_2 = read_sav(url(urlfile_1))
print(names(Insomnio_2))## [1] "Clave" "PrimerNoche" "SegundaNoche" "TerceraNoche" "Sexo"
## [6] "Edad" "Terapia1" "Terapia2" "Terapia3"print(nrow(Insomnio_2))## [1] 99Por lo que se puede ver se trata de un dataframe con 10 atributos o columnas (“Clave”, “PrimerNoche”, “SegundaNoche”, “TerceraNoche”, “Sexo”, “Edad”) y 99 observaciones o filas.
El gabinete psicológico desea determinar si existen diferencias en la cantidad media de horas dormidas por los pacientes con insomnio durante la primera noche en los tres grupos de edad.
Lo primero siempre, sempre, always, semper, es arrancar por hacer un análisis exploratorio de nuestros datos (EDA), tirando unos grafiquitos y stats. Así que:
- Realice un gráfico boxplot que permita comparar las horas dormidas durante la primera noche para los 3 grupos etarios considerados. ¿Qué puede concluir?
#hago un boxplot cruzando la información de la columna primera noche y la columna edad.
boxplot(Insomnio_2$PrimerNoche ~ Insomnio_2$Edad)
A simple vista no se ven valores atípicos (u outliers, que traen mucho “ruido” en el ANOVA), y pareciera tratarse de distribuciones asimétricas (no parecen distribuciones normales). Además, parecería haber una diferencia entre las medias de cada grupo. Vamos a hacer una exploración que nos tire unas métricas un poco más precisas:
Insomnio_2 %>%
group_by(Edad) %>% #agrupo por edad
get_summary_stats(PrimerNoche, type = "mean_sd") #pido un resumen de la media y desviación estándar de cada grupo## # A tibble: 3 × 5
## Edad variable n mean sd
## <dbl+lbl> <chr> <dbl> <dbl> <dbl>
## 1 1 [Menor a 20 años] PrimerNoche 33 5.70 1.69
## 2 2 [Entre 20 y 25 años] PrimerNoche 33 7.74 2.37
## 3 3 [Mas de 25 años] PrimerNoche 33 6.88 1.91Ahora que tenemos nuestros datos, antes de pasar a realizar nuestro ANOVA es importante revisar que se cumplan los…
TA TAN TA TANNNNNN
Supuestos para ANOVA de un factor
El ANOVA es una técnica de PdH y por ende es importante que cumpla ciertos requisitos para llevarse a cabo:
Se necesita que los datos tengan una distribución normal
Se supone que las poblaciones que proveen las muestras tienen varianzas iguales (llamada homocedasticidad de varianza σ12 = σ22
Las muestras son independientes (!), o sea no pareadas.
Verificación de supuestos
Continuación Ejemplo 1
b) Realice la verificación de los supuestos para realizar un ANOVA.(Normalidad y Homogeneidad de varianzas)
Test normalidad
Para testear la normalidad realizamos un test de Shapiro-Wilk (o el de Kolmogorov Smirnov) para los grupos seleccionados. Este test es en sí una PdH, donde mi H0 es que los datos provienen de una distribución normal (siendo la inversa la Ha).
#PdH para evaluar la normalidad: Shapiro-Wills
Insomnio_2 %>%
group_by(Edad) %>%
shapiro_test(PrimerNoche) #la función shapiro_test me permite testear si los grupos tienen una distribución normal## # A tibble: 3 × 4
## Edad variable statistic p
## <dbl+lbl> <chr> <dbl> <dbl>
## 1 1 [Menor a 20 años] PrimerNoche 0.973 0.563
## 2 2 [Entre 20 y 25 años] PrimerNoche 0.969 0.458
## 3 3 [Mas de 25 años] PrimerNoche 0.973 0.574Como para cualquier PdH, si mi p-value es “grande” debo rechazar la H0. En este caso veo que para cada uno de los tres grupos, a diferencia de lo que apreciaba en el BoxPlot, el p-value es elevado (comparándolo, por ejemplo, con un nivel de significatividad α = 0,05), por lo que no rechazamos la hipótesis nula y podemos asumir que los datos siguen una distribución normal.
Podemos visualizar esto, también, usando la función qqplot, que compara los cuantiles teóricos de la distribución normal con los cuantiles de la muestra, dando unos gráficos muy bonitos:
ggqqplot(Insomnio_2, "PrimerNoche", facet.by = "Edad")
Test homoscedasticidad
En segundo lugar, vamos a probar que las varianzas de las que provienen las muestras son homogéneas (concepto de homoscedasticidad: del latín homo: misma, cedasticidad: dispersión). Este es un supuesto que debe cumplirse para poder aplicar un test de ANOVA (y que se retomará posteriormente en los modelos de regresión lineal).
Para esto, aplicamos la llamada prueba de Levene. Como reza nuestra amada Wikipedia, el Levene test
“…es una prueba estadística inferencial utilizada para evaluar la igualdad de las varianzas para una variable calculada para dos o más grupos. Algunos procedimientos estadísticos comunes asumen que las varianzas de las poblaciones de las que se extraen diferentes muestras son iguales. La prueba de Levene evalúa este supuesto. Se pone a prueba la hipótesis nula de que las varianzas poblacionales son iguales (llamado homogeneidad de varianza ú homocedasticidad). Si el P-valor resultante de la prueba de Levene es inferior a un cierto nivel de significación (típicamente 0.05), es poco probable que las diferencias obtenidas en las variaciones de la muestra se hayan producido sobre la base de un muestreo aleatorio de una población con varianzas iguales. Por lo tanto, la hipótesis nula de igualdad de varianzas se rechaza y se concluye que hay una diferencia entre las variaciones en la población.”
Matemáticamente:
Continuación Ejemplo 1
Pero bueno, esto es sólo para curioses, ya que lo que vamos a hacer con nuestra muestra es tirarle la función levene_test(), quedando todo así:
Insomnio_2 %>% levene_test(PrimerNoche ~ as.factor(Edad))## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 96 1.18 0.313Esto nos devuelve el resultado del estadístico W (W = 1,18) que demasiado no nos importa, ya que basándonos en el p-value obtenido (p_value = 0,31) podemos decir que se encuentra bastante alejado de nuestra zona de rechazo (pensamos un nivel de significancia en base a un α = 0,05), así que no rechazamos nuestra Hipótesis Nula (las varianzas son homogéneas).
Con nuestros supuestos ya comprobados, ahora sí, podemos ir finalmente al…
Aplicación ANOVA de un factor
Como se dijo más arriba, nuestra intención de un test ANOVA es poder medir las diferencias existentes entre una determinada serie de medias muestrales Ȳi obtenidas de una población 🐮, que cumplen nuestro supuestos de normalidad y homoscedasticidad, y que antes veíamos en esa hermosa tabla de ies, yes, equises, kaes y subíndices de cosas, que ahora se nos presenta en el ejercicio que nos brinda un dataframe con datos de la pobre gente que relevó el gabinete psicológico, gente que sufre de insomnio y están buscando una solución mágica -droga- para ayudarlos.
Pero ojota: nuestra Ȳi no es ninguno de la información dada, sino que tenemos que primero construirla en base a algún factor de comparación que deseemos. Es decir, tenemos que primero agrupar nuestros datos para convertirlos en diferentes Ȳi (según lo que querramos evaluar) y después a eso le hacemos el ANOVA. Por eso el ejercicio dice:
El gabinete psicológico desea determinar si existen diferencias en la cantidad media de horas dormidas por los pacientes con insomnio durante la primera noche (Ȳi) en los tres grupos de edad. Es decir, nuestro factor es el grupo de edad (yi), nuestras observaciones (yn), y de ahí calculamos nuestras Ȳi a comparar en los k grupos existentes (en este caso k=3). Siganme para más tautologías.
Entonces el punto dice:
Continuación Ejemplo 1
c. ¿Puede afirmar que existen diferencias? Utilizar= α = 0,05
Lo que hacemos es utilizar la función aov() para hacer un test ANOVA. La misma tiene adentro dos parámetros que le pasamos: nuestra variable dependiente (“PrimeraNoche”, que corresponde a las horas dormidas en la primera noche por cada individuo, yn), agrupadas por una categoría que armamos en base al grupo etario al que corresponde cada observación y “reordenamos” la tabla para que ahora nos quede como un factor (yi).
Aclaración: en realidad esto nunca lo hicimos porque no tenemos las edades de cada uno de los pacientes sino que ya nos viene agrupada la edad (Edad: 1 =menores de 20 años, 2 = entre 20 y 25 años y 3 = más de 25 años).
Lo que haríamos si lo hiciesemos a mano sería primero armar una tabla donde nuestras x=horas dormidas y nuestra y=Edad (que si no nos la hubieran dado previamente determinado, la tendríamos que haber -atajate esta palabra que voy a inventar- intervalizado). Pero bueno, no lo voy a hacer.
Lo que si voy a hacer es crear una nueva variable del tipo factor (categórica, digamos) diferente a la original (“Edad”), para hacer las operaciones. Esto se debe a que, si bien son intervalos, al estar en números R interpreta que es una variable cuantitativa y básicamente yo quiero calcular una asociación entre diferentes variables cualitativas y cuantitativas.3
Vamos directo al código:
#creo mi variable factorial en base a los grupos etarios
Insomnio_2$Fedad = as.factor(Insomnio_2$Edad)
test_anova = aov(PrimerNoche ~ Fedad, data = Insomnio_2)
summary(test_anova) #nos tira las stats del test## Df Sum Sq Mean Sq F value Pr(>F)
## Fedad 2 69.3 34.66 8.588 0.000371 ***
## Residuals 96 387.5 4.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Y hete allí que nuestro estadístico F de Snedecor es de 8,588 y el p-value(Pr(>F)) = 0,000371. Qué quiere decir esto. Pues como dice San Triola:
“Si la P es baja, la nula debe irse”.
Y como nuestra Hipótesis Nula era que no había diferencias entre los tres rangos etáreos, debemos rechazarla ya que es menor que nuestro nivel de significancia α = 0,05.
Ergo: existen diferencias.
Peeero, no sabemos entre cuáles está la diferencia! Sólo sabemos que nuestras Ȳi son diferentes pero no sabemos entre cuáles se da la diferencia ni cómo, sólo que rechazamos nuestra H0 de que las medias del promedio de noches de cada grupo fue diferente para los pobres insomnes. No nos olvidemos que todo esto lo hacemos por les insomnes.
Así que es momento de aplicar otros tipos de test para ver entre qué Ȳi se da la diferencia: los Test Post Hoc.
c. Test Post-Hoc (la forma de la relación)
ANOVA no dice lo que realmente necesitamos (¿quién nos lo dice igual no?)
¿Qué nos dice?: Los grupos tienen medias diferentes
¿Cuáles grupos?: No se sabe
Pero ahora estamos interesados en diferencias específicas y no en rechazar una hipótesis general. Por lo tanto, es de poca utilidad si no seguimos indagando. Queremos saber entre qué factores se da la relación con más fuerza y si alguno es predominante (o no).
¿Qué se debe hacer para conocer cuáles de las medias difieren?: pues comparaciones múltiples (H0: 𝝻1 = 𝝻2 | 𝝻1 = 𝝻3 | 𝝻2 = 𝝻3 …)
¿Pero cuál fue el problema que vimos antes sobre las PdH mútiples para comparación de medias?
Me pareció que la imagen representaba mejor el problema que lo que yo escribiese.
En fin, como vimos cuando arrancamos con el ANOVA hacer multiples test de student entre medias poblacionales (𝝻) nos aumentaba la posibilidad de nuestro Error de Tipo I. Pero ahora que vemos que hay una relación entre los factores podemos aplicar una serie de cálculos 2.0 para realizar una comparación múltiple sobre la diferencia entre las medidas buscando disminuir los Errores de Tipo I y II que posibles.
Hay varios test que podemos hacer:
Pero todo depende de con qué datos cuentes previamente.
Y sabemos que los elementos que influyen en la comparación entre las medias son:
• La distribución del estadístico
• El nivel de significación fijado
• La variabilidad de cada grupo
• El tamaño muestral de cada grupo
Así que, como dije antes, hay varios métodos de test que podemos aplicar después de que descubrimos con el ANOVA que las medias eran diferentes, pero en general tienen este orden de aparición, según los supuestos que se cumplen:
¿Sobre qué trabajan cada uno de los distintos test para probar entre qué medias se producen diferencias? (cómo es la relación y cuánta)
Los más utilizados
Ahora que rechazamos nuestra H0, queremos afinar nuestra Ha, es decir, perfilar o especificar la misma. Para eso hay diferentes métodos según los supuestos que se cumplen. Estos son los más usuales (música de la Ley y el Orden)
En primera instancia se verán una serie de métodos que lo buscan la mínima distancia significativa. Es decir, parten de utilizar la t de student para la comparación de más de tres medias y operan para reducir el nivel de significancia 𝛂 que se produce en la misma en su comparación. Para esto se utilizan los Test de Bonferroni y/o el de Sidak.
Aquellas que, como el Test de Tukey, no trabajan sobre la significancia sino sobre la misma distribución del desvío estandar resultante de la comparación de las medias (el llamado rango estudentizado) de los grupos, pero ordenadas.
LSD: Tests de Mínima Diferencia Significativa: Bonferroni y Sidak
Mínima Diferencia Significativa: es a través del uso clásico del Test, recomendado para hasta tres grupos. Operan buscando reducir el nivel de significancia α resultante de la comparación de las medias a través de la t de student.
Test de Bonferroni
Test de Bonferroni: Trabaja sobre el nivel de significación, y agrega al alpha el nivel de comparaciones que se desean hacer: α = α/c donde c es el número de comparaciones a realizar. Es decir, fijo un α general de comparación y utilizo para comparar las distintas medias un error “corregido”, según las c comparaciones que quiero realizar.
Entonces, ¿cómo se dan las diferencias entre las medias en nuestro ejercicio de insomnes¡
Continuación Ejemplo 1
d) En caso de existir diferencias verificar para que grupos etarios las medias son distintas de los otros grupos (Tukey y comparar con el BoxPlot). Si bien el de Tukey es el más aplicado, antes muestro cómo son los otros.
# BONFERRONI
Insomnio_2 %>% t_test(PrimerNoche ~ Fedad, p.adjust.method = "bonferroni")## # A tibble: 3 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Primer… 1 2 33 33 -4.03 57.8 1.63e-4 4.89e-4 ***
## 2 Primer… 1 3 33 33 -2.66 63.0 1 e-2 3 e-2 *
## 3 Primer… 2 3 33 33 1.63 61.3 1.09e-1 3.27e-1 nsComo vemos, nos arroja para la comparación de cada grupo si han habido diferencias estadísticamente significativas (nos fijamos el “p.adj” contra nuestro valor de significancia α prefijado).
En este caso los del grupo etarios = 1 y 2 (Menores de 20 y Entre 20 y 25) muestran una diferencia estadísticamente significativa en la cantidad de horas dormidas, lo que se ve en el resultado del P-Value Ajustado (p.adj) = 0.000489. Un valor P-Value mucho menor al α = 0,05.
Para los grupos 1 y 3 (Menores de 20 y Mayores de 25), también se encuentra una diferencia estadísticamente significativa, aunque menor: P-Value Ajustado (p.adj) = 0.03.
Para los grupos 2 y 3 (Entre 20 y 25 y Mayores de 25), no se ve una diferencia estadísticamente significativa en las horas dormidas.
Test de Sidak
El Test de Sidak es similar a Bonferroni pero hace una alternativa un poco menos exigente: α = 1-(1-α)1/c . No tiene sentido mostrarlo (?).
HSD4: Test de Tukey, el mejor
Test de Tukey
El mejor video sobre el tema: https://www.youtube.com/watch?v=hz3fuGVm2js
Analiza cómo se comporta la varianza a través de la utilización de un rango estudentizado. No se va a utilizar la distribución t de Student para definir el término contra el cuál se va a comparar la diferencia de medias.
El Test de Tukey lo que hace es agrupar las medias en diferentes grupos y compara los grupos entre sí. Lo que busca es la diferencia extrema: compara el grupo que tiene el promedio más grande contra el grupo que tiene el promedio más chico (es decir, el máximo contra el mínimo). Esta es la principal diferencia con los LSD.
Esta diferencia ya no es una distribución del tipo t de Student ya que incluye el concepto de orden entre los promedios loque hace que la distribución cambie radicalmente y a ese tipo de distribución se la denomina rango estudentizado. A partir de la misma se calcula una nueva constante: la Diferencia Honestamente Significativa (HSD) que se usará para delimitar el invervalo de comparación para las medias de los grupos.
Es la PREFERIDA DE LES ESTADÍSTIQUES porque es la que mejor permite manejar los Errores de Tipo I y II.
Eso sí: es necesario que cada muestra contenga la misma cantidad de observaciones (n), si no, hay que usar el test de Gabriel o el de Hochberg.
Continuación Ejemplo 1
d) En caso de existir diferencias verificar para que grupos etarios las medias son distintas de los otros grupos (Tukey y comparar con el BoxPlot)
test_anova = aov(PrimerNoche ~ Fedad, data=Insomnio_2) #por las vuelvo a crear mi variable del test anova.
TukeyHSD(test_anova, "Fedad", ordered = TRUE) ## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = PrimerNoche ~ Fedad, data = Insomnio_2)
##
## $Fedad
## diff lwr upr p adj
## 3-1 1.1788718 0.001390959 2.356353 0.0496611
## 2-1 2.0416870 0.864206151 3.219168 0.0002288
## 2-3 0.8628152 -0.314665630 2.040296 0.1941436#tambien puede hacerse: aov(PrimerNoche ~ Fedad, data=Insomnio_2) %>% tukey_hsd()Como vemos, nos remite a las mismas conclusiones, aunque con un “p.adj” más afinado:
En este caso los del grupo etarios = 1 y 2 (Menores de 20 y Entre 20 y 25) muestran una diferencia estadísticamente significativa en la cantidad de horas dormidas, pero con un P-Value Ajustado = 0,000228 (a diferencia del test de Bonferroni que nos había arrojado 0.000489.)
Para los grupos 1 y 3 (Menores de 20 y Mayores de 25), el valor del P-Valuie Ajustado (p.adj = 0.0496) está mucho más cerca de nuestro nivel de significancia α = 0,05, aunque sigue cayendo en la zona de rechazo de la Hipótesis Nula por un pelito.
Para los grupos 2 y 3 (Entre 20 y 25 y Mayores de 25), no se ve tampoco una diferencia estadísticamente significativa en las horas dormidas.
Finalmente, a continuación, podemos apreciar esto mismo pero de forma gráfica:
require(graphics)
plot(TukeyHSD(test_anova, "Fedad"))
tukey_hsd(test_anova)## # A tibble: 3 × 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Fedad 1 2 0 2.04 0.864 3.22 0.000229
## 2 Fedad 1 3 0 1.18 0.00139 2.36 0.0497
## 3 Fedad 2 3 0 -0.863 -2.04 0.315 0.194
## # … with 1 more variable: p.adj.signif <chr>d. Epílogo ANOVA I Factor:
¿Qué pasa cuando todo falla y no se cumplen los supuestos?
No, no lo es! Hay algunos test ANOVA y Post Hoc que permiten resolver situaciones si hay supuestos que no se cumplen. No voy a nombrar todos, sólo dos.
ANOVA para varianzas desiguales
Se usa el test de Welch
# ANOVA PARA VARIANZAS DESIGUALES
welch_anova_test(Insomnio_2, PrimerNoche~Fedad)## # A tibble: 1 × 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche 99 8.76 2 62.9 0.000439 Welch ANOVAPost Hoc con varianzas desiguales
Se usa el test de Games-Howell (para situaciones de heterocedasticidad)
# POST HOC CON VARIANZAS DESIGUALES
games_howell_test(Insomnio_2, PrimerNoche~Fedad, conf.level = 0.95, detailed = FALSE)## # A tibble: 3 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche 1 2 2.04 0.824 3.26 0.000472 ***
## 2 PrimerNoche 1 3 1.18 0.113 2.24 0.027 *
## 3 PrimerNoche 2 3 -0.863 -2.14 0.410 0.242 nsPost Hoc contra grupo control
También a veces se da el caso de que sólo queremos comparar los diferentes grupos contra un solo grupo control (es decir, no comparar todas las medias entre sí sino todas contra nuestra Ȳ0 digamos). Para eso, usamos el test de Dunnet:
DunnettTest(x=Insomnio_2$PrimerNoche, g=Insomnio_2$Fedad) #De la librería DescTools##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`1`
## diff lwr.ci upr.ci pval
## 2-1 2.041687 0.93119563 3.152178 0.00015 ***
## 3-1 1.178872 0.06838044 2.289363 0.03567 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Acá vemos como nos toma como grupo control el primer grupo etario (menores de 20) y compara las horas dormidas de los otros dos. Esto no tiene sentido en los hechos, pero bue, sirve para mostrar como sería.
Más ejercicios
Voy a retomar los distintas ejercitaciones, empezando por la segunda parte del visto en clase:
Continuación Ejemplo 1
El gabinete psicológico desea determinar si existen diferencias en la cantidad media de horas dormidas por los pacientes con insomnio durante la primera noche en los tres grupos de edad.
e. Repetir los ítems anteriores pero considerando la tercera noche
boxplot(Insomnio_2$TerceraNoche ~ Insomnio_2$Edad)
Parecería que las horas dormidas por grupo etario son idénticas.
Test de Shapiro para evaluar si la distribución es normal.
#PdH para evaluar la normalidad: Shapiro-Wills
Insomnio_2 %>%
group_by(Edad) %>%
shapiro_test(TerceraNoche) #la función shapiro_test me permite testear si los grupos tienen una distribución normal## # A tibble: 3 × 4
## Edad variable statistic p
## <dbl+lbl> <chr> <dbl> <dbl>
## 1 1 [Menor a 20 años] TerceraNoche 0.985 0.909
## 2 2 [Entre 20 y 25 años] TerceraNoche 0.985 0.909
## 3 3 [Mas de 25 años] TerceraNoche 0.985 0.909No rechazamos la Hipótesis Nula, las distribuciones son normales para nuestro nivel de significancia. Ahora realicemos un test de igualdad de las varianzas (homocedasticidad).
Insomnio_2 %>% levene_test(TerceraNoche ~ as.factor(Edad))## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 96 4.01e-31 1El resultado es 1. Es decir, tienen la misma varianza.
Ahora, ¿existirán diferencias para un ⍺=0,05?
Insomnio_2$Fedad <- as.factor(Insomnio_2$Edad)
test_anova = aov(TerceraNoche ~ Fedad, data = Insomnio_2)
summary(test_anova) #nos tira las stats del test## Df Sum Sq Mean Sq F value Pr(>F)
## Fedad 2 0.0 0.000 0 1
## Residuals 96 552.7 5.758Evidentemente, las medias son iguales, por lo que no rechazamos H0 . No existe diferencia entre las medias, por lo que no es necesario realizar los test PostHoc.
Ejercicio 2
2) Un psicólogo clínico está interesado en evaluar la efectividad de las siguientes tres técnicas para tratar la depresión leve: reestructuración cognitiva, capacitación asertiva y un programa de ejercicio y nutrición.
Cuarenta estudiantes que sufren de depresión leve son seleccionados al azar de la lista de espera del centro de asesoría de una Universidad y se asignaron aleatoriamente a cada una de las tres técnicas mencionadas anteriormente. Los diez restantes a un grupo control placebo. Se lleva a cabo el tratamiento durante 10 semanas, tras lo cual la depresión se mide mediante el inventario de depresión de Beck (que es un cuestionario auto administrado que consta de 21 preguntas de respuesta múltiple con la finalidad de medir la severidad de una depresión).
A continuación se muestran las puntuaciones de depresión obtenidas después del tratamiento. Las puntuaciones más altas indican una mayor depresión.
Armo un Data Frame en base a la información dada.
#Primero construyo la tabla vectorizando cada atributo
placebo = c(27,16,18,26,18,28,25,20,24,26)
reest_cogn = c(10,8,14,16,18,8,12,14,9,7)
capac_cogn = c(16,18,12,15,9,13,17,20,21,19)
ejer_y_nutri = c(26,24,17,23,25,22,16,15,18,23)
#Podría armar un DF
depresion = data.frame(placebo,
reest_cogn,
capac_cogn,
ejer_y_nutri)
depresion## placebo reest_cogn capac_cogn ejer_y_nutri
## 1 27 10 16 26
## 2 16 8 18 24
## 3 18 14 12 17
## 4 26 16 15 23
## 5 18 18 9 25
## 6 28 8 13 22
## 7 25 12 17 16
## 8 20 14 20 15
## 9 24 9 21 18
## 10 26 7 19 23a) Boxplot
a) Realice un gráfico boxplot que permita comparar las puntuaciones medias para los tratamientos considerados ¿Qué puede concluir?
boxplot(depresion)
A simple vista pareciera haber dos grupos de observaciones: mientras que el grupo placebo y el del tratamiento de ejercicio y nutrición evidencian una distribución asimétrica a la derecha, aquellos que recibieron un tratamiento de reestructuración cognitiva y capacidad asertiva parecieran estar más normalmente distribuidos en sus puntajes de depresión.
Por otro lado, da la impresión de que las medias muestrales entre tres grupos son distintas entre sí (sobre todo el grupo placebo y el de ejercicio y nutrición contra los otros dos).
b) Verificación de supuestos
b) Realice la verificación de los supuestos para realizar un ANOVA- (Normalidad y Homogeneidad de varianzas)
tratamientos = c(colnames(depresion)) #armo un vector con los diferentes tratamientos y se lo paso al shapiro test
depresion %>%
shapiro_test(tratamientos)## # A tibble: 4 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 capac_cogn 0.967 0.860
## 2 ejer_y_nutri 0.901 0.225
## 3 placebo 0.887 0.157
## 4 reest_cogn 0.929 0.440El P-Value para los cuatro casos es más alto que nuestro nivel de significancia, por lo que no rechazamos nuestra hipótesis nula: las distribuciones pueden tomarse como normales (nótese de todas maneras cómo las dos variables que veíamos en el boxplot como más asimétricas presentan un p-value más bajo)
Ahora, el test de Levene para chequear homocedasticidad
id_semana = c(1:10)
depresion$id_semana = id_semana #Construyo un id de semana para el DF original
#Armo aparte un DF "largo" con la funcion melt, que me va a permitir después trabajar sobre los valores factorizados
long_depresion = melt(depresion, id.vars = 'id_semana') #le paso como id la semana, porque sino se vuelve medio loco ya que toma como id el factor a comparar.
long_depresion## id_semana variable value
## 1 1 placebo 27
## 2 2 placebo 16
## 3 3 placebo 18
## 4 4 placebo 26
## 5 5 placebo 18
## 6 6 placebo 28
## 7 7 placebo 25
## 8 8 placebo 20
## 9 9 placebo 24
## 10 10 placebo 26
## 11 1 reest_cogn 10
## 12 2 reest_cogn 8
## 13 3 reest_cogn 14
## 14 4 reest_cogn 16
## 15 5 reest_cogn 18
## 16 6 reest_cogn 8
## 17 7 reest_cogn 12
## 18 8 reest_cogn 14
## 19 9 reest_cogn 9
## 20 10 reest_cogn 7
## 21 1 capac_cogn 16
## 22 2 capac_cogn 18
## 23 3 capac_cogn 12
## 24 4 capac_cogn 15
## 25 5 capac_cogn 9
## 26 6 capac_cogn 13
## 27 7 capac_cogn 17
## 28 8 capac_cogn 20
## 29 9 capac_cogn 21
## 30 10 capac_cogn 19
## 31 1 ejer_y_nutri 26
## 32 2 ejer_y_nutri 24
## 33 3 ejer_y_nutri 17
## 34 4 ejer_y_nutri 23
## 35 5 ejer_y_nutri 25
## 36 6 ejer_y_nutri 22
## 37 7 ejer_y_nutri 16
## 38 8 ejer_y_nutri 15
## 39 9 ejer_y_nutri 18
## 40 10 ejer_y_nutri 23long_depresion %>% levene_test(value ~ variable) #test de Levene## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 36 0.110 0.954El Test de Levene demuestra que existe homocedasticidad entre las varianzas (vemos que nuestro P-Value es re mil elevado). Ahora que están cumplidos los requisitos, vamos a hacer un test de ANOVA para ver si existe una diferencia entre las medias de cada grupo.
c) Test de ANOVA
c) ¿Puede afirmar que existen diferencias? Utilizar ⍺ = 0,05
one_way_aov = aov(value ~ variable, data = long_depresion) #en mi long_depresion DF tengo sólo dos columnas: una numerica (value) y otra factorizada (variable, que corresponde a los tipos de tratamiento). Directamente pasandole esos parámetros a la fórmula y señalando donde está la Data ya resolvimos todo con una línea de código.
summary(one_way_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## variable 3 762.9 254.29 15.92 9.4e-07 ***
## Residuals 36 574.9 15.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como podemos ver el valor del estadístico F de la función de probabilidad F de Snedecor nos arroja un valor altísimo (F = 15,92) en la comparación de la varianza de los cuatro grupos de tratamiento para la depresión (pensemos que si las varianzas fueran idénticas F sería 1). Esto se ve reflejado en el P-Value que es hipermegaarchiquito (p_value = 9,4e-07), por lo que para un nivel de significancia ⍺ = 0,05, rechazamos nuestra H0, las medias son diferentes. Pero… ¿entre cuáles medias se dan esas diferencias?
d) Test Post Hoc
d) De existir diferencias realizar los test post hoc y extraer conclusiones.
Podríamos hacer cualquiera de los Test Post Hoc planteados antes (Tukey, Bonferroni, bla), pero lo cierto es que en este caso particular estamos trabajando contra un grupo control (Placebo), por lo que voy a optar por el test de Dunnet, que justamente evalúa eso.
test_dunn = DunnettTest(value ~ variable, data = long_depresion, control = 'placebo') #con el parámetro "control=" le paso al test cuál es el factor contra el que quiero comparar el resto
test_dunn##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $placebo
## diff lwr.ci upr.ci pval
## reest_cogn-placebo -11.2 -15.582618 -6.817382 8.3e-07 ***
## capac_cogn-placebo -6.8 -11.182618 -2.417382 0.0015 **
## ejer_y_nutri-placebo -1.9 -6.282618 2.482618 0.5859
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Esto ya se veía venir (?) un poco en el boxplotm pero ahora tenemos las métricas exactas. Evidentemente existe una diferencia estadísticamente significativa entre las medias de los puntajes de depresión (calculados a partir del inventario de depresión de Beck) entre aquellos pertenecientes al grupo placebo y los que recibieron un tratamiento de Reestructuración Cognitiva (el Test de Dunnet arroja un p_value = 8,3e-07) y entre el grupo placebo y y los que recibieron un tratamiento de Capacidad Asertiva (p_value = 0,0015).
Sin embargo, no se nota una diferencia estadísticamente significativa entre el grupo control y aquellos a los que se les asigno como tratamiento Ejercicio y Nutrición. Es decir, hagan terapia, no sean Ivana Nadal.
No, ta bien, hagan ejercicio y coman variado, pero posta ¿no estaría el mundo mucho mejor si algunas personas aceptaran que necesitan un poco de ayuda profesional? En fin, la hipotenusa. Y hablando de la hipotenusa, por último vamos a hacer unos gráficos para ver como es la dispersión y comparación de estas medias.
plot(test_dunn)
Pues bien, vemos que en este sencillo gráfico, de abajo para arriba, las diferencia de medias contra el grupo control (Placebo) entre los tratamientos de Ejercicio y Nutrición, el de Capacidad Asertiva y el de Reestructuración Cognitiva.
Puede visualizarse lo que decíamos antes, el tratamiento de Ejercicio y Nutrición no es estadísticamente significativo respecto al grupo control, y se expresa en que su diferencia de media está de hecho sobre el 0. Los otros dos, más a la izquierda, estarían mostrando una disminución del puntaje de depresión, respecto al grupo control.
e. ANOVA de dos o más factores
Hasta acá, lo visto aplica para una variable cuantitativa cuyo comportamiento queríamos comparar en base a situaciones agrupadas en distintas variables cualitativas.
En el ejemplo del Ejercicio 1, lo que veíamos era si entre tres grupos distintos del mismo factor (grupo etario) podía notarse diferencias estadísticamente significativas en la cantidad de horas dormidas (variable dependiente) en la Primer Noche en base a un tratamiento para el insomnio.
Ahora bien, ¿qué pasaría si quisiéramos, además de ver del efecto de la edad, existe efecto según el género en la cantidad de horas dormidas entre los participantes del estudio en la Primer Noche?
Esto nos lleva a un análisis ANOVA de dos factores, donde las preguntas a hacernos serían:
- ¿Hay efecto debido a la edad?
- ¿Hay efecto debido al género?
- Y, lo más importante e “innovador”: ¿Hay un potenciamiento o una atenuación del efecto de estos factores entre sí?
De lo que se trata entonces es de estudiar las combinaciones que se dan entre el factor 1 y el factor 2.
Si antes teníamos entonces una tablita así:
Donde cada 𝑦𝑖𝑗 podía expresarse como 𝑦𝑖𝑗 = 𝜇 + 𝛼𝑖 + 𝜀𝑖𝑗 (es decir, nuestra variable dependiente o, en este caso, horas dormidas es igual a la media general de la población insomne 𝜇 ,más el efecto del grupo etario 𝛼𝑖, más “algo” (??) propio del individuo 𝜀𝑖𝑗).
Ahora, ¿qué pasa cuando hay dos factores?
Volcamos los datos de manera similar antes, pero agregando la variable cualitativa “género” (en este caso varón y mujer).
Si queremos podemos, como en el caso de ANOVA de un factor, calcular las medias de cada uno de estos géneros para cada grupo de edad, pero lo distinto es que podemos conformar nuevos grupos que surjan de la combinación de los valores de las variables cualitativas o factores dados.
Entonces, si antes los grupos a comparar venían dados por las modalidades de ese único factor (antes, en el ejemplo, eran sólo los tres grupos etarios lo que tenía que comparar), con dos factores la definición de los grupos a comparar vienen dados de la combinación de las modalidades de las dos variables (en este caso, como tenemos dos modalidades para el factor género tendríamos seis grupos a comparar). Esta es la primera diferencia con respecto a ANOVA de un factor.
Existe una interacción entre dos factores si el efecto de uno de los factores cambia en las diferentes categorías del otro factor.5
La segunda diferencia importante respecto al ANOVA de un factor es que yo en ese entonces tenía una sola PdH contrastar (H0 = no hay diferencia entre los grupos vs. Ha = hay diferencia entre los grupos); para ANOVA de dos factores hay, en cambio, tres PdH:
PdH1: No hay diferencia debidas al efecto del factor 1 (digamo’, grupos de edad).
PdH2: No hay diferencias debidas al efecto del factor 2 (en nuestro ejemplo, género).
PdH3: No hay interacción (o diferencias debidas a la interacción) entre los factores.
Matemáticamente, si para ANOVA de un factor teníamos que 𝑦𝑖𝑗 = 𝜇 + 𝛼𝑖 + 𝜀𝑖𝑗, para ANOVA de dos factores:
𝑦𝑖𝑗 = 𝜇 + 𝛼𝑖 + 𝛽𝑘 + 𝛾𝑖𝑘 + 𝜀𝑖𝑗𝑘
Donde, siguiendo nuestro ejemplo: nuestra variable dependiente (en este caso, horas dormidas) es igual a la media general de la población insomne 𝜇 ,más el efecto del factor 1 𝛼𝑖 (grupo etario), más el efecto del género (𝛽𝑘), más la interacción entre los grupos (simbolizada con la letra gamma, 𝛾𝑖𝑘), más “algo” (??) propio del individuo 𝜀𝑖𝑗𝑘 .
Todo este choclo en realidad sintetiza la sumatoria de cuadrados que vendría a ser todo este choclo:
Que después, para realizar nuestras PdH, nos llevará a calcular distintos estadísticos F, donde tendremos para cada una de estas PdH una H0 y Ha -como se dijo antes- que evaluaremos mediante los P-Value obtenidos según los grados de libertad dados.
Supuestos para ANOVA de dos factores
El ANOVA de dos factores es un test de hipótesis paramétrico que busca evaluar hipótesis sobre el valor de medias poblacionales (parámetro) de varios grupos (muestras independientes), controlando por otra variable.
Los supuestos son básicamente los mismos que para el ANOVA de un factor:
- Para cada celda, los valores muestrales provienen de una población con distribución aproximadamente normal (este procedimiento es robusto frente a desviaciones razonables de las distribuciones normales).
- Las poblaciones tienen la misma varianza 𝞂2 (o desviación estándar 𝞂). Este procedimiento es robusto frente a desviaciones razonables del requisito de igualdad de varianzas).
- Las muestras son aleatorias simples de datos cuantitativos.
- Las muestras son independientes entre sí (este procedimiento no se aplica a muestras que no son independientes).
- Los valores muestrales se categorizan en dos factores (esta es la base del nombre del método: análisis de varianza de dos factores).
- Todas las celdas tienen el mismo número de valores muestrales (este diseño se conoce como diseño balanceado).
Sintaxis en R
Importante aclaración sobre la aplicación de ANOVA de dos factores en R, antes de pasar al ejercicio de ejemplo:
Para realizar un ANOVA de dos factores en R utilizaremos la función aov(), como se hizo previamente pero con dos salvedades: si queremos analizar los efectos de los factores por separado debemos utilizar la función aov(y ~ 𝛼 + 𝛽, data = dataframe). Por otro lado, si queremos analizar la interacción debemos usar aov(y ~ 𝛼 * 𝛽 , data = dataframe)
Aunque no tomemos en cuenta la interacción, medir el efecto de los factores por separado mediante el ANOVA de dos factores no es lo mismo que realizar dos análisis de ANOVA en paralelo, ya que la función aov(y ~ 𝛼 + 𝛽, data = dataframe) contempla también los errores
Ejemplo 1
Vamos ahora a aplicar todos nuestros relucientes conocimientos a nuestra vieja tablita de Insomnio. Voy a volver a levantarla por si se olvidaron cómo venía la mano:
urlfile_1 = 'https://github.com/oblitterator/estadistica_II_UNTREF/blob/main/Insomnio_2.sav?raw=true'
Insomnio_2 = read_sav(url(urlfile_1))
print(names(Insomnio_2))## [1] "Clave" "PrimerNoche" "SegundaNoche" "TerceraNoche" "Sexo"
## [6] "Edad" "Terapia1" "Terapia2" "Terapia3"print(nrow(Insomnio_2))## [1] 99Factorizamos las dos columnas a estudiar y, primero chequeamos el cumplimiento de los supuestos, previo a la realización del ANOVA:
Insomnio_2 %>%
group_by(Sexo,Edad) %>%
get_summary_stats(PrimerNoche, type = "mean_sd")## # A tibble: 6 × 6
## Sexo Edad variable n mean sd
## <dbl+lbl> <dbl+lbl> <chr> <dbl> <dbl> <dbl>
## 1 0 [Mujer] 1 [Menor a 20 años] PrimerNoche 17 5.60 1.76
## 2 0 [Mujer] 2 [Entre 20 y 25 años] PrimerNoche 16 8.67 1.89
## 3 0 [Mujer] 3 [Mas de 25 años] PrimerNoche 17 6.53 2.23
## 4 1 [Hombre] 1 [Menor a 20 años] PrimerNoche 16 5.80 1.65
## 5 1 [Hombre] 2 [Entre 20 y 25 años] PrimerNoche 17 6.86 2.49
## 6 1 [Hombre] 3 [Mas de 25 años] PrimerNoche 16 7.25 1.49En primer lugar, testeamos la normalidad, previo a factorizar:
# testeo de normalidad de las variables para cada grupo!
Insomnio_2 %>%
group_by(Sexo,Edad) %>%
shapiro_test(PrimerNoche)## # A tibble: 6 × 5
## Sexo Edad variable statistic p
## <dbl+lbl> <dbl+lbl> <chr> <dbl> <dbl>
## 1 0 [Mujer] 1 [Menor a 20 años] PrimerNoche 0.953 0.500
## 2 0 [Mujer] 2 [Entre 20 y 25 años] PrimerNoche 0.954 0.550
## 3 0 [Mujer] 3 [Mas de 25 años] PrimerNoche 0.958 0.600
## 4 1 [Hombre] 1 [Menor a 20 años] PrimerNoche 0.965 0.746
## 5 1 [Hombre] 2 [Entre 20 y 25 años] PrimerNoche 0.973 0.868
## 6 1 [Hombre] 3 [Mas de 25 años] PrimerNoche 0.970 0.842ggqqplot(Insomnio_2, "PrimerNoche", facet.by = "Edad*Sexo") #este último argumento es el que le pasamos a la función para que "entienda" que tiene que graficar todos los grupos
Todos los p-value nos dan arriba de 0,95, por lo que podemos estar tranquilos de que nuestra H0: todos los grupos tienen distribución normal,
Y ahora, vemos el comportamiento de las varianzas (es decir, que sean homogéneas) mediante el Test de Levene:
# testeo de homocedasticidad de varianzas
Insomnio_2$Edad = as.factor(Insomnio_2$Edad)
Insomnio_2$Sexo = as.factor(Insomnio_2$Sexo)
Insomnio_2 %>% levene_test(PrimerNoche ~ Edad*Sexo) #este último argumento ("multiplicando" los dos factores) es el que le pasamos a la función para que "entienda" que tiene que comparar## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 5 93 1.28 0.278Nuestro p-value de 0,27 nos indica que existe homocedasticidad. Ahora ya estamos preparades para realizar nuestro ANOVA de dos factores, que sería una cosa así:
################################################################
test_anova_2 = aov(PrimerNoche ~ Edad+Sexo, data = Insomnio_2)
summary(test_anova_2) #nos tira las stats del test## Df Sum Sq Mean Sq F value Pr(>F)
## Edad 2 69.3 34.66 8.546 0.000386 ***
## Sexo 1 2.2 2.17 0.536 0.465948
## Residuals 95 385.3 4.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como observamos por el p-value obtenido, el efecto de los grupos etarios es estadísticamente significativo (F= 8,546 y p-value = 0,000386)6 para explicar las diferencias entre las medias de horas dormidas entre los participantes del estudio con un nivel de significación usual (⍺ = 0,05), mientras que el factor “Sexo” no lo sería (su p-value = 0,46).
Estas, si seguimos esquemáticamente nuestra presentación del tema, serían nuestras PdH1 y PdH2. Ahora bien, ¿existe interacción entre ambos factores?
test_anova_2_cross = aov(PrimerNoche ~ Edad*Sexo, data = Insomnio_2)
summary(test_anova_2_cross) #nos tira las stats del test## Df Sum Sq Mean Sq F value Pr(>F)
## Edad 2 69.3 34.66 9.058 0.000255 ***
## Sexo 1 2.2 2.17 0.568 0.452977
## Edad:Sexo 2 29.4 14.71 3.843 0.024907 *
## Residuals 93 355.9 3.83
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Por lo que se aprecia, la combinación de los factores de grupo etario y sexo influye en las horas de sueño (p-value de Edad:sexo = 0,024907). Es decir, la presentación de un efecto está condicionada por la presencia de otro factor.
¿Qué significa la interacción entre factores?
Ahora que sabemos que existe una interacción entre los factores (que la aparición de uno de ellos hace que el otro sea diferente) y que esta produce un efecto la pregunta que debemos hacernos es: ¿se trata de un potenciamiento o de una atenuación del efecto?
En primer lugar, vamos a intentar entender esto con una representación gráfica, partiendo de nuestra H0 para la PdH3: No existe interacción entre los factores. Si esto fuese así, al graficar nuestros factores, nos daría por ejemplo:
Una ausencia de interacción en la combinación de los factores se visualiza como dos rectas paralelas (ascendentes, descendentes, curvilíneas, etc.), como se ve más arriba.
Ahora, cuando hay interacción la cosa se ve más o menos así:
Continuación Ejemplo 1
Volviendo a nuestro ejercicio, esto es lo que pasa (?) cuando graficamos nuestros factores “Edad” y “Sexo”.
# gráfico de interacción
with(Insomnio_2, interaction.plot(Edad, Sexo, PrimerNoche, fun = mean,
main = "Interaction Plot"))
La línea punteada corresponde a las mujeres y la línea continua corresponde a los varones, mientras que en el eje horizontal encontramos los tres grupos de edad de manera consecutiva.
¿Qué pasa con el promedio de horas dormidas en la primer noche? Podemos ver que, mientras que en el caso de los varones existe un aumento progresivo en las horas dormidas por grupo de edad (crece a medida que la edad de los sujetos que es mayor). Por otro lado, en las mujeres el comportamiento es distinto: si bien termina siendo creciente el promedio de horas dormidas, en el segundo grupo (de entre 20 y 25 años) que presenta un promedio muchísimo más alto con respecto a los otros grupos (un efecto pico, que en el grupo de los varones está muy atenuado).
Este gráfico muestra de manera muy clara la interacción entre dos factores. O en otras palabras, que la combinación de los factores tiene un efecto estadísticamente significativo.
Lo que nos queda ver es, entonces, cómo medimos ese efecto (el tamaño de ese efecto) que hemos encontrado como estadísticamente significativo.
f. Tamaño del efecto (la fuerza de la relación)
El tamaño del efecto puede pensarse como el grado de influencia de un factor sobre la variable dependiente en un ANOVA. Es decir, intenta ser explicativo sobre una correlación: nos indica la fuerza de la misma (no sólo si existe o su forma, brindada por los Test Post Hoc).
El tamaño del efecto se analiza mediante el denominado eta cuadrado o con el eta cuadrado parcial, depende del tipo de ANOVA que estemos haciendo.
Si estamos realizando un ANOVA de un factor utilizamos el eta cuadrado que es en sí la suma de cuadrados del efecto (SSEntre) dividido la suma de cuadrados total (SSTotal). En otras palabras sería “cuánto explica la variabilidad entre los grupos la variabilidad total”.
Por otro lado, el eta cuadrado parcial se aplica a ANOVA de más de un factor, que trae aparejado una descomposición de la variabilidad en más partes. Se llama “parcial” porque descarta la suma de cuadrados de aquellos factores o efectos que no serán considerados (sólo calculo el tamaño de los efectos que son previamente identificados como estadísticamente significativos)
¿Qué valor debe tomar eta cuadrado o eta cuadrado parcial para considerar que el efecto es “grande” (o pequeño)? Obviamente, no existe acuerdo a nivel global pero si existen unos valores de corte, aceptados en la comunidad estadística:
¿Cómo se calcula el tamaño del efecto con R?
II. Regresión Lineal
a. Regresión Simple
Ejemplo 2
Supongamos que hay un problema como el siguiente:
Un docente desea estudiar la relación entre el rendimiento escolar y la ansiedad, partiendo de algunas teorías que expresan que a determinado nivel de ansiedad hay un nivel de rendimiento óptimo pero fuera de esos márgenes el rendimiento se resiente. Para eso, toma una muestra de 10 individuos y les aplica un instrumento psicométrico de medición de ansiedad y luego les toma el examen, dando los siguientes resultados.
#creo los vectores y armo el data frame
ansiedad = c(28, 41, 35, 39, 31, 42, 50, 46, 45, 37)
puntuacion_examen = c(82, 58, 63, 89, 92, 64, 55, 70, 51, 72)
ansiedad_examen = data.frame(ansiedad,puntuacion_examen)
ansiedad_examen## ansiedad puntuacion_examen
## 1 28 82
## 2 41 58
## 3 35 63
## 4 39 89
## 5 31 92
## 6 42 64
## 7 50 55
## 8 46 70
## 9 45 51
## 10 37 72¿Cómo puede saber si existe vinculación entre la puntuación de ansiedad y la puntuación del exámen?
Este problema es un problema clásico de asociación entre dos variables cuantitativas, cuya aproximación puede darse de dos maneras:
Utilizando gráficos de dispersión
Calculando indicadores numéricos (Covarianza/Coeficiente de Regresión Lineal)
El análisis de asociaciones de variables cualitativas con cuantitativas se tratará en la próxima sección de Regresión Logística.
Gráfico de Dispersión
La resolución gráfica de este problema deriva en un gráfico de dispersión, donde cada individuo se grafica utilizando los puntajes que obtuvo en las dos variables, como si fueren las coordenadas en dos ejes cartesianos.
plot(ansiedad_examen) #con plot y el DF armo un sencillo gráfico de dispersión
En el eje Y podemos ver el rendimiento (puntaje en el examen) y en el X la puntación de ansiedad. Parecería, grosso modo, que exuste una tendencia decreciente en la ubicación relativa de los individuos: a mayor nivel de ansiedad menor rendimiento.
Este tipo de gráficos sirven para visualizar tendencias, posibles relaciones y forma de relación entre las variables.
Ahora, ¿cómo lo cuantificamos?
Covarianza y Coeficiente de Correlación de Pearson
No necesariamente la relación entre dos variables se restrinje al tipo lineal, aunque siempre en Estadística este es el punto de partida para abordar luego otro tipo de relaciones, por dos grandes motivos:
Es el más sencillo de evaluar
Muchos problemas de relaciones curvilíneas pueden aproximar su resolución a través del caso lineal
Covarianza
La covarianza, como indicador, permite evaluar la asociación de tipo lineal entre dos variables.
Debe tomarse en cuenta los siguientes recaudos:
La misma permite evaluar asociaciones del tipo lineal
Puede tomar cualquier valor
Depende (o mejor dicho, está atada) a las unidades de medida de las variables
Este valor numérico que arroja la covarianza (cov) se interpreta de la siguiente forma:
cov(X,Y) > 0 : Nos indica una asociación lineal del tipo directa (ambas crecen/decrecen juntas).
cov(X,Y) < 0 : Nos indica una asociación lineal del tipo inversa (mientras una crece la otra decrece)
cov(X,Y) = 0 : No existe asociación lineal entre las variables
Los problemas de la covarianza con indicador son los expuestos previamente: al poder tomar cualquier valor, no nos permite medir la fuerza entre dos variables (¿2, 5, un millón, veinte moles?)
El segundo problema es que, como su valor depende de la unidad de medida, si utilizamos -por ejemplo- kilogramos y metros para medir individuos y en otro caso utilizamos libras y pies, nos dan valores distintos, lo que no es adecuado para un indicado que mide asociación ¿Por qué? No hay por qué. No sí, hay por qué: porque si lo que queremos medir es asociación, esa se supone intrínseca a las variables en cuestión, y no debería generarse ruido por los valores de las unidades en cuestión.
¿Y ahora qué hacemos?
Coeficiente de Correlación de Pearson (r)
Lo que hacemos es “corregir” la covarianza (y sus buenas virtudes) con el Coeficiente de Correlación de Pearson (r).
Esto es una covarianza estandarizada (dividida por los desvíos estándar de las variables):
Este coeficiente permite evaluar asociaciones del tipo lineal entre dos variables cuantitativas, pero con la siguiente ventaja:
- Está acotado entre -1 y 1 (-1<r<1)
No se encuentra, por tanto, interferido por las unidades de medida de las variables
Continuación Ejemplo 2
En el caso de ejemplo que venimos viendo, podemos calcular el r fácilmente
cor(ansiedad_examen) #esto esta tirado de manera muy cabeza, y en caso de tener una tabla con más variables habría que hacer otros ajustes de sintaxis, pero a los hechos sirve## ansiedad puntuacion_examen
## ansiedad 1.0000000 -0.6907746
## puntuacion_examen -0.6907746 1.0000000Como se puede observar, el valor de la correlación arroja un r = -0,69, lo cuál iría de la mano de nuestra primera aproximación: existe una relación inversa lineal entre la puntuación del examen y la obtenida en el test de ansiedad.
No sólo podemos verificar ahora la forma de esta relación sino que, al mismo tiempo, podemos ver la fuerza de la misma. Por ejemplo, y si bien no existe un consenso total, en ciertos casos, un r >0,5 (ó -0,5), indica una relación “alta” entre las variables7.
Veamos otras formas podría mostrarnos diferentes valores de r.
Cuando r = 1 ó r = -1, significaría que todos los valores caen sobre la recta (están más “alineados”). A su vez, a medida que este valor se van alejando, existe mayor dispersión respecto a la recta.
Cuando r = 0, no existe relación (mucha dispersión) o todas caen sobre una recta paralela al eje de las abscisas.
Tip: a la hora de avanzar con la práctica estadística es importante hacer una buena EDA mediante gráficos y, aunque no se vea una relación lineal a simple vista, calcular el Coeficiente de Correlación.
Calcular el r es un abordaje preliminar de todas formas ya que muchas veces el interés está guiado a ver si podemos -en base a la información dada- realizar una predicción de nuestra variable dependiente.
Modelo Regresión Lineal
Digamos que, tomando como base nuestro ejemplo, el profesor de la clase quiere saber si puede predecir el puntaje que va a obtener un alumno en el exámen en base al valor de ansiedad obtenido en base al test psicométrico. O quizá quiere determinar si tomarle o no el exámen según el puntaje que sacó en el test de ansiedad. Esto es el basamento de los análisis de los modelos de regresión lineal.
Vincular las variables en función de poder predecir o explicar el comportamiento de una variable que llamamos dependiente en función de una o más variable predictoras o independientes.
A diferencia de lo visto en caso de la correlación -donde es indistinto el rol de una variable o la de otra-, en este caso estamos pensando en variables que en el estudio tienen distintos roles.
La ecuación (o recta de regresión lineal simple) para una variable es:
Este modelo vincula dos variables: una dependiente (yi) y una independiente o predictora (xi). Este vínculo está signado por dos parámetros que determinan el modelo de relación entre las variables:
⍺ = ordenada al origen o “intercepto”.
𝛃 = pendiente de la recta o “efecto” de xi sobre yi.
𝝴i = es el término de error de nuestro modelo, desconocido y que suponemos aleatorio ya que nuestra muestra se supone extraída de la misma forma, que surge de considerar que nuestro modelo no es perfecto. La media de todos los 𝝴i se supone también 08, por lo que posteriormente lo eliminamos de nuestra ecuación. Estos supuestos (normalidad), deben ser chequeados previamente a la conformación de nuestro modelo.
Tanto ⍺ como 𝛃 son parámetros poblacionales, y por tanto desconocidos, por lo que debemos estimarlos a través de los datos disponibles.
Llamamos recta de regresión estimada o ajustada a la ecuación que presentamos con los estimadores correspondientes para cada uno de los parámetros.
ŷi = a + bxi +ei
Los residuos (ei) y el método de los mínimos cuadrados para calcular a y b
Si yi es el valor verdadero de la variable en nuestra población ŷi es el resultado predicho para un determinado valor de xi , dados los parámetros estimados a y b, como se vio en la recta anterior.
¿Qué es entonces ei? El residuo o diferencia entre el valor real de la variable dependiente (yi) respecto al valor predicho (ŷi). Es decir que tendremos un ei = yi - ŷi , para cada una de las variables predichas.
¿Para qué nos sirve? ¡Para qué no nos sirve! El residuo ei describe el error del ajuste del modelo en cada punto (es decir, cuán bien ajusta el mismo). En un caso ideal, todos los ei debería tender a 0.
Pero no es sólo eso: los estimadores a y b se calculan, justamente, a través de los residuos, a través del Método de los Mínimos Cuadrados.
El método de los Mínimos Cuadrados para el cálculo de la recta
Para encontrar los valores a y b, intentamos que coincidan con el valor mínimo de la suma de los cuadrados de los residuos (en otras palabras, buscamos los estimadores que minimizan la suma de los cuadrados de los residuos). Esto se denomina “estimadores de mínimos cuadrados” (MMC).9
Gráficamente podemos verlo de la siguiente manera:
En este caso, esta línea de ajuste óptimo (con ordenada al origen a y pendiente b) sería la que nos estaría brindando el menor valor entre la sumatoria de los cuadrados de los residuos (la menor “distancia”). No debería existir otra recta que, en promedio, nos permitiera conseguir residuos más pequeños.
Testeo del nivel de significación de los estimadores a y b
Ahora bien, una vez calculados nuestros estimadores de la recta de regresión, podemos (o debemos) realizar una PdH para testear si respecto al valor de los parámetros (⍺ ó 𝛃) son estadísticamente significativos. Es decir, si :
H0: ⍺ ó 𝛃 = 0
H1: ⍺ ó 𝛃 ≠ 0
Para tal fin, se trabaja con el estadístico de prueba t en el caso del estimador de 𝛃 (con los supuesto correspondientes),siendo:
Mientras que para el caso del estimador de ⍺:
Esto nos arrojará un determinado p-value para los el estadístico t. En caso de que el p-value < 0,05, rechazaremos nuestra H0, pudiendo afirmar que los estimadores son estadísticamente significativos.10
En general, si el r previamente arrojó un valor alto, los estimadores de la recta van a dar estadísticamente significativos. Sí puede darse, en caso de muestras pequeñas, que el r sea alto pero esto no se cumpla (las PdH son muy sensibles al tamaño de la muestra!).
Testo del nivel de significación de la varianza
De la misma manera, hay un test ANOVA asociado a la regresión lineal que permite ver si la regresión es significativa o no.
∑(yi - ȳ)2 ó SST: es la suma de cuadrados de yi y refleja su variabilidad en torno a su media (es decir, la variabilidad de yi).
∑(ŷi - ȳ)2 ó SSR: es la suma de cuadrados de la regresión y refleja la cantidad de variación de los valores de yi que pueden ser explicados por el modelo.
∑(yi - ŷi )2 o SSE: es la suma de los cuadrados de los residuos vistos anteriormente.
En otras palabras: SST = SSR + SSE
Si el modelo es bueno, el SSR (la variabilidad explicada por el modelo) debería ser muy similar al SST (la variabilidad “verdadera” de y), con un SSE bajo (es decir, un error mínimo). Para comprobar esto lo que se realiza es un test del tipo ANOVA para ver si SST es similar a SSR (lo que implicaría que SSR sea mucho más grande que SSE).
Este test ANOVA se realiza, como en los casos anteriores, a través del cálculo del estadístico F:
Donde, previamente se tiene que:
Sintaxis para el Modelo de Regresión Lineal en R
lm(formula, data, subset, weights, na.action, method = “qr”, model = TRUE, x = FALSE, y =FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, …)
Continuación Ejemplo 2
Si, como se dijo antes, nuestro profesor de estadística quiere predecir en cuánto y cómo afecta la ansiedad en la puntuación del examen de sus alumnos, podría tomar la información volcada en la tabla:
ansiedad_examen## ansiedad puntuacion_examen
## 1 28 82
## 2 41 58
## 3 35 63
## 4 39 89
## 5 31 92
## 6 42 64
## 7 50 55
## 8 46 70
## 9 45 51
## 10 37 72Y utilizar la sintaxis previamente presentada para calcular los coeficientes a y b.
#Calculo los coeficientes utilizando la función lm y los guardo en variables separadas para presentar la ecuación
modelo_ansiedad = lm(ansiedad_examen$puntuacion_examen ~ ansiedad_examen$ansiedad)
est_alpha = modelo_ansiedad$coefficients[1]
est_beta = modelo_ansiedad$coefficients[2]
print(paste('En este caso nuestra recta de regresión sería: y =', round(est_alpha,2), '+',round(est_beta,2),'* x'))## [1] "En este caso nuestra recta de regresión sería: y = 125.88 + -1.43 * x"Es decir, la nota (y) disminuye -1,43 por cada punto de ansiedad (x) obtenido en el test psicométrico.
Ahora bien, ¿cómo chequeamos si nuestros estimadores son estadísticamente significativos? Hacemos un summary() de nuestro modelo.
summary(modelo_ansiedad)##
## Call:
## lm(formula = ansiedad_examen$puntuacion_examen ~ ansiedad_examen$ansiedad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.885 -7.957 -1.457 7.507 18.829
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.8830 21.1109 5.963 0.000337 ***
## ansiedad_examen$ansiedad -1.4285 0.5287 -2.702 0.026986 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.87 on 8 degrees of freedom
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4118
## F-statistic: 7.301 on 1 and 8 DF, p-value: 0.02699Como vemos, tanto el p-value de a (intercepto) como de b, son menores a 0,05, lo que nos estaría indicando que son estadísticamente significativos respecto a nuestros parámetros. De la misma manera, el p-value del estadístico F estaría indicando que nuestro modelo de regresión es significativo.
Por último, es importante ver el valor de r2 . Este último, conocido también como Coeficiente de Determinación, representa la proporción de varianza de la variable dependiente que es explicada por el modelo. Mientras nuestro r2 sea más cercano a 1, nos estaría indicando que el mismo es “mejor”ya que es mayor la proporción de la varianza de la variable dependiente que es explicada por el mismo.
b. Regresión Múltiple
Hasta ahora hemos analizado la correlación lineal entre dos variables, pero en la mayor parte de los problemas de investigación en que se aplica un análisis de regresión, se requiere más de una variable como predictor. En la realidad uno tiene una serie de variables que pueden llegar (o no) a cumplir un rol importante en el comportamiento de una variable dependiente
Ahora, como existen muchas variables uno no puede saber a priori cuáles son más útiles. Entonces las tres preguntas principales que nos deberemos hacer al adentrarnos en la Regresión Múltiple son:
¿Cómo hacer para considerar más de un predictor?
Cómo hacer para elegir el mejor conjunto de predictores.
El hecho de considerar más de una variable independiente (o predictora), ¿trae aparejado algún problema sobre el que haya que prestar atención
Triola dixit: “debido a que los cálculos requeridos son tan complicados, su realización manual no resulta práctica y constituye una amenaza para la salud mental”. Así que ya fueron advertidos, usen algún software, R, Python, pero no haga esto a mano en su casa.
Nos centramos en los siguientes dos elementos clave:
Determinación de la ecuación de regresión múltiple y;
Uso del valor de R2 ajustado y el valor P como medidas de qué tan bien se ajusta la ecuación de regresión múltiple a los datos muestrales.
Ecuación del Modelo de Regresión Lineal Múltiple
Siempre lo que estamos haciendo es intentar buscar los estimadores para los parámetros por lo que, matemáticamente, la ecuación es similar a la de Regresión Lineal Simple, sólo que para cada xk contamos con un un parámetro 𝛃k que lo acompaña.
Entonces, dadas x1 … xk , k variables independientes, el modelo de regresión lineal múltiple para los parámetros es :
yi = ⍺ (ó 𝛃0 ) + 𝛃1 x1i + … + 𝛃k xki + 𝝴i
Donde 𝝴i es el error del modelo, que se supone aleatorio y de media cero para todo i
Y como los parámetros deben estimarse a través de los datos, nuestra recta de estimación múltiple quedaría:
yi = a (ó b0 ) + b1 x1i + … + bk xki + ei
Donde, como se vio antes, ei es el residuo o diferencia entre el valor real de la variable dependiente (yi) respecto al valor predicho (ŷi),
Y a través del Método de los Mínimos Cuadrados (MMC) de estos residuos, es que calcularemos los valores óptimos de nuestros bk para el ajuste de la recta.
Aclaración: siempre puede obtenerse una recta a través de la utilización del MMC ya que esto depende de un algoritmo geométrico, no probabilístico, no es necesario ningún supuesto sobre los errores para eso. Esto es importante porque este hecho nos permitirá en otros casos comparar los modelos de regresión lineal frente a otros posibles, aunque los supuestos no se cumplan.
Peeero es partiendo de los supuestos y su comprobación -como el supuesto de normalidad- que evaluaremos el ajuste de nuestro modelo, su variabilidad, así como la significativdad de los parámetros. Y esto lo hacemos a través del análisis de los residuos, que suponemos aleatorios.
Bondad del ajuste a través del Coeficiente de Determinación (r2) ajustado
Como se mencionó antes, el r2 representa la proporción de varianza de la variable dependiente que es explicada por el modelo.
Si nuestro modelo ajusta bien, el SSR (que es la variabilidad de nuestro modelo) debería ser similar a la SST (la variabilidad total de la muestra). En dicho caso, su similitud arrojaría una proporción cercana a 1, que daría cuenta que nuestro modelo explica bien la variabilidad real.
Norsante, en este caso -al estar trabajando con más de una variable- deberemos utilizar lo que se conoce como r2 ajustado, que lo que busca es considerar en qué medida la inclusión de nuevas variables aumenta (o no) la dispersión de nuestro modelo:
Es decir, cada vez que yo agrego una variable a un modelo, el r2 tiende a “inflarse” (de manera no significativamente, en los hechos), por lo que ajustamos el mismo para poder medir bien la relación entre las varianzas. El r2 ajustado lo que hace es “penalizar” al r2 en función de las variables agregadas.
Veamos esto con un ejemplo:
Ejemplo 3
Una medida objetiva del ajuste aeróbico de una persona es el consumo de oxígeno en volumen por peso unitario del cuerpo, por unidad de tiempo.
Se utilizaron 31 individuos en un experimento con el objeto de podermodelar el consumo de oxígeno mediante las siguientes variables.
Edad.
Peso.
Tiempo en recorrer 3 km.
Pulso en reposo.
Pulso al final del ejercicio.
Pulso máximo durante el ejercicio.
urlfile_2 = 'https://raw.githubusercontent.com/oblitterator/estadistica_II_UNTREF/main/ajuste_aerobico.csv'
aerobico = read_csv(url(urlfile_2))## Rows: 31 Columns: 7## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (7): consumo_oxigeno, edad, peso, tiempo_en_recorrer_3km, pulso_en_repos...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.aerobico## # A tibble: 31 × 7
## consumo_oxigeno edad peso tiempo_en_recorrer_3… pulso_en_reposo pulso_final
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 44.6 44 89.5 11.4 62 178
## 2 45.3 40 75.1 10.1 62 185
## 3 54.3 44 85.8 8.65 45 156
## 4 59.6 42 68.2 8.17 40 166
## 5 49.9 38 89.0 9.22 55 178
## 6 44.8 47 77.4 11.6 58 176
## 7 45.7 40 76.0 12.0 70 176
## 8 49.1 43 81.2 10.8 64 162
## 9 39.4 44 81.4 13.1 63 174
## 10 60.1 38 81.9 8.63 48 170
## # … with 21 more rows, and 1 more variable: pulso_maximo <dbl>En este caso, nuestra variable a predecir (y) es el consumo de oxígeno, y el resto son nuestras posibles variables predictoras (xk).
Vamos a analizar las correlaciones ahora entre las variables:
cor(aerobico)## consumo_oxigeno edad peso
## consumo_oxigeno 1.0000000 -0.3045924 -0.16275285
## edad -0.3045924 1.0000000 -0.23353903
## peso -0.1627528 -0.2335390 1.00000000
## tiempo_en_recorrer_3km -0.8621949 0.1887453 0.14350758
## pulso_en_reposo -0.3464054 -0.1415664 0.02270099
## pulso_final -0.3979742 -0.3378703 0.18151633
## pulso_maximo -0.2367402 -0.4329159 0.24938123
## tiempo_en_recorrer_3km pulso_en_reposo pulso_final
## consumo_oxigeno -0.8621949 -0.34640540 -0.3979742
## edad 0.1887453 -0.14156640 -0.3378703
## peso 0.1435076 0.02270099 0.1815163
## tiempo_en_recorrer_3km 1.0000000 0.40053599 0.3136478
## pulso_en_reposo 0.4005360 1.00000000 0.3179732
## pulso_final 0.3136478 0.31797319 1.0000000
## pulso_maximo 0.2261030 0.25750341 0.9297538
## pulso_maximo
## consumo_oxigeno -0.2367402
## edad -0.4329159
## peso 0.2493812
## tiempo_en_recorrer_3km 0.2261030
## pulso_en_reposo 0.2575034
## pulso_final 0.9297538
## pulso_maximo 1.0000000No se entiende nada, vamos a ver si con unos graficos de dispersión se capta más la idea:
aerobico %>%
gather(-consumo_oxigeno, key = "var", value = "value") %>%
ggplot(aes(x = value,color= 'red', y = consumo_oxigeno)) +
geom_point() +
facet_wrap(~ var, scales = "free") +
theme_minimal()
En este caso, pareciera haber una cierta linealidad entre el consumo de oxígeno y el tiempo en recorrer tres kilómetros.
Pruebo hacer ahora el modelo de Regresión Lineal Múltiple, usando todas las variables:
colnames(aerobico)## [1] "consumo_oxigeno" "edad" "peso"
## [4] "tiempo_en_recorrer_3km" "pulso_en_reposo" "pulso_final"
## [7] "pulso_maximo"modelo_aerobico = lm(formula = consumo_oxigeno ~ edad + peso + tiempo_en_recorrer_3km + pulso_en_reposo + pulso_final + pulso_maximo, data = aerobico)
summary(modelo_aerobico)##
## Call:
## lm(formula = consumo_oxigeno ~ edad + peso + tiempo_en_recorrer_3km +
## pulso_en_reposo + pulso_final + pulso_maximo, data = aerobico)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4699 -0.8499 0.0141 0.9996 5.3806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.022e+02 1.245e+01 8.210 1.99e-08 ***
## edad -2.199e-01 9.959e-02 -2.208 0.03703 *
## peso -7.238e-02 5.467e-02 -1.324 0.19802
## tiempo_en_recorrer_3km -2.681e+00 3.749e-01 -7.150 2.17e-07 ***
## pulso_en_reposo -8.442e-04 5.863e-02 -0.014 0.98863
## pulso_final -3.732e-01 1.207e-01 -3.092 0.00498 **
## pulso_maximo 3.047e-01 1.372e-01 2.221 0.03606 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.322 on 24 degrees of freedom
## Multiple R-squared: 0.848, Adjusted R-squared: 0.81
## F-statistic: 22.32 on 6 and 24 DF, p-value: 1.023e-08Vemos el r común (?), pero acá lo que nos interesa acá es el r2 ajustado, con 0.81, lo que estaría indicando que buena parte de la varianza de la muestra puede ser explicada a través de la varianza del modelo. en segundo lugar. Asimismo podemos ver el p-valor del análisis de la varianza, o sea ANOVA (el F-Statistics),da un valor -en fórmula científica- muy por debajo de 0,05. Lo que indica que la varianza de nuestro modelo tiene relación con la varianza muestral.
Si me da paja ver los p-valor de todas las variables predictoras puedo mandarle esta función para que me muestre sólo los que tienen p-value menor a 0,05. O sea, los estadísticamente signifcativos
names(summary(modelo_aerobico)$coefficients[,4]<0.05)## [1] "(Intercept)" "edad" "peso"
## [4] "tiempo_en_recorrer_3km" "pulso_en_reposo" "pulso_final"
## [7] "pulso_maximo"Vemos que, además del tiempo en recorrer 3 Km que habíamos visto antes hay, junto con el que ⍺ es estadísticamente significativo, otras posibles variables predictoras.
Pero, ¿realmente nos sirven?
Bueno, saco las dos variables con las que no encontré relación estadísticamente significativa (o RES de aquí en más) y ahora pruebo diferentes formas de combinación (con todas o algunas de las variables) para intentar explorar de qué manera mi modelo se ajusta mejor. No necesariamente al aumentar variables aumentan las capacidades de mi modelo, ya que puedo aumentarle imprecisión, que las variables en cuestión solas no sean significativas, bajar el el r2 , subirlo muy poco.
Es decir, mis variables descartadas no mostraban RES en presencia de las otras. Ahora, estre las escogidas pueden darse relaciones entre ellas que las que se daban antes con las descartadas.
Para esto, utilizamos un método que nos permita elegir las mejores predictoras, las mejores dentro de mi modelo. De todas formas, siempre hay una parte subjetiva al final en la que uno define según el contexto, cuáles tienen más sentido.
Selección de Modelos
A la hora de seleccionar mi modelo tengo que considerar lo siguiente:
Puede pasar que tenga varios modelos con las mismas métricas, depende de mí (de vos) con cual te quedas. No busques la perfección.
Dado que cada vez que aumento una variable -por defecto- me sube el r2 , lo que yo tengo que ver es en realidad cuándo el aumento de ese r2 es estadísticamente significativo. Y como tengo que hacer distinto tipo de comparaciones entre ellas no lo pienso hacer a mano.
Las estrategias entonces para seleccionar automáticamente mi mejor modelo son:
1- Armar (“ajusta”) todos los posibles modelos (all possible subsets) y compararlos simultáneamente entre sí -el más exhaustivo-
2- Seleccionar de a pasos (stepwise methods), agregando o quitando variables, y chequear si al colocar una nueva variable aumenta de manera significativa mi r2 .
Importante: para cualquiera de estos métodos es importante chequear la dimensionalidad de nuestro DF, ya que comparan la cantidad de predictores con las observaciones, por lo que se verán sesgados mientras menos observaciones que features tengamos.
Métodos secuenciales de selección automática de modelo
Todas pertenecientes al paquete olsrr
Forward: ols_step_forward_p(model, penter = 0.3, progress = FALSE, details = FALSE)
Backward: ols_step_backward_p(model, prem = 0.3, progress = FALSE, details = FALSE)
Stepwise : ols_step_both_p(model, pent = 0.1, prem = 0.3, progress = FALSE, details = FALSE)
Se diferencian en la forma entre que van eligiendo los predictores y los modelos.
El Forward elije la variable con más RES, y va agregando de a una, probando si se incrementa el r2 de manera estadísticamente significativa. Puede ser que ninguna de las variables a ingresar de un valor más alto del r2, hasta que no hay posibilidad de sumar alguna variable que aumente significativamente el r2.
El Backward es al revés, arranca con todas las variables y va sacando en cada paso una, la que su inclusión o exclusión no cambia significativamente el r2, etc., hasta que la que una que saca baja significativamente el r2.
El Stepwise lo que hace es ir agregando variables y viendo si cuando agrega una no tiene que sacar alguna de las que ya estaba antes. Siempre va revisando para atrás, cuando agrega una, para ver si el ingreso o el egreso hace cambiar el r2 respecto a otras combinaciones previas. Es el mejor de los de selección de a pasos.
Ahora, ¿cómo comparamos si dos o más modelos parece que ajustan bien?
Criterios para comparación de modelos
Estos serían los siguientes criterios, recomendados en este orden.
En realidad, si no puede verse con el r2, se utiliza el Cp de Mallows. Los otros dos (a aplican más bien para Regresión Logística (ya que allí no puedo usar MMC).
El Cp de Mallows es una medida importante que indica el sesgo que se introduce al predecir la variable con un modelo mal especificado. En otras palabras, cuánto le pifias vos al predecir los valores de y usando modelos que no tienen el conjunto óptimo de predictores.
La teoría dice que un modelo óptimo tiene un Cp de Mallows igual o menor a la cantidad de parámetros del modelo. Si el modelo tiene 4 predictores, los parámetros totales son 5. Entonces si mi Cp de Mallows está cerca de cinco, es un buen modelo.
Continuación del ejemplo 3
#además de la manera en que lo calculamos previamente, podemos ver las métricas del modelo usando el paquete olsrr. En este caso, definimos de manera muy similar la sintaxis nada más que en vez de lm() usamos ols_regress()
modelo_aerobico_ols <- ols_regress(consumo_oxigeno ~ edad + peso + tiempo_en_recorrer_3km + pulso_en_reposo + pulso_final + pulso_maximo, data = aerobico)
modelo_aerobico_ols## Model Summary
## -------------------------------------------------------------
## R 0.921 RMSE 2.322
## R-Squared 0.848 Coef. Var 4.901
## Adj. R-Squared 0.810 MSE 5.392
## Pred R-Squared 0.763 MAE 1.437
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 721.974 6 120.329 22.316 0.0000
## Residual 129.407 24 5.392
## Total 851.382 30
## -------------------------------------------------------------------
##
## Parameter Estimates
## -----------------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -----------------------------------------------------------------------------------------------------
## (Intercept) 102.238 12.453 8.210 0.000 76.537 127.940
## edad -0.220 0.100 -0.215 -2.208 0.037 -0.425 -0.014
## peso -0.072 0.055 -0.113 -1.324 0.198 -0.185 0.040
## tiempo_en_recorrer_3km -2.681 0.375 -0.698 -7.150 0.000 -3.454 -1.907
## pulso_en_reposo -0.001 0.059 -0.001 -0.014 0.989 -0.122 0.120
## pulso_final -0.373 0.121 -0.718 -3.092 0.005 -0.622 -0.124
## pulso_maximo 0.305 0.137 0.524 2.221 0.036 0.022 0.588
## -----------------------------------------------------------------------------------------------------El cálculo es igual que el lm() pero sale todo mucho más chetooo. Fuera de eso, existen otras funciones interesantes del paquete olsrr(), que pasaremos a ver, pero que necesitan que le pasemos como variable el modelo construido con lm(). Por ejemplo, si queremos ver sólo las correlaciones de nuestra variable predicha con las predictoras:
ols_correlations(modelo_aerobico)## Correlations
## ---------------------------------------------------------
## Variable Zero Order Partial Part
## ---------------------------------------------------------
## edad -0.305 -0.411 -0.176
## peso -0.163 -0.261 -0.105
## tiempo_en_recorrer_3km -0.862 -0.825 -0.569
## pulso_en_reposo -0.346 -0.003 -0.001
## pulso_final -0.398 -0.534 -0.246
## pulso_maximo -0.237 0.413 0.177
## ---------------------------------------------------------#Zero order es la correlacion (el r), de cada variable independiente con la variable yAhora, pasemos a ver cómo calcular los modelos, primero comparando todos los posibles modelos simultáneamente:
all_modelo_aerobico = ols_step_all_possible(modelo_aerobico)
all_modelo_aerobico # ver como hacer para que la salida me quede más legible## Index N
## 3 1 1
## 5 2 1
## 4 3 1
## 1 4 1
## 6 5 1
## 2 6 1
## 8 7 2
## 17 8 2
## 18 9 2
## 12 10 2
## 16 11 2
## 10 12 2
## 21 13 2
## 11 14 2
## 9 15 2
## 19 16 2
## 14 17 2
## 7 18 2
## 13 19 2
## 20 20 2
## 15 21 2
## 27 22 3
## 40 23 3
## 28 24 3
## 22 25 3
## 26 26 3
## 38 27 3
## 33 28 3
## 34 29 3
## 39 30 3
## 32 31 3
## 29 32 3
## 31 33 3
## 24 34 3
## 30 35 3
## 41 36 3
## 37 37 3
## 23 38 3
## 25 39 3
## 35 40 3
## 36 41 3
## 50 42 4
## 43 43 4
## 54 44 4
## 56 45 4
## 48 46 4
## 44 47 4
## 49 48 4
## 42 49 4
## 52 50 4
## 53 51 4
## 45 52 4
## 51 53 4
## 47 54 4
## 46 55 4
## 55 56 4
## 59 57 5
## 61 58 5
## 62 59 5
## 57 60 5
## 58 61 5
## 60 62 5
## 63 63 6
## Predictors
## 3 tiempo_en_recorrer_3km
## 5 pulso_final
## 4 pulso_en_reposo
## 1 edad
## 6 pulso_maximo
## 2 peso
## 8 edad tiempo_en_recorrer_3km
## 17 tiempo_en_recorrer_3km pulso_final
## 18 tiempo_en_recorrer_3km pulso_maximo
## 12 peso tiempo_en_recorrer_3km
## 16 tiempo_en_recorrer_3km pulso_en_reposo
## 10 edad pulso_final
## 21 pulso_final pulso_maximo
## 11 edad pulso_maximo
## 9 edad pulso_en_reposo
## 19 pulso_en_reposo pulso_final
## 14 peso pulso_final
## 7 edad peso
## 13 peso pulso_en_reposo
## 20 pulso_en_reposo pulso_maximo
## 15 peso pulso_maximo
## 27 edad tiempo_en_recorrer_3km pulso_final
## 40 tiempo_en_recorrer_3km pulso_final pulso_maximo
## 28 edad tiempo_en_recorrer_3km pulso_maximo
## 22 edad peso tiempo_en_recorrer_3km
## 26 edad tiempo_en_recorrer_3km pulso_en_reposo
## 38 tiempo_en_recorrer_3km pulso_en_reposo pulso_final
## 33 peso tiempo_en_recorrer_3km pulso_final
## 34 peso tiempo_en_recorrer_3km pulso_maximo
## 39 tiempo_en_recorrer_3km pulso_en_reposo pulso_maximo
## 32 peso tiempo_en_recorrer_3km pulso_en_reposo
## 29 edad pulso_en_reposo pulso_final
## 31 edad pulso_final pulso_maximo
## 24 edad peso pulso_final
## 30 edad pulso_en_reposo pulso_maximo
## 41 pulso_en_reposo pulso_final pulso_maximo
## 37 peso pulso_final pulso_maximo
## 23 edad peso pulso_en_reposo
## 25 edad peso pulso_maximo
## 35 peso pulso_en_reposo pulso_final
## 36 peso pulso_en_reposo pulso_maximo
## 50 edad tiempo_en_recorrer_3km pulso_final pulso_maximo
## 43 edad peso tiempo_en_recorrer_3km pulso_final
## 54 peso tiempo_en_recorrer_3km pulso_final pulso_maximo
## 56 tiempo_en_recorrer_3km pulso_en_reposo pulso_final pulso_maximo
## 48 edad tiempo_en_recorrer_3km pulso_en_reposo pulso_final
## 44 edad peso tiempo_en_recorrer_3km pulso_maximo
## 49 edad tiempo_en_recorrer_3km pulso_en_reposo pulso_maximo
## 42 edad peso tiempo_en_recorrer_3km pulso_en_reposo
## 52 peso tiempo_en_recorrer_3km pulso_en_reposo pulso_final
## 53 peso tiempo_en_recorrer_3km pulso_en_reposo pulso_maximo
## 45 edad peso pulso_en_reposo pulso_final
## 51 edad pulso_en_reposo pulso_final pulso_maximo
## 47 edad peso pulso_final pulso_maximo
## 46 edad peso pulso_en_reposo pulso_maximo
## 55 peso pulso_en_reposo pulso_final pulso_maximo
## 59 edad peso tiempo_en_recorrer_3km pulso_final pulso_maximo
## 61 edad tiempo_en_recorrer_3km pulso_en_reposo pulso_final pulso_maximo
## 62 peso tiempo_en_recorrer_3km pulso_en_reposo pulso_final pulso_maximo
## 57 edad peso tiempo_en_recorrer_3km pulso_en_reposo pulso_final
## 58 edad peso tiempo_en_recorrer_3km pulso_en_reposo pulso_maximo
## 60 edad peso pulso_en_reposo pulso_final pulso_maximo
## 63 edad peso tiempo_en_recorrer_3km pulso_en_reposo pulso_final pulso_maximo
## R-Square Adj. R-Square Mallow's Cp
## 3 0.74338010 0.7345311403 13.519765
## 5 0.15838344 0.1293621761 105.889559
## 4 0.11999670 0.0896517619 111.950747
## 1 0.09277653 0.0614929635 116.248757
## 6 0.05604592 0.0234957795 122.048447
## 2 0.02648849 -0.0070808735 126.715506
## 8 0.76424693 0.7474074256 12.224935
## 17 0.76142381 0.7443826564 12.670699
## 18 0.74522106 0.7270225683 15.229081
## 12 0.74493479 0.7267158412 15.274283
## 16 0.74338145 0.7250515578 15.519551
## 10 0.37599543 0.3314236752 73.529064
## 21 0.28941948 0.2386637282 87.199232
## 11 0.25998174 0.2071232889 91.847392
## 9 0.24760963 0.1938674585 93.800923
## 19 0.21215907 0.1558847212 99.398495
## 14 0.16685536 0.1073450310 106.551859
## 7 0.15063534 0.0899664306 109.112969
## 13 0.14399971 0.0828568276 110.160721
## 20 0.14331054 0.0821184329 110.269540
## 15 0.06751590 0.0009098958 122.237360
## 27 0.81109446 0.7901049580 6.827804
## 40 0.80998844 0.7888760478 7.002442
## 28 0.78173017 0.7574779672 11.464366
## 22 0.77083060 0.7453673364 13.185386
## 26 0.76562507 0.7395834084 14.007329
## 38 0.76227983 0.7358664753 14.535536
## 33 0.76182904 0.7353655984 14.606715
## 34 0.74615485 0.7179498380 17.081637
## 39 0.74526986 0.7169665163 17.221375
## 32 0.74494195 0.7166021675 17.273152
## 29 0.43845650 0.3760627779 65.666587
## 31 0.42227346 0.3580816224 68.221856
## 24 0.40912553 0.3434728125 70.297888
## 30 0.35682183 0.2853575841 78.556538
## 41 0.32686542 0.2520726885 83.286594
## 37 0.32077932 0.2453103606 84.247576
## 23 0.30753079 0.2305897655 86.339493
## 25 0.29021246 0.2113471762 89.074022
## 35 0.22232446 0.1359160617 99.793401
## 36 0.15778804 0.0642089332 109.983571
## 50 0.83681815 0.8117132468 4.766086
## 43 0.81649255 0.7882606383 7.975456
## 54 0.81584902 0.7875181054 8.077068
## 56 0.81162125 0.7826399053 8.744625
## 48 0.81115901 0.7821065501 8.817612
## 44 0.78622430 0.7533357276 12.754753
## 49 0.78244195 0.7489714787 13.351978
## 42 0.77297910 0.7380528027 14.846143
## 52 0.76260533 0.7260830738 16.484140
## 53 0.74617451 0.7071244322 19.078534
## 45 0.47593322 0.3953075564 61.749089
## 51 0.47235647 0.3911805443 62.313850
## 47 0.47171966 0.3904457653 62.414401
## 46 0.39278110 0.2993628100 74.878640
## 55 0.35917422 0.2605856441 80.185098
## 59 0.84800181 0.8176021762 5.000207
## 61 0.83690359 0.8042843054 6.752595
## 62 0.81712216 0.7805465875 9.876043
## 57 0.81677090 0.7801250814 9.931505
## 58 0.78744812 0.7449377489 14.561513
## 60 0.52421007 0.4290520785 56.126272
## 63 0.84800313 0.8100039082 7.000000Aquí me salen todos los posibles modelos (63 en este caso) ordenados de mayor r2 a menor r2 . Pensar que con pocas variables esto puede funcionar pero con muchos datos y muchas variables puede llegar a demandar mucho capacidad de procesamiento.
En ese sentido, la función ols_step_best_subset(“modelo”), nos da el mejor modelo para cada k cantidad de variabes.
ols_step_best_subset(modelo_aerobico)## Best Subsets Regression
## ----------------------------------------------------------------------------------------
## Model Index Predictors
## ----------------------------------------------------------------------------------------
## 1 tiempo_en_recorrer_3km
## 2 edad tiempo_en_recorrer_3km
## 3 edad tiempo_en_recorrer_3km pulso_final
## 4 edad tiempo_en_recorrer_3km pulso_final pulso_maximo
## 5 edad peso tiempo_en_recorrer_3km pulso_final pulso_maximo
## 6 edad peso tiempo_en_recorrer_3km pulso_en_reposo pulso_final pulso_maximo
## ----------------------------------------------------------------------------------------
##
## Subsets Regression Summary
## ---------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ---------------------------------------------------------------------------------------------------------------------------------
## 1 0.7434 0.7345 0.7053 13.5198 154.5083 65.4840 158.8103 233.5838 8.0199 0.2691 0.2920
## 2 0.7642 0.7474 0.7054 12.2249 153.8792 64.8457 159.6151 222.5380 7.8621 0.2655 0.2863
## 3 0.8111 0.7901 0.7591 6.8278 149.0115 61.3486 156.1814 185.1748 6.7253 0.2291 0.2449
## 4 0.8368 0.8117 0.7785 4.7661 146.4737 60.4483 155.0776 166.3575 6.2053 0.2137 0.2259
## 5 0.8480 0.8176 0.7867 5.0002 146.2728 61.6283 156.3107 161.4127 6.1782 0.2157 0.2250
## 6 0.8480 0.8100 0.7628 7.0000 148.2725 64.2115 159.7444 168.4292 6.6095 0.2344 0.2407
## ---------------------------------------------------------------------------------------------------------------------------------
## 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 CriteriaY ahora, probamos ahora el mejor método para seleccionar (Stepwise):
ols_step_both_p(modelo_aerobico)##
## Stepwise Selection Summary
## ---------------------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## ---------------------------------------------------------------------------------------------------
## 1 tiempo_en_recorrer_3km addition 0.743 0.735 13.5200 154.5083 2.7448
## ---------------------------------------------------------------------------------------------------Aquí vemos que nos selecciona como el mejor modelo el que tiene sólo la variable que antes habíamos encontrado, gráficamente, como estadísticamente representativa (el tiempo_en_recorrer_3km). Es decir, sólo con esa variable puedo predecir de manera óptima, considerando las otras variables que tengo.
¿Cómo continuamos? Antes de quedarse simplemente con el modelo que seleccionó el Stepwise, uno podría fijarse que pasa, en este caso, con el mejor modelo antes encontrado para dos variables, y chequear con un summary de dicho modelo, cuán estadísticamente significativas resultan los predictores juntos.
Entonces, en el caso anterior nos había señalado que el mejor modelo con dos variables sería con tiempo_en_recorrer_3km y edad. Armemos para probar una lm() y veamos el summary:
summary(lm(consumo_oxigeno~tiempo_en_recorrer_3km + edad, data = aerobico))##
## Call:
## lm(formula = consumo_oxigeno ~ tiempo_en_recorrer_3km + edad,
## data = aerobico)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8744 -1.5443 0.1273 1.5575 4.9567
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 88.46229 5.37264 16.465 6.21e-16 ***
## tiempo_en_recorrer_3km -3.20395 0.35877 -8.930 1.10e-09 ***
## edad -0.15037 0.09551 -1.574 0.127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.677 on 28 degrees of freedom
## Multiple R-squared: 0.7642, Adjusted R-squared: 0.7474
## F-statistic: 45.38 on 2 and 28 DF, p-value: 1.638e-09Como podemos ver, en este caso, la edad junto con el tiempo… no es estadísticamente significativa, por lo que nos quedamos con nuestro modelo seleccionado por el método stepwise.
Multicolinealidad
¿Por qué variables tan altamente correlacionadas con la variable dependiente no son consideradas buenos predictores en mi modelo?
Esas variables con alta correlación no sólo tienen una alta correlación con la dependiente sino que tienen alta correlación entre sí. Y en realidad, lo que hacen, es explicar “lo mismo”11 . Es decir una de las variables predictoras ya es suficiente -agregar la otra, sólo incrementa mi dispersión.
Cuando las variables dependientes están altamente correlacionadas se habla del fenómeno de colinealidad (si son sólo dos) o multicolinealidad.
Debido a los problemas mencionados, es necesario analizar si existe la multicolinealidad entre variables en nuestro modelo para evitar los efectos adversos que produce (limita el r2 , dificulta la determinación de la importancia de los predictores, incremento de la varianza de los estimadores)
Hay diferentes formas de diagnosticar esto (correlación entre predictores), pero los principales son:
Variance Inflation Factor (ó VIF): indica cuándo un predictor tiene una fuerte asociación lineal con el resto de los predictores. Un VIF > 10, indica multicolinealidad. Su fórmula es:
Tolerancia (ó TOL): opera de la manera inversa (1/VIF). Un TOL < 0.10 indica multicolinealidad.
Estos indicadores se calculan para cada una de las variables predictoras que uno está incluyendo en el modelo. Lo que se hace es calcular para cada variable su regresión respecto de los otros predictores tenidos en cuenta. Ese r2, se coloca en el VIF y se obtiene la métrica para cada variable.
La búsqueda es sacar esos predictores que pueden ser predichos por otra variable independiente.
La función que nos permite apreciar esto es ols_vif_tol(“modelo”)
ols_vif_tol(modelo_aerobico)## Variables Tolerance VIF
## 1 edad 0.6672145 1.498768
## 2 peso 0.8668311 1.153627
## 3 tiempo_en_recorrer_3km 0.6643874 1.505146
## 4 pulso_en_reposo 0.7599631 1.315853
## 5 pulso_final 0.1174186 8.516537
## 6 pulso_maximo 0.1136534 8.798680En este caso tanto la TOL como VIF pareciera dar bien para nuestras variables, así que no estamos en presencia de multcolinealidad (TOL>0,10 y VIF <10 en todos los casos). Si fuera de otra manera, deberíamos quitar las variables que presentan multicolinealidad.12
Valores atípicos y puntos de influencia
Los modelos de Regresión tienen una debilidad: son altamente vulnerables a la presencia de valores atípicos ¿Qué es un valor atípico o outlier?
Es un valor que tiene una probabilidad muy baja de ocurrir, de ser obtenido. En ciertos casos, puede deberse a un error en la toma de la muestra. Estos valores deben ser detectados y analizados en forma puntual ¿por qué? Porque la estadística es una disciplina cuantitativa, por lo que todos los métodos tiene por objetivo la búsqueda, interpretación y explicación de las regularidades poblacionales (no lo irregular).
Debido a que el método de regresión es un procedimiento de maximización, se ve fuertemente influenciado por la presencia de casos atípicos; de hecho, pocos valores atípicos pueden producir severos trastornos en los modelos estimados. Veamos por ejemplo, esta situación:
¿Cómo afecta a nuestro modelo este tipo de outliers?
Una media se ve desplazada por el valor atípico, deja de apuntar al centro de la distribución.
Una varianza/desvío estándar son muy afectades por el valor atípico
Los outliers pueden estar presentes tanto en las variable dependiente como en las independientes. Para analizarlos, es necesario primero definir sobre qué tipo de variable estamos haciendo el diagnóstico.
Para las variables dependientes analizaremos los residuos.
Para las variables dependientes, o predictores, se analizarán los puntos de influencia.
¿Qué hacer con los valores atípicos?
Existen dos maneras de trabajar los valores atípicos:
Detección y eliminación de los mismos.
Utilización de Regresiones Robustas: generar modificaciones en las técnicas de estimación de los parámetros de la regresión a los efectos de que las mismas no se vean influenciadas por los outliers.
En la práctica, sin embargo, la regresión robusta no es de la más utilizada, por la complejidad de sus métodos.
Valores atípicos en la variable dependiente
Supongamos que tenemos el siguiente gráfico de residuos (o sea, de la variable dependiente)
Como se ve en el gráfico previo, un sólo outlier hace que nuestro modelo pierda ajuste (vemos como se incrementa el r2 en el segundo gráfico frente al primero)
Valores atípicos en la variable independiente
Más importante son los variables atípicos en las variables independientes. En la siguiente imagen, vemos como en el segundo gráfico todos los valores de x caen dentro del rango, con la excepción del último -alejado a la derecha y arriba. Eso podría llegar a ser un outlier. Ahora, es necesario analizar su valor de y, porque quizá -si bien está fuera de los valores del rango- se adecúa a la recta. Como se ve en el ejemplo, la recta azul incluye dicho dato y la roja, no. Sin embargo, el valor la pendiente y el intercepto parecieran ser similares, por lo que no se trata de un ouliter.
Sin embargo, en el gráfico de arriba, el valor de y asociado a ese x alejado es distinto y vemos cómo cambia la pendiente según la recta de regresión de MMC, justamente porque para minimizar el cuadrado hace que la recta cambie la pendiente. Por eso, se denomina a estos puntos como “puntos de influencia” (o leverage points)
En resumen: pueden haber valores atípicos entre los predictores que no generen modificaciones sustanciales a la pendiente de la recta y el ajuste del modelo. Pero, un valor atípico en las variables independientes también puede llegar a producir un cambio notorio en el modelo.
¿Cómo hacemos para identificarlos y trabajar con ellos?
Pues a través del criterio de ajustar los parámetros de un modelo con y sin los casos, para ver si los mismos tienen influencia en la estimación de los parámetros. En otras palabras, un punto no influyente debería producir diferencias mínimas en la estimación de los parámetros. Para observar teóricamente esto, se recurre a la parte de álgebra matricial del modelo de regresión, que no voy a abordar.
¿Entonces qué hacemos?
Cumplimiento de supuestos y análisis de los residuos
Siempre debemos chequear los supuestos (antes de hacer cualquier cosa en tu vida) pero sobre todo, antes de armar un Modelo de Regresión Lineal Múltiple. Los mismos son:
Y además! Que no se vea afectado por outliers!
Así que para esto, hacemos un análisis de los residuos
Así que vamos a hacer unos graficos y métricas y esas cosas para ver si se cumple o tengo que trabajar sobre los valores atípicos. Para esto, tenemos que hacer una estandarización de los residuos (ese hii es lo que es saca de algebra matricial que no explique) y evaluar su varianza (si hay un outlier aumenta mucho la varianza).
Entonces lo que se hace es estandarizar el residuo pero con un estimador del desvío estándar que se calcule utilizando todos los datos excepto aquel individuo cuyo residuo está anaizando ¿Con esto qué me evito? Me evito que si ese individuo tiene un outlier, quede incluido en el calculo de la variabilidad. Entonces voy a tener una variabilidad no influida por ese caso.
El análisis de residuos más recomendado es el Residuo R de Student.
Chequeo de supuestos: gráficos de residuos para chequeo de linealidad, normalidad y homocedasticidad.
¿Cómo se grafican y analizan los residuos? De la siguiente forma: en el eje y, de ordenadas, los residuos y en las abscisas:
a) Si tengo un único predictor, grafico residuos contra el predictor
b) Si tengo una regresión con más de un predictor, mi gráfico tiene que tener los residuos en la ordenada, mientras que en las abcisas colocar el y_hat (el predictor de y).
Con estas consideraciones, si se cumplen los supuestos de la regresión, el gráfico de los residuos con varianza constante para un modelo lineal deberia verse de esta manera para el caso a) de un sólo predictor.
y de esta manera para el caso b), con más de un predictor:
Si está todo bien (es decir, se cumplen nuestros supuestos de linealidad y homocedasticidad), los residuos van a quedar dentro de un rango de -3 y 3, centrado en 0. No debería haber irregularidades que nos haga pensar que la relación con los predictores no es lineal o puntos fuera de la banda (valores atípicos).
Si fuese que los residuos dan una forma como la que se ve en los cuadrantes izquierdos
Significaría que la relación no es lineal, y habría que pensar en operaciones sobre nuestras variables para linealizarlas, agregar predictores, agregar la multiplicación de los mismos, etc. (exceden este apunte). Ahora bien, en el gráfico anterior no se estaría cumpliendo el supuesto de linealidad, pero aparentemente sí el de homocedasticidad.
¿Cómo se ve entonces una varianza de residuos no homogénea? De la manera que aparece del lado izquierdo de la siguiente placa:
Algo así como un “efecto embudo”, no de forma paralela (como el gráfico de la esquina inferior derecha). Para resolverlo, en algunos casos, puede realizarse alguna transformación de los datos (por ejemplo, un logaritmo).
Por último, el tema de la independencia de las observaciones está generalmente más vinculado al trabajo con series temporales:
No lo veremos particularmente tampoco.
Otros supuestos: análisis y corrección de valores atípicos y puntos de influencia en base a la distancia de Cook
Como se dijo previamente, uno de los supuestos que debemos chequear antes de darle el ok a nuestro modelo es que no se vea significativamente influenciado por valores atípicos o puntos de influencia. Para chequear esto, utilizamos la distancia de Cook.
Como vemos, nos da información del residuo (variable y) pero también si hay también algún valor alocado entre las predictoras x (posibles puntos de influencia). El valor de Cook, que opera sobre los casos individuales, nos sirve para ubicar valores atípicos:
Otra posibilidad son DFBetas
El DFBeta calcula y compara el estimador del parámetro con el caso y sin el caso (estandarizado obviamente). Si hay un outlier, el estimador del parámetro tiene que cambiar mucho en la comparación. Nos identifica los puntos de influencia.
El DFits es lo mismo, pero analiza el valor de la variable predicha (o sea, los residuos y valores atípicos)
Para ambos hay valores de corte, que permiten determinar cuáles son los individuos que tienen quizá situaciones que deban ser analizadas puntualmente, porque podrían estar teniendo una influencia no deseada en el ajuste del modelo.
Continuación Ejemplo 3
Gráfico de los residuos
Este sería el Residuo del Rango Studentizado
ols_plot_resid_stud_fit(modelo_aerobico, print_plot = TRUE) # Esta es la mejor forma de graficarlos, después hay otras también: ols_plot_resid_fit(model, print_plot = TRUE), ols_plot_resid_stand(model, print_plot = TRUE)
Si la distribución es normal, tenemos el 95% de las variaciones a lo largo de la curva de -2 y 2 (como vemos, solo tres puntos quedan afuera). Quizá esos puntos rojos, si son más o mayores debería ver en particular esos casos.
Gráfico que detecta outliers y puntos influyentes
Para ver los casos de los valores de x que puedan tener cierta influencia en el ajuste del modelo, se calcula el leverage (punto de influencia, ese hii que no explique pero te calcula valores atípicos en x, básicamente)
ols_plot_resid_lev(modelo_aerobico, print_plot = TRUE)
#Este es otro ols_leverage(model)En este caso, sería el valor 10, como vemos.
Gráficos de la distancia de Cook
ols_plot_cooksd_bar(modelo_aerobico, print_plot = TRUE)
Esos pueden llegar a ser outliers (y) o leverage points (x) que puede estar teniendo impacto en la estimación. Coincide con el análisis de los residuos. Es otra forma de representar lo anterior.
DFBeta y DFfits
DFfits es la comparación valor predicho y incluyendo y excluyendo el caso:
ols_plot_dffits(modelo_aerobico, print_plot = TRUE)
Acá lo que dice es que el predictor de y cambia mucho si los casos 10, 15 y 20, están o no están (el resto quedan dentro de
DFBetas es la comparación de los estimadores de los parámetros (𝛃) incluyendo y excluyendo los casos:
ols_plot_dfbetas(modelo_aerobico, print_plot = TRUE)
Podemos ver tanto la influencia entre los Betas como en Alpha (el intercepto). Podemos interpretar que los casos 4 (para el intercepto), 10, 15 y 20 (para los beta) están teniendo influencia en la configuración de los estimadores de los parámetros.
¿Ahora qué haría en este caso? Pues bien, parecería ser que, entre todos, el 10, 15 y el 20 son los que más influencia están teniendo. Podemos
Correr una regresión robusta (hoy no mi ciela)
O excluir esos casos y ver que pasa con el ajuste, si mejora o no.
De todas formas, puede ocurrir que la exclusión, de todas formas, no haga variar demasiado el modelo.
Hay que tener en cuenta que si saco una observación, puede que nuevamente me aparezca otra que crea que es valor atípico ¿Entonces me la paso sacando observaciones que simulan ser outliers para toda la eternidad? No! Porque así llegás al superajuste (u overfiting), que lo que hace es que tu modelo sea genia… pero sólo para esa muestra.
Nosotros necesitamos modelos que sean aplicables a cualquier muestra de la población objetivo que tengamos que analizar. Hay que ver si se justifica mucho sacar esos casos, medir cuánto cambia el r2 sacando ese caso (si es 0,0algo ni me gasto), o se hace la cross-validation: el famoso test and train!
Voy a dejar para después probar los residuos del modelo sacando esas observaciones. Para eso tengo que modificar el dataframe primero, después tirarle el lm() y luego el Coso de Cook o alguno de los DFBets o DFfits.
aerobico_sin_outliers = aerobico %>% slice(-c(10, 15, 20),)Y ahora calculamos de nuevo el lm()
modelo_aerobico_sin_out= lm(consumo_oxigeno ~ edad + peso + tiempo_en_recorrer_3km + pulso_en_reposo + pulso_final + pulso_maximo, data = aerobico_sin_outliers)
summary(modelo_aerobico_sin_out)##
## Call:
## lm(formula = consumo_oxigeno ~ edad + peso + tiempo_en_recorrer_3km +
## pulso_en_reposo + pulso_final + pulso_maximo, data = aerobico_sin_outliers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9407 -0.9089 0.3082 0.9311 2.7368
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 107.69268 10.47799 10.278 1.20e-09 ***
## edad -0.29586 0.08758 -3.378 0.00284 **
## peso -0.10266 0.04627 -2.219 0.03766 *
## tiempo_en_recorrer_3km -2.29987 0.33530 -6.859 8.84e-07 ***
## pulso_en_reposo -0.06778 0.06204 -1.093 0.28698
## pulso_final -0.35335 0.12717 -2.779 0.01126 *
## pulso_maximo 0.28332 0.14998 1.889 0.07278 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.913 on 21 degrees of freedom
## Multiple R-squared: 0.8812, Adjusted R-squared: 0.8473
## F-statistic: 25.97 on 6 and 21 DF, p-value: 1.099e-08ols_step_both_p(modelo_aerobico_sin_out)##
## Stepwise Selection Summary
## ---------------------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## ---------------------------------------------------------------------------------------------------
## 1 tiempo_en_recorrer_3km addition 0.758 0.749 18.7620 133.6354 2.4530
## 2 pulso_final addition 0.787 0.770 15.5880 132.0245 2.3454
## 3 edad addition 0.839 0.819 8.4500 126.2248 2.0825
## 4 peso addition 0.858 0.833 7.1470 124.7697 2.0000
## ---------------------------------------------------------------------------------------------------ols_step_both_p(modelo_aerobico)##
## Stepwise Selection Summary
## ---------------------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## ---------------------------------------------------------------------------------------------------
## 1 tiempo_en_recorrer_3km addition 0.743 0.735 13.5200 154.5083 2.7448
## ---------------------------------------------------------------------------------------------------Y bueno, ahora sacando los outliers pareciera que tengo cuatro variables estadísticamente significativas y me aumenta el r2 del modelo en 0,1. Habría que testear un poco más y ver cuán estadísticamente significativas son cada una de esas cuatro variables dentro del modelo nuevo. Pero yo a priori me quedaría con el otro que es más simple y no es tanta la diferencia del r2.
Y ahora, ejercicios:
Ejercicio 1
Por medio de este ejercicio se propone evaluar la capacidad de análisis de asociación y de predicción entre variables cuantitativas.
Para este fin se utilizará el conjunto de datos que se describe a continuación: Archivo CORDOBA.SAV
El archivo CORDOBA.SAV contiene información de algunas características de vivienda, hogar y población de departamentos de la provincia de Córdoba. Este archivo consta de las siguientes variables13:
DEPART: Departamento provincial
NBI1: Porcentaje de hogares con Necesidades Básicas Insatisfechas sobre el total de hogares de cada departamento.
NBI2: Porcentaje de población en hogares con Necesidades Básicas Insatisfechas sobre el total de población de cada departamento.
CALMAT1: porcentaje de viviendas que presentan materiales resistentes y sólidos en todos los paramentos (pisos, paredes o techos) e incorpora todos los elementos de aislación y terminación.
CALMAT2: porcentaje de viviendas que presentan materiales resistentes y sólidos en todos los paramentos pero le faltan elementos de aislación o terminación al menos en uno de sus componentes (pisos, paredes, techos).
CALMAT3: porcentaje de viviendas que presentan materiales resistentes y sólidos en todos los paramentos pero le faltan elementos de aislación o terminación en todos sus componentes, o bien presenta techos de chapa de metal o fibrocemento.
CALMAT4: porcentaje de viviendas que presentan materiales no resistentes ni sólidos o de desecho al menos en uno de los paramentos.
ALFAB: Porcentaje de alfabetos
COB: Porcentaje de personas que tienen Obra social y/o plan de salud privado o mutual.
urlfile_1 = 'https://github.com/oblitterator/estadistica_II_UNTREF/blob/main/cordoba.sav?raw=true'
cordoba_nbi = read_sav(url(urlfile_1))
cordoba_nbi## # A tibble: 26 × 9
## depart nbi1 nbi2 calmat1 calmat2 calmat3 calmat4 alfab cob
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Calamuchita 13.7 15.1 47.5 40.2 9.15 3.21 97.1 49.9
## 2 Capital 9.6 12.2 57.3 35.8 6 0.85 98.7 55.3
## 3 Colón 13.1 16.4 41.2 45.6 11.2 2.02 98.0 49.3
## 4 Cruz del Eje 24.6 29.2 32.8 45.4 10.6 11.1 94.8 37.2
## 5 General Roca 13.8 15.3 59.7 28.7 9.61 1.96 96.6 53.3
## 6 General San Martín 9.4 10 53.1 40.1 5.97 0.82 98.0 60.3
## 7 Ischilín 18.1 21.1 38.5 46.8 9.81 4.9 96.2 42.7
## 8 Juárez Celman 10.6 9.8 49.0 42.7 7.37 0.95 97.2 62.5
## 9 Marcos Juárez 9.5 8.8 53.3 40.1 5.83 0.79 97.7 62.3
## 10 Minas 35.9 39.5 19.3 41.6 12.7 26.3 95.4 33.8
## # … with 16 more rowscolnames(cordoba_nbi)## [1] "depart" "nbi1" "nbi2" "calmat1" "calmat2" "calmat3" "calmat4"
## [8] "alfab" "cob"Se requiere realizar los siguientes análisis estadísticos sobre este conjunto de variables:
1. Realizar un análisis exploratorio de todas las variables de interés.
a. Mediante los gráficos de Boxplot e Histograma, evaluar simetría, presencia de observaciones atípicas, dispersión de cada variable.
b. Calcular indicadores de tendencia central y de dispersión, indicando cuál/cuáles son los más adecuados a las características de las variables.
2. Realizar gráficos de dispersión de la variable NBI con cada una de las restantes variables.
3. Calcular correlaciones lineales de la variable NBI con cada una de las variables disponibles.
4. Calcular la recta de regresión utilizando como variable dependiente el NBI y como independiente cada una de las restantes variables
1. Realizar un análisis exploratorio de todas las variables de interés.
a. Mediante los gráficos de Boxplot e Histograma, evaluar simetría, presencia de observaciones atípicas, dispersión de cada variable.
a. Boxplot e Histogramas
Primero los Boxplot:
melt_cba = melt(cordoba_nbi, id.var = 'depart')## Warning: attributes are not identical across measure variables; they will be
## droppedboxplot(value ~ variable, data = melt_cba)
Ahora los histogramas:
#Esta función pertenece al paquete Hmisc
#Revisar por q rompe
hist.data.frame(x = cordoba_nbi[-1], nclass = 26) # con [-1] empiezo el calculo desde la segunda varaible para que no me grafique las etiquetas de los departamentos
b. Calcular indicadores de tendencia central y de dispersión, indicando cuál/cuáles son los más adecuados a las características de las variables.
b. Indicadores de tendencia central y de dispersión
summary(cordoba_nbi[-1])## nbi1 nbi2 calmat1 calmat2
## Min. : 9.00 Min. : 8.80 Min. :10.70 Min. :27.38
## 1st Qu.:10.53 1st Qu.:10.72 1st Qu.:32.58 1st Qu.:40.10
## Median :13.70 Median :15.70 Median :46.23 Median :41.95
## Mean :17.01 Mean :18.77 Mean :40.80 Mean :43.00
## 3rd Qu.:21.48 3rd Qu.:25.23 3rd Qu.:52.43 3rd Qu.:46.50
## Max. :39.90 Max. :40.80 Max. :64.98 Max. :62.75
## calmat3 calmat4 alfab cob
## Min. : 5.030 Min. : 0.5600 Min. :91.66 Min. :29.93
## 1st Qu.: 6.400 1st Qu.: 0.9675 1st Qu.:95.66 1st Qu.:37.26
## Median : 9.380 Median : 2.0950 Median :97.08 Median :49.62
## Mean : 9.364 Mean : 6.8381 Mean :96.47 Mean :47.77
## 3rd Qu.:11.645 3rd Qu.:10.7250 3rd Qu.:97.51 3rd Qu.:57.55
## Max. :16.100 Max. :34.8900 Max. :98.67 Max. :62.532. Gráficos de dispersión
2. Realizar gráficos de dispersión de la variable NBI con cada una de las restantes variables.
cordoba_nbi %>%
gather(-nbi1, -nbi2, key = "var", value = "value") %>%
ggplot(aes(x = value,color= 'red', y = nbi1)) +
geom_point() +
facet_wrap(~ var, scales = "free") +
theme_minimal() #le saco nbi2 porque son multicolineales## Warning: attributes are not identical across measure variables;
## they will be dropped
3. Correlaciones lineales
3. Calcular correlaciones lineales de la variable NBI con cada una de las variables disponibles.
Coeficiente de Correlación (r)
cor(cordoba_nbi[-1])## nbi1 nbi2 calmat1 calmat2 calmat3 calmat4
## nbi1 1.0000000 0.9885581 -0.8924089 0.3893973 0.6928712 0.9644483
## nbi2 0.9885581 1.0000000 -0.8862140 0.3798820 0.7446440 0.9436909
## calmat1 -0.8924089 -0.8862140 1.0000000 -0.7458662 -0.8127386 -0.8096014
## calmat2 0.3893973 0.3798820 -0.7458662 1.0000000 0.5955555 0.2370857
## calmat3 0.6928712 0.7446440 -0.8127386 0.5955555 1.0000000 0.5483462
## calmat4 0.9644483 0.9436909 -0.8096014 0.2370857 0.5483462 1.0000000
## alfab -0.8165181 -0.8076804 0.9034763 -0.6971058 -0.8078057 -0.6863838
## cob -0.8794067 -0.9190894 0.8370384 -0.3913652 -0.8007534 -0.8301653
## alfab cob
## nbi1 -0.8165181 -0.8794067
## nbi2 -0.8076804 -0.9190894
## calmat1 0.9034763 0.8370384
## calmat2 -0.6971058 -0.3913652
## calmat3 -0.8078057 -0.8007534
## calmat4 -0.6863838 -0.8301653
## alfab 1.0000000 0.7524168
## cob 0.7524168 1.0000000#ols_correlations(modelo_cba_nbi) es otra posibilidad, más linda4. Recta de Regresión
4. Calcular la recta de regresión utilizando como variable dependiente el NBI y como independiente cada una de las restantes variables
modelo_cba_nbi = lm(formula = nbi1 ~ calmat1 + calmat2 + calmat3 + calmat4 + alfab + cob, data = cordoba_nbi)
#modelo_cba_nbi = ols_regress(formula = nbi1 ~ calmat1 + calmat2 + calmat3 + calmat4 + alfab + cob, data = cordoba_nbi)
#ols_vif_tol(modelo_cba_nbi) para ver multicolonealidad entre variables predictorasHago el summary de mi modelo:
summary(modelo_cba_nbi)##
## Call:
## lm(formula = nbi1 ~ calmat1 + calmat2 + calmat3 + calmat4 + alfab +
## cob, data = cordoba_nbi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2743 -0.5281 0.0653 0.4753 2.1315
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5588.52605 5375.98653 1.040 0.312
## calmat1 -54.73675 53.73513 -1.019 0.321
## calmat2 -54.75060 53.75323 -1.019 0.321
## calmat3 -54.51372 53.77170 -1.014 0.323
## calmat4 -54.00628 53.73910 -1.005 0.328
## alfab -1.06366 0.41079 -2.589 0.018 *
## cob -0.03639 0.07081 -0.514 0.613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.392 on 19 degrees of freedom
## Multiple R-squared: 0.9806, Adjusted R-squared: 0.9745
## F-statistic: 160.4 on 6 and 19 DF, p-value: 3.1e-15Según los p-value obtenidos, el grado de alfabetismo (p_value = 0,018) sería la variable estadísticamente significativa para predecir el índice de NBI en Córdoba, con un Coeficiente de Determinación Ajustado (R2 ajustado) del 0,9745. Aclaración: podría hacer esto con la función ols_regress() del paquete olsrr, pero voy a utilizarlo sólo para ir ajustando mi modelo en este caso.
Para esto, voy a utilizar el método stepwise:
ols_step_both_p(modelo_cba_nbi)##
## Stepwise Selection Summary
## -------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------------------
## 1 calmat4 addition 0.930 0.927 46.5370 122.1731 2.3518
## 2 alfab addition 0.975 0.973 4.2240 97.1323 1.4282
## 3 calmat3 addition 0.979 0.977 2.1540 94.3497 1.3320
## -------------------------------------------------------------------------------------Entiendo que, en función de esto, al tomar esas tres variables mi modelo predice mejor las NBI (llega a su óptimo R2), así que lo que debería hacer es comporar mis distintos modelos con cada una de esas variables predictoras y volver a construir el mismo.
modelo_final = lm(formula = nbi1 ~ calmat4 + alfab + calmat3, data = cordoba_nbi)
summary(modelo_final)##
## Call:
## lm(formula = nbi1 ~ calmat4 + alfab + calmat3, data = cordoba_nbi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.06711 -0.54161 0.02508 0.52491 2.28025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.18838 31.11206 3.574 0.00170 **
## calmat4 0.74971 0.04117 18.211 9.36e-15 ***
## alfab -1.05915 0.31091 -3.407 0.00253 **
## calmat3 0.30670 0.14551 2.108 0.04667 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.332 on 22 degrees of freedom
## Multiple R-squared: 0.9795, Adjusted R-squared: 0.9767
## F-statistic: 349.7 on 3 and 22 DF, p-value: < 2.2e-16Y, finalmente, grafico:
plot(modelo_final)
Ejercicio 3
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 aquí):
| Variable | Descripción |
|---|---|
| FTP | Cantidad de policias full time cada 100.000 habitantes |
| UEMP | Porcentaje de desempleados |
| M | Cantidad de trabajadores en la Industria Manufacturera (en miles) |
| LIC | Cantidad de licencias de portación de armas cada 100.000 habitantes |
| GR | Cantidad de armas de fuego registradas cada 100.000 habitantes |
| CLEAR | Porcentaje de homicidios resueltos con arresto del responsable |
| W | Cantidad de hombres blancos en la población |
| NMAN | Cantidad de trabajadores fuera de la industria manufacturera |
| G | Cantidad de trabajadores del gobierno |
| HE | Ingreso promedio por hora |
| WE | Ingreso promedio semanal |
| H | Tasa de homicidios cada 100.000 habitantes |
Se requiere construir una ecuación de regresión que relacione la Tasa de Homicidios (variable H), con el resto de las variables disponibles, y determinar si estas variables son útiles para predecir la Tasa de Homicidios.
Antes que nada, levanto mi Data Frame
urlfile_2 = "https://github.com/oblitterator/estadistica_II_UNTREF/blob/main/p301.sav?raw=true"
crimenes_detroit = read_sav(url(urlfile_2))
crimenes_detroit = data.frame(crimenes_detroit)
colnames(crimenes_detroit)## [1] "YEAR" "FTP" "UNEMP" "M" "LIC" "GR" "CLEAR" "W" "NMAN"
## [10] "G" "HE" "WE" "H"Realizo unos ploteos para ver relación entre variables:
crimenes_detroit$YEAR = as.factor(crimenes_detroit$YEAR)
melt_detroit = melt(crimenes_detroit, id.var = 'YEAR')## Warning: attributes are not identical across measure variables; they will be
## droppedboxplot(value ~ variable, data = melt_detroit)
Y ahora calculo las correlaciones:
modelo_detroit = lm(formula = H ~ FTP + UNEMP + M + LIC + GR + CLEAR + W + NMAN + G + HE + WE, data = crimenes_detroit)
summary(modelo_detroit)##
## Call:
## lm(formula = H ~ FTP + UNEMP + M + LIC + GR + CLEAR + W + NMAN +
## G + HE + WE, data = crimenes_detroit)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.004254 -0.037388 0.252667 -0.174622 -0.169070 0.170520 -0.043951 0.036741
## 9 10 11 12 13
## -0.041456 -0.003869 -0.044030 0.003702 0.055010
## attr(,"label")
## [1] "H"
## attr(,"format.spss")
## [1] "F9.2"
## attr(,"display_width")
## [1] 9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.227e+02 6.681e+01 -1.837 0.3173
## FTP 4.011e-02 1.608e-02 2.495 0.2427
## UNEMP 3.980e-01 3.643e-01 1.093 0.4719
## M -6.186e-02 2.472e-02 -2.502 0.2420
## LIC 1.916e-02 1.856e-03 10.326 0.0615 .
## GR -2.063e-03 1.494e-03 -1.381 0.3989
## CLEAR -6.854e-02 8.864e-02 -0.773 0.5810
## W 1.191e-04 7.825e-05 1.522 0.3700
## NMAN 8.545e-02 4.327e-02 1.975 0.2984
## G 1.495e-01 7.022e-02 2.129 0.2796
## HE -5.890e+00 2.571e+00 -2.291 0.2620
## WE 2.829e-01 6.329e-02 4.470 0.1401
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4042 on 1 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9994
## F-statistic: 1792 on 11 and 1 DF, p-value: 0.01842cor(crimenes_detroit[-1])## FTP UNEMP M LIC GR CLEAR
## FTP 1.0000000 0.292549420 0.4182908 0.5685212 0.70236680 -0.9742599
## UNEMP 0.2925494 1.000000000 -0.6517190 -0.1671449 -0.03287467 -0.3056631
## M 0.4182908 -0.651718959 1.0000000 0.6983145 0.62778120 -0.4292949
## LIC 0.5685212 -0.167144946 0.6983145 1.0000000 0.90379302 -0.5548320
## GR 0.7023668 -0.032874670 0.6277812 0.9037930 1.00000000 -0.6831611
## CLEAR -0.9742599 -0.305663099 -0.4292949 -0.5548320 -0.68316112 1.0000000
## W -0.8799107 0.090013214 -0.7652255 -0.7815460 -0.84782368 0.8821079
## NMAN 0.8818185 -0.038947279 0.7503422 0.7848606 0.84329238 -0.8915961
## G 0.8793086 0.007512467 0.7097715 0.8037000 0.86242889 -0.8934356
## HE 0.9365055 0.230809707 0.4536647 0.4216629 0.57352189 -0.9574269
## WE 0.9222516 0.131281687 0.5023350 0.3910047 0.55589549 -0.9362843
## H 0.9640610 0.210141728 0.5464227 0.7262896 0.81628728 -0.9684603
## W NMAN G HE WE H
## FTP -0.87991074 0.88181847 0.879308571 0.9365055 0.9222516 0.9640610
## UNEMP 0.09001321 -0.03894728 0.007512467 0.2308097 0.1312817 0.2101417
## M -0.76522553 0.75034215 0.709771506 0.4536647 0.5023350 0.5464227
## LIC -0.78154603 0.78486061 0.803700034 0.4216629 0.3910047 0.7262896
## GR -0.84782368 0.84329238 0.862428887 0.5735219 0.5558955 0.8162873
## CLEAR 0.88210786 -0.89159606 -0.893435576 -0.9574269 -0.9362843 -0.9684603
## W 1.00000000 -0.99537633 -0.985670491 -0.8622352 -0.8585848 -0.9469213
## NMAN -0.99537633 1.00000000 0.990035941 0.8695054 0.8555694 0.9559347
## G -0.98567049 0.99003594 1.000000000 0.8572888 0.8264266 0.9580545
## HE -0.86223522 0.86950540 0.857288766 1.0000000 0.9827635 0.9134063
## WE -0.85858484 0.85556944 0.826426612 0.9827635 1.0000000 0.8881526
## H -0.94692129 0.95593473 0.958054543 0.9134063 0.8881526 1.0000000III. Regresión Logística
Ahora bien, en ciertos casos, puede ocurrir que uno tenga como variable dependiente una variable no cuantitativa con otro nivel de medición, por ejemplo dicotómica.
Este, por ejemplo, es un caso muy usual en medicina, donde quizá uno quiera saber la probabilidad de supervivencia de un paciente o la efectividad de un tratamiento en términos de variables predictoras14.
Estos son problemas de predicción y se selección de predictores. Pero ¿cómo aplicar lo visto para Regresión Lineal si mi variable dependiente es dicotómica?
Veremos ahora las diferencias que aparecen y las decisiones que deberíanmos tomar para realizar el ajuste de nuestro modelo. Vamos con un ejemplo clásico del libro de Hosmer y Lemeshow15:
Se quiere determinar si la edad está vinculada a la presencia de enfermedades coronarias significativas (CHD)
Se cuenta con una muestra de 100 pacientes.
Ejemplo 4
En el ejemplo del libro, donde lo que interesa analizar es si existe una vinculación o asociación entre la edad y la posibilidad de padecer problemas coronarios. Para esto se cuenta con una tabla de 100 observaciones con las siguientes variables:
ID: Identificador de la persona
AGE: Edad de la persona
ARGP: Variable de agrupamiento de la edad (grupos etarios), con un total de 8 agrupamientos.
CHD: Variable dependiente. Cuando CHD=1 implica presencia de afecciones coronarias y CHD=0, ausencia.
Levanto el archivo con estos 100 casos:
urlfile_3 = 'https://github.com/oblitterator/estadistica_II_UNTREF/raw/main/chd_age.csv'
df_chd = read_csv(url(urlfile_3))## Rows: 100 Columns: 4## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (4): ID, AGE, AGRP, CHD##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.df_chd## # A tibble: 100 × 4
## ID AGE AGRP CHD
## <dbl> <dbl> <dbl> <dbl>
## 1 1 20 1 0
## 2 2 23 1 0
## 3 3 24 1 0
## 4 4 25 1 0
## 5 5 25 1 1
## 6 6 26 1 0
## 7 7 26 1 0
## 8 8 28 1 0
## 9 9 28 1 0
## 10 10 29 1 0
## # … with 90 more rowsLo primero que hago es realizar un gráfico de dispersión para mi variable dependiente (CHD).
plot(df_chd$CHD)
Como se ve, a diferencia de lo que se apreciaba con variables cuantitativas, aquí nos aparecen dos líneas paralelas (una donde y=1 y otra donde y=0). Esta gráfica está lejos de mostrarnos a nosotros una línea de tendencia lineal o de cualquier tipo. Es decir, parecería que la linealidad aquí no puede ser aplicada.
Una opción, entonces, es agrupar los datos. En este caso podemos utilizar los grupos de edad que nos arrojaban los datos crudos y conformar una nueva tabla con esos intervalos:
#Armo un nuevo df pero con las frecuencias agrupadas
#Tabla de frecuencias de grupos etarios por prevalencia de enfermedades coronarias
df_agg_chd = data.frame(count(df_chd$AGRP),
aggregate(CHD ~ AGRP, data=df_chd,mean)[2],
aggregate(CHD ~ AGRP, data=df_chd, sum)[2])
colnames(df_agg_chd) = c('AGRP', 'n', 'mean_CHD', 'CHD_1')
df_agg_chd$CHD_0 = df_agg_chd$n-df_agg_chd$CHD_1Además del agrupamiento, he calculado aquí la media por grupo etario de presencia de enfermedades coronarias. Esta media puede pensarse, también, como la probabilidad (o no) de tener un problema cardíaco dado que uno pertenece a un determinado grupo etario.
Mirando a simple vista las proporciones, parecería ver que hay un crecimiento de la probabilidad de tener un evento cardíaco y la edad.
Veamos que pasa si hago un gráfico ahora con las probabilidades de tener eventos coronarios, sumandole a mi tabla original los valores calculados previamente.
df_chd = merge(df_chd, df_agg_chd[1:3])
plot(mean_CHD ~ AGE, data = df_chd)
Como puede apreciarse, si bien ya hay una curva más notoria, no parece haber una relación lineal sino más bien una suerte de “S”, entre mi variable predictora (la edad) y mi variable predicha (CHD).
Cuando la variable dependiente es dicotómica la relación no puede hacerse directamente entre la variable dependiente y la variable independiente sino que la relación debe buscarse entre la probabilidad de la variable dependiente y el o las variables predictoras que uno quiera considerar.
El modelo que plantee la relación entre la probabilidad del evento y las variables predictoras no puede ser lineal, dada la forma de la distribución. En probabilidades esta forma de distribución de una variable dicotómica es bastante usual; lo que necesitamos ahora es una función que nos permita modelar este comportamiento.
Modelo de Regresión Logística
Modelos logit y probit
Buscar alguna forma de predecir el comportamiento de la probabilidad del evento, en lugar de predecir los valores de la variable dicotómica, usando alguna función de probabilidad adecuada. Los dos más usuales son:
• Modelos logit: utilizan la función logística
• Modelos probit: utilizan la inversa de la distribución normal
Es decir, se diferencian en la función que utilizan para plantear la relación entre la probabilidad del evento y las varibles predictoras. En este curso se utilizará la primera opción.
Los modelos de regresión logística buscan también plantear alguna forma de relación lineal entre los predictores, no con la probabilidad sino justamente con una función aplicada a la probabilidad que se denomina logit.
Si denominamos 𝛑(x) a la probabilidad de que la variable dependiente valga 1 (o sea, que se produzca el evento) dado un determinado valor del predictor(y), podemos modelar dicha relación para un sólo predictor con la función logística, quedándonos la siguiente expresión:
Donde, como en otros casos: 𝛃0 (ó ⍺ ) = Intercepto; 𝛃1 = Pendiente. Vemos ahora que, sin que importe qué valores tome x siempre vamos a predecir valores dentro del rango \[0,1\].
Esto puede “linealizarse” aplicando algunas operaciones matemáticas:
¿Pero qué significa esta función y sus términos?
Bien, en primer lugar tenemos que descomponer los términos de la ecuación, para entender qué es lo que estamos “linealizando”.
Odds-Ratio (o Razón de Probabilidades u Oportunidad)16
En principio, podemos decir que el término que encontramos del lado izquierdo de nuestra ecuación es el odds-ratio:
La cantidad 𝛑(x)/1-𝛑(x) se denomina odds-ratio y lo que expresa es la relación entre la probabilidad de que y=1 (𝛑(x)) y la probabilidad de que y=0 (1-𝛑(X)).
El odds-ratio toma valores entre 0 e infinito.
Si hacemos el logaritmo del mismo vemos que presenta una relación lineal con x:
Intentaremos explicar qué mide el odds-ratio a través de un ejemplo.
Ejemplo 5
En medicina, el odds-ratio está muy generalizado para interpretar el riesgo de ocurrencia o el efecto de una determinada variable de una enfermedad sobre un paciente. En este caso, se evaluó la asociación entre obesidad e hipertensión arterial en escolares, empleando un estudio transversal.
La muestra estuvo constituida por 2.980 escolares entre 6 y 14 años, de los cuales 622 eran obesos (162 hipertensos y 460 no hipertensos) y 1.593 eran eutróficos (142 hipertensos y 1.451 no hipertensos). El resto de los escolares (n = 765) presentaban otros categorías de diagnóstico nutricional (bajo peso y sobrepeso) y no fueron considerados.
obesos = c(162,460)
no_obesos = c(142,1451)
hta_obesidad = data.frame(obesos, no_obesos)
rownames(hta_obesidad) = c("presente", "ausente")
hta_obesidad## obesos no_obesos
## presente 162 142
## ausente 460 1451Viendo sólo los números absolutos de hipertensos, pareciera ser que la cantidad es la misma para los dos estados nutricionales (obesos = 162, no_obesos = 142).
Pero rara vez es recomendable las frecuencias absolutas para hacer comparaciones, por eso es que tenemos que definir indicadores que nos permitan determinar si la proporción de hipertensos es más alta entre los obesos que entre los no obesos. Esto puede hacerse de dos maneras distintas:
El Riesgo o probabilidad se define como la cantidad de hipertensos en este caso. Sería básicamente la ocurrencia de nuestro evento sobre los casos. Para este ejemplo:
hta_obesidad[1]/sum(hta_obesidad$obesos)## obesos
## presente 0.2604502
## ausente 0.7395498Y ahí decís: listo, no pasa naranja, re baja prevalencia. Pero después chequeas entre los no obesos y:
hta_obesidad[2]/sum(hta_obesidad$no_obesos)## no_obesos
## presente 0.08913999
## ausente 0.91086001Apa. Qué pasa acá. Entre el 91% de los no obesos no hay enfermedades coronarias ¿Algo estará sucediendo estadísticamente?
#DF con los riesgos sobre cada uno de sus denominadores
hta_riesgo = data.frame(c(hta_obesidad[1]/sum(hta_obesidad$obesos)),
c(hta_obesidad[2]/sum(hta_obesidad$no_obesos))
)
hta_riesgo ## obesos no_obesos
## 1 0.2604502 0.08913999
## 2 0.7395498 0.91086001Pues bien en el riesgo relativo uno analiza la probabilidad de tener hipertensión siendo obeso, sobre tener hipertensión no siendo obeso. Eso, da -y ya dije muchas veces eso- un riesgo relativo de:
#Resultado de casos del evento presentes en obesos sobre casos presentes en no obesos:
hta_riesgo$obesos[1]/hta_riesgo$no_obesos[1]## [1] 2.921811Esto dice que si uno es obeso la probabilidad, o riesgo relativo, de que uno sea hipertenso es 2.92 veces más alta que si uno no fuera obeso.
Si calculamos la oportunidad (odds-ratio) de ser hipertenso siendo obeso sería:
#Hago el Odds Ratio ahora
(hta_riesgo$obesos[1]/hta_riesgo$obesos[2])/
(hta_riesgo$no_obesos[1]/hta_riesgo$no_obesos[2])## [1] 3.598622Es decir, que la “oportunidad” de tener HTA siendo obeso se triplica, o es 3,6 veces más alta, que siendo no obeso.
El Odds-Ratio tiene la capacidad entonces de cuantificar el efecto de una variable independiente (estado nutricional) sobre una variable dependiente (HTA), como se ve en el ejemplo anterior.
El Odds es en sí la probabilidad de la ocurrencia del evento. El logaritmo del odds-ratio es lo que suponemos relacionado linealmente con el predictor (aclaración: este caso para predecir una “x” dicotómica”).
La importancia del los Odds Ratio reside en que los mismos permiten interpretar los resultados del Modelo de Regresión Logística, ¿cómo? Pues bien:
El Odds es un cociente entre dos situaciones: la del numerador que representa la influencia del factor que uno quiere analizar y la del denominador es la situación de ausencia de ese factor. Entonces:
- Un Odds-Ratio de 1 estaría mostrando un valor “neutro” entre las dos variables (estaría mostrando una no relación o no efecto entre las variables)
- Un Odds-Ratio mayor a 1 significa que ese factor incrementa el riesgo. En presencia de obesidad uno tiene más chance de ser hipertenso, en el ejemplo anterior.
- Por el contrario, un Odds-Ratio menor a 1 implica que la presencia de ese factor lo que hace es disminuir el riesgo.
Esta relación lineal que surge del Odds es la base teórica del Modelo de Regresión Logística, pero obviamente su potencialidad parte de poder relacionar varios predictores. Si sólo se tratase de dos variables con el Odds-Ratio sería suficiente.
Retomando la Regresión Logística
Veamos ahora un ejemplo con distintas variables predictoras, pero primero tomando la relación entre dos dicotómicas para ejemplificar la lectura del Odds.
Ejemplo 6: Bajo peso al nacer
El bajo peso al nacer, definido por un peso al nacer inferior a 2500 gr., ha sido una preocupación de los médicos durante años debido a que tanto las tasas de mortalidad como la de nacimientos defectuosos son muy altas para los niños con bajo peso al nacer.
El comportamiento de la mujer durante el embarazo (incluyendo la dieta, los hábitos tabáquicos y los cuidados prenatales) pueden alterar las chances de un parto de un niño con bajo peso.
Las variables con las que contamos en la tabla son:
• 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 (ó “x”, la variable a predecir)
Cargamos el DF:
url_4 = 'https://raw.githubusercontent.com/oblitterator/estadistica_II_UNTREF/main/lowbwt.csv'
lowbwt = read.csv(url(url_4))
head(lowbwt)## ID LOW AGE LWT RACE SMOKE PTL HT UI FTV BWT
## 1 4 1 28 120 3 1 1 0 1 0 709
## 2 10 1 29 130 1 0 0 0 1 2 1021
## 3 11 1 34 187 2 1 0 1 0 0 1135
## 4 13 1 25 105 3 0 1 1 0 0 1330
## 5 15 1 25 85 3 0 0 0 1 0 1474
## 6 16 1 27 150 3 0 0 0 0 0 1588En primer lugar, vamos a calcular el Odds sólo para una varible predictora: SMOKE (si la madre fumó durante el embarazo).
odds_bajo_peso = lowbwt %>%
mutate(bajo_peso = ifelse(BWT <=2500, 1, 0)) %>%
dplyr::select(SMOKE, bajo_peso)
odds_bajo_peso = as.data.frame.matrix(with(odds_bajo_peso,
table(bajo_peso,SMOKE)),
row.names = c("no_bp","si_bp"))
odds_bajo_peso = odds_bajo_peso %>%
setNames(c("no_fumo", "si_fumo"))
odds_bajo_peso## no_fumo si_fumo
## no_bp 86 44
## si_bp 29 30(odds_bajo_peso$no_fumo[1]/odds_bajo_peso$no_fumo[2])/ #Ahora calculanos el Odds
(odds_bajo_peso$si_fumo[1]/odds_bajo_peso$si_fumo[2])## [1] 2.021944Como puede apreciarse, tomando sólo la variable dicotómica “SMOKE” podemos decir que para las madres fumadoras se incrementa en 2.02 (“duplica”) la chance de tener hijos con bajo peso.
Esto también podría expresarse (si linealizamos el Odds en base al logaritmo), de la siguiente manera:
Este es el ejemplo más simple de una regresión logísitca en base a una variable dicotómica (SMOKE, en el ejemplo). Como se ve en la ecuación lineal, un parámetro 𝛃1 distinto de 0 nos indica en este casi que SMOKE es un predictor de bajo peso (si 𝛃1 no cumpliese el supuesto de ser distinto de 0, no sería estadísticamente significativo cosa que se analizará más adelante con R)
En el modelo de regresión lineal, los 𝜷i eran el cambio promedio en y ante un cambio unitario en x. En la regresión logística, incrementar una unidad en x, cambia el logaritmo del odds-ratio en 𝜷i. O, lo que es lo mismo, multiplica el odds por e𝜷i La relación entre 𝛑(x) y x no es una línea recta: cuánto cambia 𝛑(x) ante un cambio unitario en x depende de los valores de X. Aún así, el signo de 𝜷i expresa la dirección de cambio en 𝛑(x) (independientemente del valor de x).
Reg. Logística para variable dicotómica
Bien, ahora que conocemos los fundamentos teóricos y formuleicos del asunto, vamos a utilizar la función glm() para construir nuestro modelo de regresión logística, y analizar los coeficientes, antes de avanzar en la inclusión de más variables predictoras:
logi_fumar = glm(LOW~SMOKE,data=lowbwt,family="binomial") #Se aclara que la familia es "binomial" ya que estamos trabajando con varaibles dicotómicas. Para otros casos consultar otros parámetros.
summary(logi_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: 4Como se ve aquí arriba -además de arrojarnos el intercepto- la función glm() nos indica que SMOKE es estadísticamente significativo como predictor del bajo peso al nacer (p_value = 0,0276)
OddsRatio(logi_fumar) #Con esta función, obtenemos los Odds Ratio, así como el intervalo de confianza de los mismos## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ SMOKE, family = "binomial", data = lowbwt)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.337 0.218 0.507 4.14e-07 ***
## SMOKE 2.022 1.082 3.801 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.209 Nagelkerke R2: 0.036#####
#Tambien puede calcularse el OddsRatio de la siguiente forma:
# library(oddsratio)
# or_glm(data = LOWBWT, model = fumar, incr = list(SMOKE=1))Como podemos observar en la tabla precedente, el Odds de Fumar es prácticamente el mismo que calculado a mano previamente. Lo interesante es que esta función aplicada al modelo nos devuelve el intervalo de confianza del Odds: entre 1,082 en su límite inferior y 3,801 para su limite superior. Siendo nuestra H0= 1, hecho de que el valor “1” no esté contemplado dentro de nuestro intervalo (que, como se dijo antes, indicaría un efecto neutral de la variable) es lo que nos permite asumir la significatividad de este parámetro en nuestro modelo.
Reg. Logística para variable cualitativa no dicotómica
Para aquellos casos donde tomamos para la construcción de un modelo de regresión logística variables cualitativas no dicotómicas (es decir, más de dos posibles valores no ordinales) es necesario convertir en “dummies” la misma. Esta es una forma de dicotomizar nuestra variable, aumentando la cantidad de features de nuestro data frame.
Supongamos que, para predecir el bajo peso al nacer, queremos utilizar como predictor la variable “RACE” (no pienso ni hacer un comentario sobre la fijación de los que están del hemisferio norte con este tema), que tiene tres posibles valores no ordinales:
1 = White
2 = Black
3 = Other
(““No ordinales”“)
Pues bien, lo que tenemos que hacer es aplicar una función para obtener los dummies en base a dicha columna y pasarsela a nuestro DF. Es importante recalcar que, para no producir multicolinealidad, siempre la cantidad de columnas dummies debe ser el total de variables menos 1.
#install.packages('fastDummies') #instalo el paquete
#library(fastDummies) # lo llamo
lowbwt_race = lowbwt %>%
dplyr::select(LOW, RACE) %>%
fastDummies::dummy_cols('RACE') #En realidad, debería agregar la sintaxis ",remove_first_dummy = TRUE", ya que, como se dijo antes, se busca evitar la multicolinealidad. En este caso voy a dejar el primer para armar el modelo y ver los coeficientes para cada combinatoria de dummies.Y ahora, conformo un primer modelo como antes pero tomando mi variable LOW y mis dos dummies generadas (RACE_1 y RACE_2):
logi_raza_1 = glm(LOW~RACE_1+RACE_2,data=lowbwt_race,family="binomial") #Se aclara que la familia es "binomial" ya que estamos trabajando con varaibles dicotómicas. Para otros casos consultar otros parámetros.
summary(logi_raza_1)##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_2, family = "binomial", data = lowbwt_race)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0489 -0.9665 -0.7401 1.4042 1.6905
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5188 0.2526 -2.054 0.0400 *
## RACE_1 -0.6362 0.3478 -1.829 0.0674 .
## RACE_2 0.2086 0.4705 0.443 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4OddsRatio(logi_raza_1)## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_2, family = "binomial", data = lowbwt_race)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.595 0.358 0.969 0.0400 *
## RACE_1 0.529 0.266 1.045 0.0674 .
## RACE_2 1.232 0.483 3.093 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.209 Nagelkerke R2: 0.037Con LOW y RACE_1 y RACE_3
logi_raza_2 = glm(LOW~RACE_1+RACE_3,data=lowbwt_race,family="binomial")
summary(logi_raza_2)##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_3, family = "binomial", data = lowbwt_race)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0489 -0.9665 -0.7401 1.4042 1.6905
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3102 0.3970 -0.781 0.4346
## RACE_1 -0.8448 0.4634 -1.823 0.0683 .
## RACE_3 -0.2086 0.4705 -0.443 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4OddsRatio(logi_raza_2)## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_3, family = "binomial", data = lowbwt_race)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.733 0.328 1.587 0.4346
## RACE_1 0.430 0.173 1.080 0.0683 .
## RACE_3 0.812 0.323 2.072 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.209 Nagelkerke R2: 0.037Y finalmente, con RACE_2 y RACE_3
logi_raza_3 = glm(LOW~RACE_2+RACE_3,data=lowbwt_race,family="binomial")
summary(logi_raza_3)##
## Call:
## glm(formula = LOW ~ RACE_2 + RACE_3, family = "binomial", data = lowbwt_race)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0489 -0.9665 -0.7401 1.4042 1.6905
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1550 0.2391 -4.830 1.36e-06 ***
## RACE_2 0.8448 0.4634 1.823 0.0683 .
## RACE_3 0.6362 0.3478 1.829 0.0674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4OddsRatio(logi_raza_3)## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ RACE_2 + RACE_3, family = "binomial", data = lowbwt_race)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.315 0.193 0.495 1.36e-06 ***
## RACE_2 2.328 0.926 5.775 0.0683 .
## RACE_3 1.889 0.957 3.758 0.0674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.209 Nagelkerke R2: 0.037Como vemos, cualquier combinatoria de dummies de RACE no es estadísticamente significativo para predecir el bajo peso al nacer. Chequeemos los Odds Ratio de alguno de ellos (es indistinto):
OddsRatio(logi_raza_1)## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_2, family = "binomial", data = lowbwt_race)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.595 0.358 0.969 0.0400 *
## RACE_1 0.529 0.266 1.045 0.0674 .
## RACE_2 1.232 0.483 3.093 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.209 Nagelkerke R2: 0.037En ambos casos, para un intervalo de confianza del 95% no rechazamos la Hipótesis Nula, ya que entre el límite inferior y superior del Odds Ratio está incluido el valor 1, la neutralidad del efecto.
En resumen: lo que hacen las dummies es simular el comportamiento de las variables cualitativas no dicotómicas, para poder utilizar la misma ecuación de cálculo vista anteriormente. Como vemos, cualquier combinación de las mismas
Jackpot: en realidad no es necesario hacer todo ese preprocesamiento con la función de dummies que hice antes. Simplemente tenemos que pasarle a glm() las variables dentro del modelo como factor(x):
glm(LOW~factor(RACE),data=lowbwt_race,family="binomial")##
## Call: glm(formula = LOW ~ factor(RACE), family = "binomial", data = lowbwt_race)
##
## Coefficients:
## (Intercept) factor(RACE)2 factor(RACE)3
## -1.1550 0.8448 0.6362
##
## Degrees of Freedom: 188 Total (i.e. Null); 186 Residual
## Null Deviance: 234.7
## Residual Deviance: 229.7 AIC: 235.7summary(logi_raza_2)##
## Call:
## glm(formula = LOW ~ RACE_1 + RACE_3, family = "binomial", data = lowbwt_race)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0489 -0.9665 -0.7401 1.4042 1.6905
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3102 0.3970 -0.781 0.4346
## RACE_1 -0.8448 0.4634 -1.823 0.0683 .
## RACE_3 -0.2086 0.4705 -0.443 0.6575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4Como vemos, ya nos arma las dummies el propio modelo (RACE)2 y (RACE)3. Es indistinto que categoría se elimina para hacer la dicotomización, ya que igual se encuentra expresada por la combinación de las otras variable, pero si es importante después ver el efecto (la chance o el riesgo de ocurriencia) entre los Odds de las variables predictoras y las predichas.
Lo que si podemos medir es el tamaño del efecto, en caso de que las mismas sean significativas, sobre la variable a predecir. Por ejemplo, cuanto ser de “Raza Negra” (INADI Alert) aumenta la chance de que el niño tenga bajo peso al nacer.
Reg. Logística para variable cuantitativa
Para aquellos casos donde tomamos para la construcción de un modelo de regresión logística variables cuantitativas, es importante primer comprobar que las mismas tengan
a) Relación lineal con la distribución logística y
b) Estandarizarlas para que la unidad de medida de las mismas no influya en sus pesos dentro del modelo. Esto para casos de una sola variable en realidad no es necesario estrictamente, pero sí cuando tenemos dentro del mismo dos o más.
Entonces, si queremos trabajar con las variables cuantitativas continuas AGE (edad de la madre) y LWT (peso de la madre), primero comprobamos su correlación con la distribución.
library(ltm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:olsrr':
##
## cement## The following objects are masked from 'package:raster':
##
## area, select## The following object is masked from 'package:rstatix':
##
## select## The following object is masked from 'package:dplyr':
##
## select## Loading required package: msm## Loading required package: polycorbiserial.cor(lowbwt$AGE, lowbwt$LOW, use = c("all.obs", "complete.obs"))## [1] 0.1189393biserial.cor(lowbwt$LWT, lowbwt$LOW, use = c("all.obs", "complete.obs"))## [1] 0.1696269Parecería ser que la relación entre nuestras variable cuantitativas predictoras y la predicha es baja, aunque mayor entre la chance de nacer con bajo peso y el peso de la madre. Ahora, calculamos los coeficientes utilizando las dos variables:
sum(is.na(lowbwt$AGE)) #chequeo si hay valores nulos## [1] 0sum(is.na(lowbwt$LWT)) ## [1] 0lowbwt = lowbwt %>% mutate_each_(funs(scale(.) %>% as.vector),
vars=c("AGE","LWT")) #esta hermosa pieza de código lo que me permite es elegir con scale() sólo las variables a escalar, en vez de todas.## 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))logi_edad_peso_madre = glm(LOW ~ AGE+LWT,data=lowbwt,family="binomial") #Se aclara que la familia es "binomial" ya que estamos trabajando con varaibles dicotómicas. Para otros casos consultar otros parámetros.
summary(logi_edad_peso_madre)##
## Call:
## glm(formula = LOW ~ AGE + LWT, family = "binomial", data = lowbwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1352 -0.9088 -0.7480 1.3392 2.0595
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8343 0.1636 -5.099 3.42e-07 ***
## AGE -0.2108 0.1711 -1.232 0.2178
## LWT -0.3907 0.1899 -2.057 0.0397 *
## ---
## 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: 227.12 on 186 degrees of freedom
## AIC: 233.12
##
## Number of Fisher Scoring iterations: 4Bueno, si bien antes nos había dado una correlación que parecía “baja” entre LWT y LOW (r=0,17), el summary de nuestro modelo nos estaría mostrando una cierta significatividad de la variable en nuestro modelo. Veamos que pasa si es ruido producido por la variable AGE (claramente, con un p value bajo) o si realmente es así:
logi_peso_madre = glm(LOW ~ LWT,data=lowbwt,family="binomial") #Se aclara que la familia es "binomial" ya que estamos trabajando con varaibles dicotómicas. Para otros casos consultar otros parámetros.
summary(logi_peso_madre)##
## Call:
## glm(formula = LOW ~ LWT, family = "binomial", data = lowbwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0951 -0.9022 -0.8018 1.3609 1.9821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8267 0.1625 -5.087 3.64e-07 ***
## LWT -0.4299 0.1887 -2.279 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 228.69 on 187 degrees of freedom
## AIC: 232.69
##
## Number of Fisher Scoring iterations: 4Bueno, aumenta el nivel de significatividad del parámetro del estimador, por lo que podemos decir que LWT tiene cierto “efecto” sobre el bajo peso al nacer. Finalmente, calculemos el Odds:
OddsRatio(logi_peso_madre)## Waiting for profiling to be done...##
## Call:
## glm(formula = LOW ~ LWT, family = "binomial", data = lowbwt)
##
## Odds Ratios:
## or or.lci or.uci Pr(>|z|)
## (Intercept) 0.438 0.315 0.597 3.64e-07 ***
## LWT 0.651 0.438 0.922 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Brier Score: 0.208 Nagelkerke R2: 0.044El Odds Ratio entre una variable dicotómica y una cuantitativa lo que indica es un aumento o disminución de la probabilidad (sentido, no tamaño), por lo que menos de uno indica disminución y valores mayores a uno, aumento.
Lo que nos estaría diciendo el Odds para un intervalo de confianza del 95% es que a medida que disminuye el peso de la madre disminuye la chance de bajo peso. Y acá es donde hay que quizá pensar un poco más allá de los números, creo.
Habría que ver quizá si esto aumentando un poco el tamaño de la muestra también resulta estadísticamente significativo, ya que los estudios médicos indican que a mayor peso de la madre menos chance de bajo peso.
En fin, la
Selección de mejor modelo logístico de predicción
Así como en el caso de la Regresión Lineal Múltiple el método de selección del mejor modelo es a través del Método de Mínimos Cuadrados (encontrar la recta a partir de la minimización de los residuos entre el valor real y el predicho), en el caso de la Logística se utiliza el Método de Máxima Verosimilitud.17
Teóricamente, como estamos trabajando con chances, nuestra búsqueda está puesta en maximizar la probabilidad conjunta comparando las distintas combinaciones entre las variables (o features) de una muestra obtenida.
Entonces, busquemos los parámetros de la recta que haga que la probabilidad conjunta de esa combinación de variables de nuestra muestra sea máxima (de alguna forma, siempre lo que se hace en estadística es encontrar los extremos y medir lo que hay entre ellos).
Para medir esto, lo que se hace es calcular la deviancia de las funciones de máxima verosimilitud en cada modelo y compararlas entre sí.
Entonces, esto lo que nos permite es automatizar la selección de nuestros modelos midiendo las deviancias (que siguen una forma de distribución del tipo chi cuadrado) entre los modelos. Con eso podemos ver si la relación de cada combinatoria es estadísticiamente significativa para predecir nuestra variable dicotómica.
Por ejemplo: tenemos un modelo sólo con nuestra variable SMOKE y otro modelo con SMOKE y RACE (es decir, el segundo incluiría al primero). Lo que realizamos es calcular y medir la relación entre las deviancias de cada modelo. Estas, están distribuidas como una función de chi cuadrado, con k-1 grados de libertad, que corresponden al total de variables (𝛃k) de cada modelo menos uno.
Es decir, utilizamos la función de probabilidad Chi Cuadrado para ver cuán estadísticamente significativo es el valor de deviancia de cada modelo. Esta es la base para el método de selección automática de los predictores para los modelos de regresión logística.
Capacidad de clasificación
Es posible clasificar a los individuos en función de los estimadores del modelo
Como es posible calcular la probabilidad de cada sujeto, se puede “predecir” un valor de Y en términos de esa probabilidad, y calcular los aciertos
Para esto, es necesario decidir a partir de qué probabilidad se considera al individuo como que tuvo el evento.
Entonces, volviendo, nosotros tenemos una variable dependiente que es dicotómica. Ajustamos un modelo (estimando los parámetros de 𝛃0, 𝛃1… 𝛃k) que permite predecir la probabilidad de que la variable dependiente sea 1 para cada individuo (x).
Así, al poder predecir la probabilidad de que un individuo valga 1, yo puedo también predecir cuando vale 0, a partir de fijar un valor dentro de esa probabilidad que crea conveniente para cada caso. Es decir, un punto o probabilidad de corte.
La gran mayoría de los software toman como punto de corte 0,5 (50%): si la probabilidad es mayor a 0,5 se clasifica el individuo (x) como 1 (Y=1); caso contrario, como 0 (Y=0). Este punto de corte puede ser optimizado (como se verá más adelante) pero lo importante es comprender que en base a este valor que predice nuestro modelo y el valor observado de la dependiente podemos armar una tabla de contingencia conocida como “Matriz de Confusión”.
Matriz de confusión
“En el campo de la inteligencia artificial y en especial en el problema de la clasificación estadística, una matriz de confusión es una herramienta que permite la visualización del desempeño de un algoritmo que se emplea en aprendizaje supervisado. Cada columna de la matriz representa el número de predicciones de cada clase, mientras que cada fila representa a las instancias en la clase real. Uno de los beneficios de las matrices de confusión es que facilitan ver si el sistema está confundiendo dos clases.”18
Una matriz de confusión puede expresarse, entonces, de la siguiente manera:
| Éxito (predicho) | Fracaso (predicho) | ||
|---|---|---|---|
| Éxito (observado) | A | B | A+B |
| Fracaso (observado) | C | D | C+D |
| A+C | B+D | Total (A+B+C+D) |
Entonces, la frecuencias de A y D nos indica los casos en los cuales nuestro modelo predijo Éxito/Fracaso y dicha predicción se corresponde con el valor observado (son frecuencias absolutas). Asimismo, un modelo de regresión logística con buen ajuste, tendería a valores bajos de C y B. En otras palabras, los valores de C y B son los valores de “mala” clasificación, es decir, donde el modelo predijo algo distinto a lo que se observo.
La mensura de estas dos situaciones reciben el nombre de Sensitividad y Especificidad:
Sensitividad: falsos positivos (suma de casos positivos que mi modelo predijo como ocurrencia del evento según su chance, sobre la frecuencia absoluta de los positivos en la muestra)
Especificidad: falsos negativos (suma de casos negativos que mi modelo predijo como no ocurrencia del evento, sobre la frecuencia absoluta de la no ocurrencia del suceso en la muestra)
Es importante calcular ambas métricas, ya que se espera que ambos valores sean altos, y por ende el error de clasificación del modelo sea bajo.
La curva ROC
La curva ROC surge de visualizar el área entre el valor de la Sensibilidad (los verdaderos positivos) y el valor de 1-Especificidad (los falsos positivos). La situación ideal sería Sensitividad Máxima (igual a 1) y 1-Especificidad sea 0 (o sea, que la Especificidad sea 1).
Es decir, que sea un rectángulo de 1 de área (solo por motivos de visualización se colocan el valor de la proporción de la Especificidad invertida en el eje de las x).
Obviamente, un caso de 1 para ambos valores es imposible (es decir, un valor del área de la curva ROC = 1), y si sucede es necesario evaluar si se estás overfitieando como un campéon. Por otro lado, para tener un valor de referencia una curva de 0,5 estaría indicando que nuestro modelo no distingue mejor que el azar.
Pero no sólo nos permite esto, sino que nos permite definir nuestro punto de corte. Es decir, la probabilidad de predicción que tiene mayor balance entre Sensibilidad y Especificidad, entre verdaderos positivos y falsos positivos. Obviamente (?) también podemos privilegiar una de las dos situaciones por motivos empíricos. Nos preocupa más que un test se quede más largo que corto. O viceversa.
Todo esto en general igual está más enfocado a ámbitos de ciencias naturales o medicina que a ciencias sociales; no obstante, es muy importante a la hora de trabajar en la realización de Modelos de Aprendizaje Supervisado y no Supervisado.
Selección automática de predictores en R
Bueno, para finalizar veamos -así como se analizó en los modelos de regresión lineal- de que manera realizar una selección automática de predictores con R. Estos estarán basados en los valores antes analizados de la deviancia, el AIC y sus variantes.
En primer lugar, se corre un modelo de regresión logística sin predictores, asignándolo a una variable para guardar los resultados.
birthwt.glm <- glm(LOW ~ 1, # al colocar ~ 1 instanciamos el modelo primero sin los predictores
family = binomial,
data = lowbwt) En segunda instancia, con la función stepAIC(), lo que hacemos es determinar en primer lugar nuestro modelo de “mínima” (que es aquel instanciado previamente, sin ningún predictor); en segundo lugar pasamos nuestro modelo de “máxima” (toda la sintaxis que viene después de “upper =”), con todos las predictoras para las que queremos que analice las métricas antes mencionadas.
birthwt.step = stepAIC(birthwt.glm, 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.46Extraído de: https://www.cienciadedatos.net/documentos/19_anova.html↩︎
Esta es la misma fórmula, a grandes rasgos, utilizada para el cálculo del Coeficiente Razón de Correlación entre una variable cuantitativa y una cualitativa, recordad que allí vimos que:
- VI (o SSIntra) = Suma de los cuadrados de los desvíos de los valores respecto a la media de cada uno de los grupos o estratos que integra.
- VE (o SSEntre) = Suma de los cuadrados de los desvíos de las medias de cada grupo respecto a la media total.
- VT (o SSTotal) = Suma de los cuadrados de los desvíos de los valores respecto a la media total o general (SSTotal = SSEntre+ SSIntra).
Lo ideal sería no tener que crear una columna aparte pero la joda es si utilizo la función as.factor() después no me deja fitearle el modelo para hacer el análisis. En fin. Seguro hay alguna forma, todavía no la encontré.↩︎
HSD viene de Honestly Significant Differences. Un velo de republicanismo hipócrita recorre este apunte.↩︎
Triola p. 643 (XXXX)↩︎
Pueden probar que el resultado que nos arroja el estadístico F en este caso de ANOVA de dos factores para el factor “Fedad” es distinto que cuando lo hacíamos para un solo factor chequeando el resultado que nos dio el mismo más arriba.↩︎
Esto igual depende de qué tipo de universo y variables estemos analizando. En el caso de las ciencias físicas puede ser requerido un valor mucho mayor, pero en ciencias sociales el r=0,5 es considerado “alto”.↩︎
La explicación de porqué se considera 0 parte en realidad de las leyes de la probabilidad. Si entendemos que el error se debe a causas aleatorias, el mismo se distribuye como una normal a larga escala, por lo que termina anulándose, en la teoría.↩︎
Se hace al cuadrado ya que, como algunos residuos pueden ser negativos y otros positivos, no queremos que se sustraigan entre sí. Por eso antes de sumarlos los elevamos para que tengan el mismo signo antes de sumarlos.↩︎
La PdH para chequear si los parámetros ⍺ y/o 𝛽 son iguales o diferentes a 0 (es decir, si son estadísticamente significativos), tiene por objeto eliminar variables de mi fórmula del modelo de regresión. En caso de que ⍺ (la ordenada al origen) sea cero, directamente hace que no tome esta (o estas) variable en cuenta para el cálculo de mis correlaciones (me ayuda a construir el modelo más parsimonioso).↩︎
Por ejemplo, la predicción de pobreza estructural -nbi1- respecto a las variables dependientes calmat3 y calmat4, o en el caso más claro, nbi2, que aparece en el Ejercicio 1, sobre Córdoba.↩︎
Otras posibilidadades son: combinar las variables altamente correlacionadas, hacer componentes principales y usar como predictores las nuevas variables o utilizar la Regresión de Ridge.↩︎
Las Necesidades Básicas Insatisfechas fueron definidas según la metodología utilizada en “La pobreza en la Argentina” (Serie Estudios INDEC. N° 1, Buenos Aires, 1984).
Los hogares con Necesidades Básicas Insatisfechas (NBI) son los hogares que presentan al menos uno de los siguientes indicadores de privación:
1- Hacinamiento: hogares que tuvieran más de tres personas por cuarto.
2- Vivienda: hogares en una vivienda de tipo inconveniente (pieza de inquilinato, vivienda precaria u otro tipo, lo que excluye casa, departamento y rancho).
3- Condiciones sanitarias: hogares que no tuvieran ningún tipo de retrete.
4- Asistencia escolar: hogares que tuvieran algún niño en edad escolar (6 a 12 años) que no asistiera a la escuela.
5- Capacidad de subsistencia: hogares que tuvieran cuatro o más personas por miembro ocupado y, además, cuyo jefe no haya completado tercer grado de escolaridad primaria.↩︎
En medicina se conocen como “factores de riesgo”, pero también en otras ciencias sociales se busca a veces predecir si una persona puede o no conseguir trabajo según determinadas características personales o de su entorno.↩︎
Hosmer & Lemeshow (2000), Applied Logistic Regression Wiley-. Interscience … McGraw-Hill / Interamericana de España,S.A.U..↩︎
Increíblemente también conocida como Razón de Momios↩︎
Este también podría utilizarse para la Regresión Lineal (no al revés), pero el MMC tiene una visualización más directa de la relación.↩︎