Son todos lo mismo

7 minute read

Es muy probable que si tuvieron cursos de estadística en la facultad estén familiarizados (al menos de nombre) con el test t, la regresión y el análisis de la varianza (ANOVA). En estos cursos usualmente se presentan como técnicas para diferentes situaciones y que no tienen nada en común. Por eso puede resultarles extraño si les digo que en realidad son todos el mismo análisis.

corea

¿Cómo es eso? ¿estás tomando desde temprano, Ger? no, en serio son todos casos particulares de algo que se llama modelo lineal. Pero no necesitás creer en mi palabra lo podemos ver con un ejemplo práctico.

Regresión

En este problema lo que buscamos es predecir cuanto le va a costar la cobertura médica a un ciudadano de Estados Unidos. El dataset se puede descargar desde aquí.

> seguro <- read.csv("insurance.csv")

¿Qué información tenemos? chusmeemos un poco:

> str(seguro)
'data.frame':   557 obs. of  7 variables:
 $ age     : int  19 33 60 62 27 19 23 30 34 59 ...
 $ sex     : Factor w/ 2 levels "female","male": 1 2 1 1 2 2 2 2 1 1 ...
 $ bmi     : num  27.9 22.7 25.8 26.3 42.1 ...
 $ children: int  0 0 0 0 0 1 0 0 1 3 ...
 $ smoker  : Factor w/ 2 levels "no","yes": 2 1 1 2 2 1 1 2 2 1 ...
 $ region  : Factor w/ 4 levels "northeast","northwest",..: 4 2 2 3 3 4 1 4 1 3 ...
 $ charges : num  16885 21984 28923 27809 39612 ...

> head(seguro)
   age    sex    bmi children smoker    region   charges
1   19 female 27.900        0    yes southwest 16884.924
4   33   male 22.705        0     no northwest 21984.471
10  60 female 25.840        0     no northwest 28923.137
12  62 female 26.290        0    yes southeast 27808.725
15  27   male 42.130        0    yes southeast 39611.758
16  19   male 24.600        1     no southwest  1837.237

Contamos con las siguientes variables predictivas:

  • sexo: femenino o masculino
  • bmi: índice de masa corporal (peso / altura^2)
  • children: número de hijos
  • smoker: si fuma
  • region: región de EE.UU. donde vive

En este ejemplo solo utilizaremos el BMI, por cuestiones de simplicidad. Lo primero que nos interesa es hacer una regresión entre charges (costos de la cobertura médica) y el BMI. Nuestra hipótesis es que a mayor BMI aumentarán los costos para el asegurado, ya que se incrementan los riesgos de sufrir enfermedades. Veamos si es así.

> modelo1 <- lm(charges ~ bmi, data=seguro)

El nombre de la función ya nos está diciendo algo: lm es el acrónimo de linear model y es la misma función que utilizaremos para el ANOVA. Para ver los resultados usamos la función summary:

> summary(modelo1)                                                                                                                              
                                                                                                                                                
Call:                                                                                                                                           
lm(formula = charges ~ bmi, data = seguro)                                                                                                      
                                                                                                                                                
Residuals:                                                                                                                                      
   Min     1Q Median     3Q    Max                                                                                                              
-25407  -5916   -582   5784  35621                                                                                                              
                                                                                                                                                
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -31617.21    1628.28  -19.42   <2e-16 ***
bmi           1929.82      57.96   33.29   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8209 on 555 degrees of freedom
Multiple R-squared:  0.6664,    Adjusted R-squared:  0.6658 
F-statistic:  1109 on 1 and 555 DF,  p-value: < 2.2e-16

En este caso vemos que existe una relación lineal entre los gastos en cobertura médica y el índice de masa corporal del individuo (p<0.01).

¿Que pasaría si la variable explicativa fuera categórica en lugar de continua? Siguiendo con el ejemplo de la cobertura médica, discreticemos la variable bmi en intervalos:

  • BMI < 25: normal
  • BMI >= 25 y <= 29: sobrepeso
  • BMI > 29: obeso

En R lo podemos hacer así:

> seguro$bmicat <- cut(seguro$bmi, breaks=c(0,25,29,100),
                      labels=c("normal", "sobrepeso", "obeso"))

¿Cómo analizamos ahora la relación entre el costo y el bmi discretizado? una opción sería hacer una regresión con variables dummy. Las variables dummy son variables (valga la redundancia) que solo aceptan dos valores: 1 (se cumple la condición) o 0 (no se cumple). La cantidad de dummies siempre será igual a la cantidad de categorías - 1, en este caso 2. La forma más sencilla de crearlas en R es con la función model.matrix:

> seguro <- cbind(seguro, with(seguro, model.matrix(~bmicat))[,2:3])
> head(seguro)
  age    sex    bmi smoker   charges     bmicat bmicatobeso bmicatsobrepreso
1  19 female 27.900    yes 16884.924 sobrepreso           0                1
2  33   male 22.705     no 21984.471     normal           0                0
3  60 female 25.840     no 28923.137 sobrepreso           0                1
4  62 female 26.290    yes 27808.725 sobrepreso           0                1
5  27   male 42.130    yes 39611.758      obeso           1                0
6  19   male 24.600     no  1837.237     normal           0                0

¿Cómo codificaríamos a un sujeto normal, a uno con sobrepeso y a uno obeso con dummies? Pues de esta manera:

  bmicatobeso bmicatsobrepreso
normal 0 0
sobrepeso 0 1
obeso 1 0

Ahora podemos proceder a correr la regresión para esta variable dummy:

> modelo2 <- lm(charges ~ bmicatsobrepeso + bmicatobeso, data=seguro)
> summary(modelo2)

Call:
lm(formula = charges ~ bmicatsobrepeso + bmicatobeso, data = seguro)

Residuals:
     Min       1Q   Median       3Q      Max 
-21191.9  -6159.7   -429.2   5226.7  25572.1 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      10284.3      503.4   20.43   <2e-16 ***
bmicatsobrepeso   9002.1      880.4   10.22   <2e-16 ***
bmicatobeso      26736.4      763.5   35.02   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7912 on 554 degrees of freedom
Multiple R-squared:  0.6906,    Adjusted R-squared:  0.6895 
F-statistic: 618.3 on 2 and 554 DF,  p-value: < 2.2e-16

ANOVA

Una alternativa al análisis con variables dummies que hicimos recién es hacer un ANOVA usando como variable explicativa el bmi categorizado (bmicat) que defimos antes:

> modelo3 <- lm(charges ~ bmicat, data=seguro)
> summary(modelo3)
Call:
lm(formula = charges ~ bmicat, data = seguro)

Residuals:
     Min       1Q   Median       3Q      Max 
-21191.9  -6159.7   -429.2   5226.7  25572.1 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      10284.3      503.4   20.43   <2e-16 ***
bmicatsobrepeso   9002.1      880.4   10.22   <2e-16 ***
bmicatobeso      26736.4      763.5   35.02   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7912 on 554 degrees of freedom
Multiple R-squared:  0.6906,    Adjusted R-squared:  0.6895 
F-statistic: 618.3 on 2 and 554 DF,  p-value: < 2.2e-16

Como ya dijimos, la función que se usa para ajustar el modelo es la misma que en el caso de la regresión: lm. Entonces, pasando en limpio: la regresión y el ANOVA son casos particulares del modelo lineal que se diferencian únicamente en la naturaleza de las explicativas. Si las explicativas son continuas estamos en presencia de una regresión, mientras que si son categóricas estamos ante un ANOVA. Cuando hacemos un ANOVA implícitamente estamos haciendo una regresión con variables dummy, por eso los resultados que obtenemos en ambos casos son exactamente iguales (pueden chequearlo).

Prueba t de Student

Ahora que está clara la relación entre ANOVA y regresión…

pito
¿Qué pito toca el test-t en todo esto?

Recordarán que el test-t para dos muestras sirve para comparar las medias de dos poblaciones con distribución normal. Es lo que utilizaríamos, por ejemplo, si queremos comparar el efecto de un medicamento con un placebo, la cantidad de masa muscular ganada mediante dos planes de entrenamiento diferentes o la cantidad promedio de queso que llevan las pizzas italianas y las argentinas.

pizzas

Me ofrezco para este experimento

Siguiendo con nuestro dataset, podemos usar un test-t considerando solo a los sujetos normales y a los con sobrepreso (descartaremos de este análisis a los obesos).

> seguro2 <- seguro[seguro$bmicat != "obeso", ]

Ahora comparemos los resultados del test-t y de un ANOVA usando lm:

> t.test(charges ~ bmicat, data=seguro2, var.equal=TRUE)
        Two Sample t-test

data:  charges by bmicat
t = -11.203, df = 365, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10582.212  -7421.957
sample estimates:
   mean in group normal mean in group sobrepeso 
               10284.29                19286.37 

> summary(lm(charges ~ bmicat, data=seguro2))
Call:
lm(formula = charges ~ bmicat, data = seguro2)

Residuals:
   Min     1Q Median     3Q    Max 
 -9189  -6266  -1344   4358  24785 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      10284.3      459.5   22.38   <2e-16 ***
bmicatsobrepeso   9002.1      803.5   11.20   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7221 on 365 degrees of freedom
Multiple R-squared:  0.2559,    Adjusted R-squared:  0.2538 
F-statistic: 125.5 on 1 and 365 DF,  p-value: < 2.2e-16

Vemos que coinciden los valores reportados (t = 11.03, df = 365, p-value < 2.2e-16). Las medias estimadas también son idénticas 10284.29 para el grupo normal y 19286.37 para el grupo con sobrepeso (en el caso del lm este número sale de sumar los estimados del intercepto y el sobrepeso, 10284.3 + 9002.1 = 19286.37).

¿Por qué sucede esto? porque el test-t no es otra cosa que un caso particular del ANOVA (y por lo tanto del modelo lineal) donde la cantidad de tratamientos es igual a dos.

Como dice Vox Dei “todo concluye al fin nada puede escapar, todo tiene un final, todo termina” y así como quien no quiere la cosa llegamos al final. Nos vemos la próxima y si les gustó el post recuerden comentar / compartir / megustear.

Desde aquí pueden descargar el script completo para R.

Leave a comment