INTRODUCCIÓN

¿A qué llamamos Modelo Estadístico?

Podríamos decir que el modelado estadístico es una forma simplificada de entender (o “modelar”) un fenómeno a partir de un conjunto de datos. Hay diversas técnicas que nos permiten modelar datos y van desde las regresiones lineales que veremos a continuación (las más simples) hasta las redes neuronales (las más complejas).

Un modelo siempre se compone por una sola variable “dependiente” y una o más variables “independientes/explicativas” (o predictoras). Existen 2 tipos de modelos: explicativos y predictivos. El primero nos permiten analizar/explicar las relaciones entre las variables y el segundo a predecir comportamientos a futuro.

¿Qué es una Regresión Lineal?

La regresión lineal, donde la variable a explicar depende linealmente de las variables explicativas, es quizás el modelo estadístico más utilizado debido a su sencilla ejecución e interpretación en comparación a otros. Sin embargo, no siempre es el más adecuado y su uso depende de la complejidad del fenómeno a modelar.

La regresión lineal utiliza el método de Mínimos Cuadrados Ordinarios (MCO), uno de los más comunes a la hora de analizar la relación entre la variable dependiente y las explicativas. Lo que hace es encontrar la recta que mejor se ajuste a nuestros datos minimizando el error entre todas las observaciones.

Pero, ¿Cómo se encuentra esta recta? Esto se logra minimizando la suma del error cuadrático (suma de errores entre la predicción y el valor observado, elevado al cuadrado).

Hay 2 tipos:

  1. Regresión Lineal Simple, con 1 sola variable explicativa.

  2. Regresión Lineal Múltiple, con 2 o más variables explicativas.

Cabe destacar que, en ambos casos, las variables explicativas pueden ser numéricas o categóricas; sin embargo la variable dependiente es siempre numérica.

REGRESIÓN LINEAL SIMPLE

En R podemos hacer regresiones lineales simples con la función lm() de la librería tidyverse:

#install.packages(tidyverse)
library(tidyverse)

Con el objetivo de entender cómo funciona una regresión lineal simple y cómo se interpretan los coeficientes arrojados, vamos a trabajar con datos descargados de https://ourworldindata.org/covid-deaths que contienen información relacionada a fallecimientos por COVID-19.

Sin embargo, para facilitar el trabajo durante la clase, el dataset que utilizaremos fue previamente procesado y contiene información actualizada al 1 de marzo de 2021. Pueden descargarlo desde este link

Comencemos cargando nuestro dataset:

covid_data <- read.csv("data/covid_data.csv")

Y veamos que tiene:

dim(covid_data)
## [1] 173   9

Vemos que tenemos 173 registros y 9 columnas, pero ¿Qué información tiene cada columna? Usemos names para saberlo:

names(covid_data)
## [1] "iso_code"                 "location"                
## [3] "continent"                "total_cases"             
## [5] "total_cases_per_million"  "total_deaths"            
## [7] "total_deaths_per_million" "population"              
## [9] "aged_70_older"

Bien, podemos ver que cada registro pertenece a un país y que la información asociada refiere a: cantidad de casos de COVID, cantidad de muertes por COVID, porcentaje de adultos mayores de 70, y población.

Ahora veamos que tipo de dato contiene cada columna:

str(covid_data)
## 'data.frame':    173 obs. of  9 variables:
##  $ iso_code                : Factor w/ 173 levels "AFG","AGO","ALB",..: 1 3 48 2 7 5 6 8 9 10 ...
##  $ location                : Factor w/ 173 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ continent               : Factor w/ 6 levels "Africa","Asia",..: 2 3 1 1 4 6 2 5 3 2 ...
##  $ total_cases             : int  55733 107931 113255 20854 769 2112023 172216 28986 460849 234662 ...
##  $ total_cases_per_million : num  1432 37505 2583 635 7853 ...
##  $ total_deaths            : int  2444 1816 2987 508 14 52077 3195 909 8574 3223 ...
##  $ total_deaths_per_million: num  62.8 631 68.1 15.5 143 ...
##  $ population              : int  38928341 2877800 43851043 32866268 97928 45195777 2963234 25499881 9006400 10139175 ...
##  $ aged_70_older           : num  1.34 8.64 3.86 1.36 4.63 ...

Para analizar relaciones entre variables, comencemos trabajando con las variables numéricas “total_deaths_per_million” y “aged_70_older” e intentemos responder la siguiente pregunta:

¿Cuál es la relación que existe entre fallecimientos por COVID por millón de habitantes y porcentaje de población mayor de 70 años por país?

Hagamos un primer gráfico que relacione ambas variables: el scatter plot que vimos la clase anterior!

ggplot(covid_data)+
  geom_point(aes(x=aged_70_older, y=total_deaths_per_million))

A simple vista parece existir una relación positiva entre ambas variables, es decir que a mayor % de población +70 años en un país, mayor cantidad de fallecidos por millón de habitantes. Pero mejoremos el gráfico sumándole etiquetas y colores que nos ayuden en la interpretación:

ggplot(covid_data)+
  geom_point(aes(x=aged_70_older, y=total_deaths_per_million, colour=continent), size=2)+
  labs(title = "¿Cuál es la relación entre ambas variables?",
       subtitle = "Muertes por millón de hab. vs % población +70",
       x= "% población +70",
       y= "Muertes por millón de hab.",
       colour= "Continente")+
  theme_light()

Por ejemplo, ahora podemos ver que además de la relación positiva entre variables, en términos generales los países de Europa son los que presentan los valores más altos, y los de África los valores más bajos.

Pero dejemos las suposiciones de lado y armemos nuestro primer modelo con la función lm() que nos permitirá medir la relación existente entre nuestras variables.

Ojo! tengamos cuidado con las conclusiones que sacamos al ver los resultados porque estos modelos tienen una limitación: pueden probar correlación pero no causalidad entre 2 variables.

lm_simple <- lm(data=covid_data, total_deaths_per_million ~ aged_70_older)

Como verán, para hacer un modelo lineal con la función lm(), lo primero que debemos escribir es el dataset con el que vamos a hacer el modelo, en nuestro caso data=covid_data.

Luego debemos escribir una fórmula donde indicamos del lado izquierdo cuál es nuestra variable dependiente (a predecir) y del lado derecho cuál es nuestra variable independiente/explicativa. Ambas variables se separan con el símbolo ~.

Finalmente, el resultado que podemos ver dentro de lm_simple no es el clásico dataframe al que estamos acostumbrados, sino que es una lista que guardó distintos atributos estadísticos del modelo. Veamos de que se trata todo esto y aprendamos a interpretar los resultados:

lm_simple
## 
## Call:
## lm(formula = total_deaths_per_million ~ aged_70_older, data = covid_data)
## 
## Coefficients:
##   (Intercept)  aged_70_older  
##        -17.98          79.30

Bien, acá lo que nos interesa analizar son los “Coefficients” que arrojó nuestro modelo y nos indican el efecto esperado en la variable dependiente cuando varía una unidad de la variable independiente.

Entonces, el coeficiente que arroja la variable aged_70_older es 79,3 y esto significa que, para cada punto adicional de % de población mayor de 70 años por país, se registrarán 79 nuevos fallecidos por millón de habitantes.

A su vez, dentro de los coeficientes, aparece una variable llamada “(Intercept)” que hace referencia a la ordenada al origen o intersección con el eje Y de la línea recta que estamos calculando en nuestro modelo.

Pero, ¿Qué quiere decir esto? ¿Cómo hace la función lm() para calcular estos coeficientes? Esto ocurre porque, tal como su nombre lo indica, la regresión lineal encuentra la línea recta que mejor ajuste a todos los puntos/observaciones, disminuyendo la distancia a cada uno.

Para que esta definición sea más visual, probemos graficar la recta con geom_abline() indicando los coeficientes que arrojó el modelo: intersección (-17.98) y pendiente (79.30).

ggplot(covid_data)+
  geom_point(aes(x=aged_70_older, y=total_deaths_per_million, colour=continent), size=2)+
  geom_abline(aes(intercept = -17.98, slope = 79.30), color = "red", size=0.75)+
  labs(title = "¿Cuál es la relación entre ambas variables?",
       subtitle = "Muertes por millón de hab. vs % población +70",
       x= "% población +70",
       y= "Muertes por millón de hab.",
       colour= "Continente")+
  theme_light()

Notarán que los puntos son los mismos que en el gráfico anterior, pero la novedad es que pudimos agregar la línea recta que mejor ajusta a todas las observaciones, disminuyendo la distancia con cada una.

Pero recuerden que al principio dijimos que nuestro modelo también nos iba a ayudar a predecir/estimar nuevos valores. Veamos un ejemplo: Si quisiésemos entender rápidamente y de forma visual cuántos fallecidos por millón de habitantes tendría un país con un 30% de población mayor de 70 años, con solo extender los límites del gráfico podríamos averiguarlo. Para esto utilicemos xlim() y ylim():

ggplot(covid_data)+
  geom_point(aes(x=aged_70_older, y=total_deaths_per_million, colour=continent), size=2)+
  geom_abline(aes(intercept = -17.98, slope = 79.30), color = "red", size=0.75)+
  labs(title = "¿Cuál es la relación entre ambas variables?",
       subtitle = "Muertes por millón de hab. vs % población +70",
       x= "% población +70",
       y= "Muertes por millón de hab.",
       colour= "Continente")+
  theme_light() +
    xlim(c(0, 30)) +
    ylim(c(0, 3000))

Entonces, según nuestro modelo, si algún país tuviese un 30% de su población mayor de 70 años, los fallecidos por millón de habitantes serían aproximadamente 2400.

Y ¿Cómo calculo este valor de forma precisa? Bueno, pensemos que si cada un 1% de población mayor de 70 años, hay 79,3 nuevos fallecidos por millón de habitantes, con hacer (79,3*30)-17,98 llegaríamos al resultado: 2361. Sin embargo, también podemos calcularlo de forma rápida a partir de la función predict():

test_30 <- data.frame(aged_70_older = 30)

predict(lm_simple, test_30)
##        1 
## 2360.918

Hasta acá, graficamos la línea asignando los coeficientes de forma manual. Sin embargo, ggplot() también tiene una función llamada geom_smooth() que nos permite hacerlo de forma automática, solo asignando la variable dependiente y la independiente.

ggplot(covid_data)+
  geom_point(aes(x=aged_70_older, y=total_deaths_per_million, colour=continent), size=2)+
  geom_smooth(aes(x=aged_70_older, y=total_deaths_per_million), method = "lm", colour="red") +
  labs(title = "¿Cuál es la relación entre ambas variables?",
       subtitle = "Muertes por millón de hab. vs % población +70",
       x= "% población +70",
       y= "Muertes por millón de hab.",
       colour= "Continente")+
  theme_light()
## `geom_smooth()` using formula 'y ~ x'

Como podrán ver, la línea recta es exactamente la misma pero se calculó automáticamente. Sin embargo verán que la función también graficó un área grisada que muestra el intervalo de confiaza del 95% para la línea recta de regresión. Es decir que, que dentro del área gris se encuentra el valor real con un 95% de certeza.

Ahora analicemos los residuos del modelo

En las regresiones lineales, es muy importante analizar los residuos ya que hacen referencia a las diferencias que existen entre los valores que observamos y los que predecimos a partir de la recta. Es decir que cuanto menores sean los residuos, mejor se ajustará nuestro modelo. Para calcularlos en R, usamos la función residuals() de la siguiente forma:

residuos_lm_simple <- residuals(lm_simple)

Revisemos el resultado con head():

head(residuos_lm_simple)
##          1          2          3          4          5          6 
##  -25.25625  -36.34162 -219.74882  -74.56367 -206.27944  580.18896

Lo que obtuvimos fue una lista con todos los residuos (distancia a la recta) de cada una de las observaciones. Pero para poder analizarlos, calculemos estos valores directamente en una nueva columna del dataset original:

covid_data <- covid_data %>%
  mutate(residuos_lm = residuals(lm_simple))

Revisemos como quedan:

head(covid_data)
##   iso_code            location     continent total_cases
## 1      AFG         Afghanistan          Asia       55733
## 2      ALB             Albania        Europe      107931
## 3      DZA             Algeria        Africa      113255
## 4      AGO              Angola        Africa       20854
## 5      ATG Antigua and Barbuda North America         769
## 6      ARG           Argentina South America     2112023
##   total_cases_per_million total_deaths total_deaths_per_million population
## 1                1431.682         2444                   62.782   38928341
## 2               37504.691         1816                  631.038    2877800
## 3                2582.721         2987                   68.117   43851043
## 4                 634.511          508                   15.457   32866268
## 5                7852.708           14                  142.962      97928
## 6               46730.539        52077                 1152.254   45195777
##   aged_70_older residuos_lm
## 1         1.337   -25.25625
## 2         8.643   -36.34162
## 3         3.857  -219.74882
## 4         1.362   -74.56367
## 5         4.631  -206.27944
## 6         7.441   580.18896

Efectivamente, en la última columna se asignó un valor residual para cada registro (país). Cabe destacar que, mientras más cerca de 0 estén estos residuos, mejor se ajustará el modelo. Veamos esto gráficamente:

ggplot(covid_data) +
    geom_point(aes(x = aged_70_older, y = residuos_lm, colour=continent), size=2) +
    geom_hline(yintercept=0, colour="red", size=0.75) +
  labs(title = "Análisis de Residuos",
       x= "% mayores de 70",
       y= "Residuos del modelo")+
  theme_light()

A simple vista, podríamos decir que los países que más se ajustan a la recta son los de África y que los que están más alejados son los de Europa.

Pero no se asusten, es normal que las predicciones difieran de los valores reales. Sin embargo hay un aspecto a tener en cuenta: los residuos deben distribuirse de forma normal (con una media de cero) sobre la recta de regresión. Para revisar esto, hagamos un histograma:

ggplot(covid_data) +
  geom_histogram(aes(x=residuos_lm)) +
  labs(title = "Análisis de Residuos",
       x= "Residuos del Modelo")+
  theme_light()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

El histograma no reproduce una perfecta distribución normal, pero podemos justificar esto con el tamaño de la muestra que estamos trabajando: 173 observaciones. Probablemente, si la muestra fuese mayor, la distribución se ajustaría mejor.

Hasta acá solo analizamos la relación entre 2 variables numéricas, pero como indicamos al principio, también es posible utilizar regresiones lineales para analizar la relación entre una variable numérica y una categórica. Probemos que ocurre con la variables numérica “total_deaths_per_million” y la categórica “continente” e intentemos responder la siguiente pregunta:

¿Cuánto incide el continente en la cantidad de muertes por millón de habitantes de cada país?

Empecemos graficando un scatter plot:

ggplot(data = covid_data) +
    geom_point(aes(x = continent, y = total_deaths_per_million, color = continent)) +
    labs(title = "Fallecimientos por millón por continente",
         y = "Muertes por millón")

Rápidamente podemos notar que los países de Europa tienen la mayor cantidad de muertes por millón de habitantes. Pero, a simple vista parecería ser que un scatter plot no es el mejor gráfico para analizar esto, así que probemos con un boxplot o gráfico de cajas que sirve para comparar la distribución y tendencia central de una variable numérica para cada nivel de una variable categórica.

ggplot(data = covid_data) +
    geom_boxplot(aes(x=continent, y=total_deaths_per_million, fill=continent), show.legend = FALSE) +
    labs(title = "Fallecimientos por millón de hab. por continente",
         y = "Muertes por millón") +
  scale_fill_brewer(palette = "Set3")

Ahora podemos percibir con mayor detalle la relación entre los fallecimientos y el continente. Pero con este gráfico no nos alcanza para sacar conclusiones precisas así que hagamos un nuevo modelo lineal que incluya la variable continent como explicativa:

lm_simple_continente <- lm(data=covid_data, total_deaths_per_million ~ continent)

lm_simple_continente
## 
## Call:
## lm(formula = total_deaths_per_million ~ continent, data = covid_data)
## 
## Coefficients:
##            (Intercept)           continentAsia         continentEurope  
##                  91.72                  108.24                  961.24  
## continentNorth America        continentOceania  continentSouth America  
##                 334.76                  -80.52                  669.65

Tal como en el caso anterior, mi modelo calculó el coeficiente de la variable explicativa, pero ¿Por qué aparecen 6 valores en vez de 2? Bueno, porque la variable explicativa en este caso es categórica y tiene 6 niveles diferentes: África, América del Norte, América del Sur, Asia, Europa y Oceanía. Entonces mi modelo ya no calcula un único coeficiente, sino que calcula un coeficiente para cada nivel de la variable.

Si observan con detalle, verán que el Intercept fue calculado a partir del promedio de fallecidos por millón de habitantes en los países de África (en orden alfabético era la primer categoría). Y el resto de coeficientes se calcularon a partir de la diferencia de fallecidos por millón de habitantes respecto a África. Es decir que, por ejemplo, en promedio los países de Asia tienen 108 fallecidos por millón de habitantes más que los países de África; en cambio Oceanía tiene en promedio 80 fallecidos menos.

Calculemos nuevamente los residuos:

covid_data <- covid_data %>%
  mutate(residuos_lm_continente = residuals(lm_simple_continente))
ggplot(covid_data) +
    geom_point(aes(x = continent, y = residuos_lm_continente, colour=continent), size=2) +
    geom_hline(yintercept=0, colour="red", size=0.75) +
  labs(title = "Análisis de Residuos",
       x= "Continente",
       y= "Residuos del modelo")+
  theme_light()

Podemos ver que en Norte América hay algunos paises que tiene un valor muy por encima de la recta estimada, y en Europa hay algunos que están muy por debajo. Veamos de que países se trata:

covid_data %>% 
    filter(continent == "North America") %>% 
    arrange(desc(residuos_lm_continente)) %>% 
    head(3)
##   iso_code      location     continent total_cases total_cases_per_million
## 1      USA United States North America    28637235                86516.63
## 2      MEX        Mexico North America     2089281                16204.42
## 3      PAN        Panama North America      341420                79128.24
##   total_deaths total_deaths_per_million population aged_70_older residuos_lm
## 1       514810                 1555.305  331002647         9.732    801.5713
## 2       186152                 1443.791  128932753         4.321   1119.1315
## 3         5858                 1357.663    4314768         5.030    976.7822
##   residuos_lm_continente
## 1              1128.8189
## 2              1017.3049
## 3               931.1769

Los valores más extremos de Norte América son Estados Unidos, México y Panamá.

covid_data %>% 
    filter(continent == "Europe") %>% 
    arrange(residuos_lm_continente) %>% 
    head(3)
##   iso_code location continent total_cases total_cases_per_million total_deaths
## 1      ISL  Iceland    Europe        6054                17740.66           29
## 2      NOR   Norway    Europe       71735                13232.21          623
## 3      FIN  Finland    Europe       58064                10479.51          750
##   total_deaths_per_million population aged_70_older residuos_lm
## 1                   84.982     341250         9.207   -627.1209
## 2                  114.918    5421242        10.813   -724.5354
## 3                  135.362    5540718        13.264   -898.4475
##   residuos_lm_continente
## 1              -967.9782
## 2              -938.0422
## 3              -917.5982

Los valores más extremos de Europa son Islandia, Noruega y Finlandia.

Bueno hasta acá, hicimos 2 modelos lineales simples donde pudimos ver como por separado las variables aged_70_older y continent logran explicar la variable total_deaths_per_million. ¿Pero qué ocurre si en un mismo modelo incluimos a la vez 2 variables explicativas?

REGRESIÓN LINEAL MÚLTIPLE

Este tipo de regresiones nos permite generar un modelo lineal (como el que vimos en el apartado anterior) pero la diferencia está en que, en vez de 1, utiliza 2 o más variables explicativas/predictoras. ¿Por qué haríamos esto? Para poder representar de la mejor forma posible las diferentes cuestiones que inciden en el caso de estudio que en este caso es la cantidad de fallecidos por COVID19 por millón de habitantes.

Entonces, teniendo esto en mente, veamos cuanto incide el % de población mayor de 70 y el continente en la cantidad de fallecidos. Para esto, sigamos utilizando lm() pero ahora agreguemos un signo + entre ambas variables explicativas. La fórmula sería algo así: Y ~ X1 + X2

¿Cuánto incide el porcentaje de población mayor de 70 años y la cantidad de contagios por millón de habitantes en la cantidad de muertes por millón de habitantes de cada país?

lm_multiple <- lm(data=covid_data, total_deaths_per_million ~ aged_70_older + total_cases_per_million)

lm_multiple
## 
## Call:
## lm(formula = total_deaths_per_million ~ aged_70_older + total_cases_per_million, 
##     data = covid_data)
## 
## Coefficients:
##             (Intercept)            aged_70_older  total_cases_per_million  
##               -54.81533                 28.08490                  0.01428

Como veniamos viendo en la regresión simple, primero tenemos un coeficiente de intersección o coordenada al origen que no nos brinda demasiada información del modelo y luego vemos coeficientes asignados a ambas variables explicativas: aged_70_older y continent. ¿Cómo podemos leer estos coeficientes?

  • aged_70_older: En países con cantidades similares de contagios por millón de habitantes, el incremento de un 1% de población mayor de 70 años, generará 28 nuevas muertes por millón de habitantes.

  • total_cases_per_million: En países con porcentajes similares de población mayor a 70 años, el incremento de un 1 nuevo contagio por millón de habitantes, generarán 0,01428 muertes por millón de habitantes. O dicho en otras palabras, cada 70 nuevos contagios por millón de habitantes habrá 1 fallecido por millón de habitantes.

Hasta acá solo aprendimos a interpretar coeficientes, pero si queremos entender cuales de las variables son buenas predictoras y cuales no, debemos usar la función summary():

summary(lm_multiple)
## 
## Call:
## lm(formula = total_deaths_per_million ~ aged_70_older + total_cases_per_million, 
##     data = covid_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -773.10 -106.66   -2.35   49.25 1145.87 
## 
## Coefficients:
##                            Estimate  Std. Error t value             Pr(>|t|)
## (Intercept)             -54.8153314  34.2127398  -1.602                0.111
## aged_70_older            28.0849003   6.0502825   4.642           0.00000687
## total_cases_per_million   0.0142792   0.0009931  14.378 < 0.0000000000000002
##                            
## (Intercept)                
## aged_70_older           ***
## total_cases_per_million ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.5 on 170 degrees of freedom
## Multiple R-squared:  0.7342, Adjusted R-squared:  0.7311 
## F-statistic: 234.8 on 2 and 170 DF,  p-value: < 0.00000000000000022

La función summary() nos muestra un resumen estadístico de nuestro objeto (un modelo lineal múltiple) que incluye:

  • Valores relacionados a los residuos del modelo: mínimo, máximo, mediana, cuartiles. Como vimos anteriormente, los residuos exponen la diferencia entre la predicción del modelo y los resultados reales, por lo tanto, mientras más bajos sean, mejor es nuestro modelo.

En nuestro modelo podemos ver que los residuos se encuentran en un rango de -773,10 a 1145,87, y que la mediana está muy cerca de 0.

  • Coeficientes calculados para cada variable: estimación del efecto, error estándar, valor t y valor p.

Con error estándar (Std. Error) nos referimos a la desviación estándar de cada una de las estimaciones. Con este valor, que siempre es positivo, vamos a medir la precisión con la que el modelo realizó las estimaciones. Es decir que, cuanto menor sea el error estándar, más precisa será la estimación.

Luego, a partir de dividir el valor estimado con el error estándar se genera el valor t (t value).

Y por último si está el valor p o Pr(>|t|) que puede ir de 0 a 1. Si es menor a 0,5 (valor default), vamos a considerar que la capacidad predictiva de la variable es estadísticamente significativa; y si es mayor, vamos a considerar que no es estadísticamente significativa.

Pero más allá de saber interpretar el valor p, R nos brinda la información a partir de asteriscos: Si se fijan con atención, para cada valor p hay una cantidad de asteriscos asignada que va de " " o “.” (alto valor p y baja significancia estadística) a *** (bajo valor p y alta significancia estadística).

Cabe destacar que, la presencia de un solo asterisco ya nos muestra significancia estadística debido a que nos está indicando que el valor p es igual o mayor al 95% de confianza (valor default).

Puntualmente, y hablando de nuestro modelo lineal múltiple, podemos ver que tanto el porcentaje de población mayor de 70 años, como la cantidad de contagios por millón de habitantes son estadísticamente significativos por tener asignados 3 * y presentar un valor p < 0,05.

Pero esto no es todo, también el resumen estadístico nos arroja 4 medidas más: el “Residual standard error”, el “Multiple R-squared”, el “Adjusted R-squared” y la “F-statistic”:

El error estándar residual se calcula con la desviación estándar de los residuos, y mantiene la misma lógica: cuanto más pequeño, mejor.

El R2 múltiple va de 0 a 1 y explica que tan bien se ajusta el modelo a los datos. Cuanto mayor sea el R2, mejor se ajustará el modelo a los datos. Se utiliza principalmente para regresiones lineales simples (una sola variable indepentiente).

El R2 ajustado normaliza el R2 múltiple, teniendo en cuenta el número de variables independientes. Por lo tanto, es útil en regresiones lineales múltiples. En nuestro modelo el valor es bastante alto: 0,73, lo cual indica que el modelo se ajusta bastante bien a los datos.

En nuestro caso, tenemos un R2 de 0,73 por lo que podemos concluir que la regresión lineal múltiple realizada se ajusta de forma aceptable a nuestros datos.

La estadística F tiene en cuenta la cantidad de variables y observaciones utilizadas y verifica si al menos una de las variables independientes tiene significancia estadística.

El p valor permite medir la significancia estadística del análisis completo. Por defecto, suele realizarse bajo un nivel de confianza del 95%, o dicho en otras palabras, una significancia por debajo de 0,05.

Hasta acá estamos utilizando 2 variables explicativas, pero aprovechando que no nos cuesta esfuerzo adicional, probemos incluyendo una 3er variable explicativa a nuestro modelo e intentemos responder la siguiente pregunta:

¿Cuánto incide el porcentaje de población mayor de 70 años, el total de contagios por millón de habitantes y el continente en la cantidad de muertes por millón de habitantes de cada país?

En este caso, la fórmula del modelo lineal sería algo así: Y ~ X1 + X2 + X3

lm_multiple_2 <- lm(data=covid_data, total_deaths_per_million ~ aged_70_older + total_cases_per_million + continent)

summary(lm_multiple_2)
## 
## Call:
## lm(formula = total_deaths_per_million ~ aged_70_older + total_cases_per_million + 
##     continent, data = covid_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -605.14 -101.37   -8.68   83.26 1092.81 
## 
## Coefficients:
##                            Estimate  Std. Error t value             Pr(>|t|)
## (Intercept)               -4.525738   39.290119  -0.115             0.908436
## aged_70_older             14.836814    8.479488   1.750             0.082024
## total_cases_per_million    0.013321    0.001021  13.053 < 0.0000000000000002
## continentAsia           -112.878924   55.887493  -2.020             0.045027
## continentEurope          191.216048  104.876706   1.823             0.070076
## continentNorth America    75.540210   72.604564   1.040             0.299661
## continentOceania         -84.173506  136.430645  -0.617             0.538105
## continentSouth America   319.683075   87.478330   3.654             0.000346
##                            
## (Intercept)                
## aged_70_older           .  
## total_cases_per_million ***
## continentAsia           *  
## continentEurope         .  
## continentNorth America     
## continentOceania           
## continentSouth America  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 254.1 on 165 degrees of freedom
## Multiple R-squared:  0.7774, Adjusted R-squared:  0.768 
## F-statistic: 82.32 on 7 and 165 DF,  p-value: < 0.00000000000000022

Al incluir variables categóricas como predictoras, tenemos que tener en cuenta que una categoría se considera como referencia (Intercept) y el resto de categorías se comparan con él (en este caso es África).

Por un lado, vemos que la regresión múltiple realizada, a diferencia de la anterior, no considera la variable aged_70_older como estadísticamente significativa. ¿A qué se debe esto? A que mi modelo está “más completo” y tiene en cuenta más variables, por lo tanto los valores de la variable en cuestión ya no son significativos luego de conocer el valor de la variable continente. Es decir que, como mi modelo conoció otros datos que pueden tener influencia en la cantidad de fallecimientos por millón de habitantes, dejó de considerar estadísticamente significativa a la variable aged_70_older.

Por otro lado, observamos que tanto la variable total_cases_per_million como las categorías Asia y Sudamérica incluidas en la variable categórica continente son significativas. Esto último quiere decir que, con las demas variables fijas, los países de Sudamérica tendrán en promedio 319 más muertes por millón de habitantes que África, mientras que los de Asia tendrán 112 muertes menos.