Hasta aquí hemos usado la regresión lineal para hacer explícita la relación entre una variable resultante \(y\) y una única variable predictiva o explicativa \(x\). En algunos de nuestros resultados pudimos intuir que el agregado de alguna variable explicativa adicional podría mejorar nuestras predicciones. De eso se trata la regresión lineal múltiple: incorporar una cantidad arbitraria de variables al modelo, buscando representar las múltiples dinámicas que inciden en el fenómeno estudiado.
Una buena noticia es que, en general, agregar variables a nuestro modelo estadístico no requiere mucho esfuerzo adicional. En la época en que los cálculos matemáticos debían hacerse sin la ayuda de una computadora, sumar variables sin ton ni son debía tener poca gracia, debido a la creciente cantidad de cálculos a resolver. Para nosotros que dejamos la tarea en manos de software especializado, el problema es el opuesto. Es tan fácil sumar variables al modelo, que debemos evitar la tentación de arrojar todo dentro de la fórmula de regresión lineal y decidir luego que parece importante y que no.
Pasemos a la práctica. Vamos a modelar la expectativa como resultante de la población y del PBI per cápita de los países, usando los datos más reciente (tomados en 2007). La única diferencia respecto a una regresión lineal simple es que usamos + para agregar variables en la fórmula de lm()
##
## Call:
## lm(formula = expVida ~ pobl + PBI_PC, data = data_mundial_2007)
##
## Coefficients:
## (Intercept) pobl PBI_PC
## 59.205198140717 0.000000007001 0.000641608517¿Cómo interpretamos esos resultados? Más o menos de la misma manera que con la regresión simple. Como antes, tenemos un coeficiente para la intersección, al que no prestamos mucha atención porque no nos dice nada de la relación entre las variables. Lo que cambia es que esta vez tenemos dos variables predictoras en lugar a una, cada una con su coeficiente. Los coeficientes positivos indican que la relación de la población con la expectativa de vida es de correlación positiva (cuando una crece la otra tiende a crecer también), y lo mismo ocurre con el PBI. La magnitud de los coeficientes es pequeña (minúscula en el caso de la población), lo cual dificulta “narrar” los resultados, pero podemos hacerlo así:
Pensemos un poco si los resultados tienen sentido. La correlación positiva entre PBI y longevidad es de lo más razonable. No nos extraña que los países de mayores ingresos tiendan a ser aquellos cuyos habitantes viven más tiempo. La correlación con la población es quizás inesperada. Si la longevidad se incrementa junto a la cantidad de habitantes, ¿acaso no deberíamos encontrar a varios de los países más populosos entre los más longevos?
Veamos el top ten de países más poblados:
## pais continente anio expVida pobl PBI_PC residuo_ml
## 1 Japan Asia 2007 82.603 127467972 31656.07 11.87452
## 2 Hong Kong, China Asia 2007 82.208 6980412 39724.98 11.47952
## 3 Iceland Europe 2007 81.757 301931 36180.79 4.10840
## 4 Switzerland Europe 2007 81.701 7554661 37506.42 4.05240
## 5 Australia Oceania 2007 81.235 20434176 34435.37 0.51550
## 6 Spain Europe 2007 80.941 40448191 28821.06 3.29240
## 7 Sweden Europe 2007 80.884 9031088 33859.75 3.23540
## 8 Israel Asia 2007 80.745 6426679 25523.28 10.01652
## 9 France Europe 2007 80.657 61083916 30470.02 3.00840
## 10 Canada Americas 2007 80.653 33390141 36319.24 7.04488y el de países con mayor expectativa de vida:
## pais continente anio expVida pobl PBI_PC residuo_ml
## 1 China Asia 2007 72.961 1318683096 4959.115 2.23251515
## 2 India Asia 2007 64.698 1110396331 2452.210 -6.03048485
## 3 United States Americas 2007 78.242 301139947 42951.653 4.63388000
## 4 Indonesia Asia 2007 70.650 223547000 3540.652 -0.07848485
## 5 Brazil Americas 2007 72.390 190010647 9065.801 -1.21812000
## 6 Pakistan Asia 2007 65.483 169270617 2605.948 -5.24548485
## 7 Bangladesh Asia 2007 64.062 150448339 1391.254 -6.66648485
## 8 Nigeria Africa 2007 46.859 135031164 2013.977 -7.94703846
## 9 Japan Asia 2007 82.603 127467972 31656.068 11.87451515
## 10 Mexico Americas 2007 76.195 108700891 11977.575 2.58688000El único país presente en ambas listas es Japón. Ni nuestro conocimiento del mundo, ni los datos parecen apoyar la noción de que población y longevidad van juntos. Ya hemos usado cor() para obtener una medida de la intensidad de la correlación entre dos variables. Veamos que pasa con longevidad vs. población:
## [1] 0.04755312Recordemos que la intensidad de una correlación es su valor absoluto, que toma un máximo de 1, mientras que el signo (positivo o negativo) indica si la relación entre variables es directa o inversa. Aquí obtuvimos un valor bien bajo, cercano a cero: la correlación es nula. Entonces ¿Por qué aparece en nuestro modelo de regresión lineal?
En resumidas cuentas, aparece porque nosotros le pedimos que aparezca. Es decir, instruimos en forma específica a lm() para que incorpore a la población en el modelo. El caso es que población no es un buen predictor de longevidad (la correlación es bajísima), pero si lo pedimos, lo tenemos: el coeficiente nos indica el valor que minimiza las discrepancias entre valores observado y valores predichos trazando una línea recta. Lo que no indica por si solo es el grado en el cual podemos confiar en esa variable para darnos buenas predicciones o estimados.
Sería muy útil que el resultado de lm() indique cuáles variables son buenas predictoras y cuáles no. Y por suerte, lo hace cuando lo interrogamos con summary(), la misma función que hemos estado usando para obtener el resumen de un dataframe. Cuando la usamos con un objeto de R que contiene un modelo estadístico, lo que obtenemos son sus detalles:
##
## Call:
## lm(formula = expVida ~ pobl + PBI_PC, data = data_mundial_2007)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.496 -6.119 1.899 7.018 13.383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.205198140717 1.040398672164 56.906 <0.0000000000000002 ***
## pobl 0.000000007001 0.000000005068 1.381 0.169
## PBI_PC 0.000641608517 0.000058176209 11.029 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.87 on 139 degrees of freedom
## Multiple R-squared: 0.4679, Adjusted R-squared: 0.4602
## F-statistic: 61.11 on 2 and 139 DF, p-value: < 0.00000000000000022El resumen incluye los parámetros que definieron al modelo, los valores por cuartil de los residuos, y una tabla con variables numéricas. En esa tabla, bajo la columna Estimate tenemos el “efecto” estimado de cada variable explicativa sobre la dependiente. Es decir, los coeficientes que ya conocemos. Luego aparecen tres columnas con atributos estadísticos: Std. Error, t value, y Pr(>|t|). En castellano las llamaríamos, respectivamente, error estándar, valor t y valor p. Interpretar estos valores cae fuera de nuestros objetivos, pero podemos señalar que el más famoso entre ellos es el valor p, porque se usa como medida: si vale menos de 0,5, se considera que la capacidad de predicción de la variable asociada es significativa. Para interpretar todo esto de manera sencilla, una vez más vamos a confiar en R para guiarnos. He aquí la curiosa forma de determinar si una variable es buena predictora o no: contar estrellitas. Junto a cada fila aparecen, a veces, de uno a tres asteriscos. Son la forma de R de decirnos cuales son las variables explicativas que muestran una relación “estadísticamente significativa” con nuestra variable dependiente. Cuanto más bajo el valor p, más significativa es la relación y más estrellitas aparecen:
. o nada: No se encuentra una relación entre esta variable y la que queremos predecir.*: Es muy probable que esta variable tenga una relación con la que queremos predecir. Ya podemos publicar estos resultados en un paper científico.**: Es muy, pero muy probable que esta variable tenga una relación con la que queremos predecir. 99% seguro.***: Juramos que las variables están relacionadas. Más no se puede pedir.Lo de un asterisco/estrella (*) indicando que los resultados ya alcanzan rigor científico no es broma. El asterisco solitario indica que, a nivel estadístico, se supera el 95% de confianza en que la relación existe en la realidad y no es producto de una casualidad en los datos. Pasando ese umbral se considera que los datos son “estadísticamente significativos”, y desde hace muchos años encontrar un valor p menor a 0,05 es la meta dorada de los investigadores que emplean análisis estadístico. ¿Porqué un 95% de confianza alcanza? ¿Porqué no relajar el límite a 90%, o quizás mejor, exigir al menos un 99 o 99,9% de seguridad? La verdad es que no hay ninguna razón trascendental. El 95% de certeza es tan sólo un umbral arbitrario que en algún momento se volvió estándar. Es importante aclarar que en los últimos años ha crecido una reacción de rechazo a esta norma arbitraria, dentro de la propia comunidad científica. Quienes siguen confiando en los valores p son llamados “frecuentistas”; los que proponen cuantificar de otra forma nuestro grado de certeza son llamados “bayesianos”. Google mediante, quien quiera saber más sobre la apasionante rivalidad tendrá horas de diversión aseguradas.
En lo que a nosotros respecta, por ahora vamos a aceptar el enfoque frecuentista, y cuando veamos una estrella diremos que la variable asociada es un buen predictor. O para ser más precisos, que su relación con la variable dependiente es estadísticamente significativa.
Volvamos a nuestros modelos. Cuando hicimos regresiones simples no sabíamos aún de valores p, y no revisamos la significancia de las variables predictoras. Hagámoslo ahora con el modelo de expectativa de vida en Argentina vs. PBI :
##
## Call:
## lm(formula = expVida ~ anio, data = data_arg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53006 -0.13516 -0.01219 0.14228 0.55202
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -389.606345 9.677730 -40.26 0.000000000002140 ***
## anio 0.231708 0.004889 47.40 0.000000000000422 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2923 on 10 degrees of freedom
## Multiple R-squared: 0.9956, Adjusted R-squared: 0.9951
## F-statistic: 2246 on 1 and 10 DF, p-value: 0.0000000000004216Las tres estrellitas, distinción máxima, indican que sin dudas el año está relacionado con la expectativa de vida. Esto no es una sorpresa: la linea de la regresión lineal se ajusta con tanta precisión a los valores observados, que no podía ser de otra manera.
Continuando con las regresiones múltiples, intentemos un modelo con tres variables predictoras. A población y PBI, las que ya teníamos en cuenta, vamos a agregar una variable categórica: el continente.
modelo_exp_multiple <- lm(expVida ~ pobl + PBI_PC + continente, data = data_mundial_2007)
summary(modelo_exp_multiple)##
## Call:
## lm(formula = expVida ~ pobl + PBI_PC + continente, data = data_mundial_2007)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.8199 -2.8905 0.1574 2.9046 20.0585
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 53.7141900516204 0.9355709763972 57.413 < 0.0000000000000002 ***
## pobl 0.0000000009586 0.0000000039259 0.244 0.80747
## PBI_PC 0.0003479123814 0.0000571704015 6.086 0.00000001127738 ***
## continenteAmericas 16.0313726693021 1.6713252557392 9.592 < 0.0000000000000002 ***
## continenteAsia 12.5640427449841 1.6209815371922 7.751 0.00000000000197 ***
## continenteEurope 15.1989177617593 1.9662500363509 7.730 0.00000000000220 ***
## continenteOceania 16.6222095573924 4.9925674316223 3.329 0.00112 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.597 on 135 degrees of freedom
## Multiple R-squared: 0.7141, Adjusted R-squared: 0.7014
## F-statistic: 56.2 on 6 and 135 DF, p-value: < 0.00000000000000022Observamos que la variable categórica es significativa. Con las demás variables fijas -es decir, en países de similar PBI y población- el continente de origen explica en gran medida las diferencias en expectativa de vida en cada país, y con un efecto estimado enorme - ¡de 12 a 16 años!-. Notemos de todos modos que el coeficiente de la variable continente había sido mayor en el modelo simple, llegando a casi 26 años para Oceanía. ¿Porqué es menor ahora? Porque nuestro modelo es más completo, y tiene en cuenta más variables. Cuando lo único que teníamos para comparar países era su continente, era era la única variable a la que atribuir diferencias. Ahora que consideramos múltiples variables para explicar las diferencias, notamos la parte de la influencia que se lleva el PBI, reduciendo la del continente.