Cómo programar una regresión lineal en R

Existen muchos modelos que podemos programar en R, pero entre ellos el más básico, y uno de los más importantes, es la regresión lineal. Seguro que has oído hablar de ella, pero… ¿conoces exactamente cómo funciona y las presunciones que realiza?

No te preocupes, que en este post vamos a aprender cómo programar una regresión lineal en R desde 0, aprendiendo todo lo que tienes que saber sobre este modelo. ¡Vamos a ello!

Entendiendo el funcionamiento de la regresión lineal (con R)

Supongamos que un país no conoce su esperanza de vida y la quiere predecir. Para ello, utilizará el PIB per Cápita ya que, es un valor que sí puede calcular . Además, a primera vista parece que son variables relacionadas: es de suponer que países con mayor PIB per Cápita (más ricos), tendrán mayor esperanza de vida (seguramente derivado de otros factores como la sanidad), por lo que puede ser un buen estimador.

Para este ejemplo, vamos a usar los datos del 2007 de la librería gapminder. Asimismo, en vez de usar el Pib per Cápita, usaré el logaritmo del PIB per Cápita (lo explico más adelante en el post). En cualquier caso, veamos si se da esa relación:

library(gapminder)

year = 2007
data = gapminder[gapminder$year == year, ]
data$loggdpPercap = log(data$gdpPercap)

plot(data$lifeExp ~  data$loggdpPercap)

Efectivamente parece que hay una relación: visualmente vemos cómo podríamos dibujar una línea en la que, nos devuelva el valor de la esperanza de vida para diferentes valores del logaritmo del PIB per Cáptia.

Sin embargo, aunque podamos dibujar muchas líneas, solo hay una que se ajusta a los datos de la mejor manera posible. Y, qué mejor manera de verlo, que creando tres funciones lineales y viendo como no vale cualquier recta:


# Creo las funciones
regres_azul = function(x){67.00742 + x*0}
regres_rojo = function(x){96.666 - 3.333*x}
regres_amarillo = function(x){15.833 + 5.833*x}

# Calculo los datos predichos por dichas funciones
datos_azul = sapply(c(5:11), regres_azul)
datos_verde = sapply(c(5:11), regres_amarillo)
datos_rojo = sapply(c(5:11), regres_rojo)

# Creo el gráfico
plot(data$lifeExp ~ data$loggdpPercap)
lines(c(5:11), datos_azul, type = "l", col = "blue", lty = 2)
lines(c(5:11), datos_rojo, type = "l", col = "red", lty = 2)
lines(c(5:11), datos_amarillo, type = "l", col = "yellow3")

Como vemos en el ejemplo, las líneas azul y rojas no parece que serían muy buenas a la hora de predecir la esperanza de vida, mientras que la línea naranja sí que parece sería mejor.

Vemos claro que hay ciertas regresiones que se ajustan más a nuestros datos que otras. De hecho, hay una regresión que se ajusta mejor que cualquier otra regresión que hagamos. Y es que, en la regresión lineal, como en cualquier otro modelo, siempre hay, al menos, un óptimo.

Pero, ¿cómo decide el algoritmo cómo pintar la línea? De hecho, ¿cómo sabemos qué regresión es mejor? Para ello, primero hay que conocer los residuos de nuestra regresión. ¡Vamos a ello!

Calcular los residuos de la regresión lineal

La forma de saber qué regresión funciona mejor, es calcular cómo de cerca está la regresión de nuestros datos. A esa distancia entre el valor real y el que predice la regresión, se le llama residuo. Como es lógico, cuanto mejor se ajuste la regresión a los datos, menor será la distancia entre el valor predicho y el real. Es decir, que a mejor ajuste del modelo, menor serán los residuos.

Veamos un ejemplo gráfico de a qué me refiero mostrando los residuos de la regresión azul.

library(ggplot2)

data$predicted_azul = sapply(data$loggdpPercap, regres_azul)

ggplot(data, aes(loggdpPercap, lifeExp)) +
  geom_point() + 
  geom_segment(aes(xend = loggdpPercap, yend = predicted_azul)) +
  geom_point(aes(y = predicted_azul), shape = 1) +
  theme_minimal()
Residuos de un modelo de regresion lineal simple hecho en R utilizando la media.

Como vemos, la regresión azul tiene unos residuos muy grandes: hay mucha distancia entre los valores predichos (sin relleno negro) y los reales (con relleno negro).

Podríamos hacer lo mismo con la regresión amarilla. Antes hemos visto como, visualmente, se ajustaba mejor a los datos, por lo que es de esperar que los residuos deben ser menores:

data$predicted_amarillo = sapply(data$loggdpPercap, regres_amarillo)

ggplot(data, aes(loggdpPercap, lifeExp)) +
  geom_point() + 
  geom_segment(aes(xend = loggdpPercap, yend = predicted_amarillo)) +
  geom_point(aes(y = predicted_amarillo), shape = 1) +
  theme_minimal()
Residuos de un modelo de regresion lineal con un buen ajuste a nivel visual.

Efectivamente, visualmente vemos como los residuos son más bajos. Pero, ¿cómo sabemos esto matemáticamente? Para ello, una opción podría ser sumar todos los residuos de la regresión, ¿no? Veamos.

cat(" Residual Sum Azul: ", sum(data$lifeExp - data$predicted_azul), "\n",
"Residual Sum Amarillo: ",sum(data$lifeExp - data$predicted_amarillo))
 Residual Sum Azul:  0.00036 
 Residual Sum Amarillo:  130.4315

¡Problema! Nuestra función amarilla ajusta mejor que la azul, pero… ¡tiene más error! ¿Qué está pasando?

La respuesta es simple: la suma de residuos no es una buena forma de calcular el ajuste del modelo, ya que los errores positivos se compensan con los negativos.

Una forma de evitar este problema, es elevar los residuos al cuadrado antes de sumarlos. Así, los valores no se compensaran, haciendo que sea un calculo más certero. De hecho, es un valor muy utilizado, el cual se le conoce como la suma de residuos al cuadrado, o RSS (Residual Sum of Squares).

Veamos si la suma de errores al cuadrado indica lo que visualmente vemos: que la regresión amarilla ajusta mejor que la regresión azul.

rss = function(real, predicted){
  sum((real - predicted)^2)
}

cat(" RSS Azul: ", rss(data$lifeExp, data$predicted_azul), "\n",
"RSS Amarillo: ", rss(data$lifeExp, data$predicted_amarillo))
 RSS Azul:  20551.85 
 RSS Amarillo:  7707.969

¡Efectivamente! La suma de residuos al cuadrado sí sirve para conocer qué modelo ajusta mejor. De hecho, ya sabemos qué es lo que tenemos que buscar: la recta que minimiza la suma de residuos al cuadrado. Pero… ¿cómo la encontramos? Pues… en la realidad hay muchos métodos, pero el más utilizado es el método de Mínimos Cuadrados Ordinarios. ¡Veamos cómo funciona!

Mínimos Cuadrados Ordinarios (MCO)

El método de Mínimos Cuadrados Ordinales es un método que permite estimar los parámetros de una regresión lineal bajo ciertos supuestos (que veremos más adelante). Si se cumplen esos supuestos, el método de Mínimos Cuadrados Ordinaros proporcionará un estimador insesgado de varianza mínima. Vamos, lo que estamos buscando.

Según este método, dada una regresión \(Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\), podemos estimar los parámetros de nuestra regresión de la siguiente manera:

  • Pendiente de la Recta: \(\hat{\beta}_1 = \frac{\sum(X_i – \bar{X}) (Y_i – \bar{Y})} {\sum(X_i – \bar{X})^2}\). Se refiere a cuánto aumenta (o disminuye) el valor predicho con cada incremento en 1 de nuestra variable independiente. En nuestro caso, se refiere a cuánto aumenta (o disminuye) la esperanza de vida, con el incremento en 1 del logaritmo del PIB per Capita.
  • Ordenada (o intercept): \(\hat{\beta}_0 = \bar{Y} – \hat{\beta}_1 \bar{X}\). Se refiere cuánto vale la variable predicha cuando la variable independiente es cero. En nuestro caso, se refiere a cuál sería la esperanza de vida de un país cuyo logaritmo del PIB per Capita es 0.

Así pues, vamos a crear una función para estimar los parámetros. Empezamos con la pendiente.

pred_slope = function(x, y){
  x_mean = mean(x)
  y_mean = mean(y)

  sum((x - x_mean)*(y - y_mean))/sum((x - x_mean)^2)
}

slope = pred_slope(data$loggdpPercap, data$lifeExp)
slope
[1] 7.202802

Ahora que ya hemos estimada la pendiente, podemos estimar la ordenada:

pred_intercept = function(x, y, slope){
  y_mean = mean(y)
  x_mean = mean(x)

  y_mean - slope*x_mean

}

intercept = pred_intercept(data$loggdpPercap, data$lifeExp, slope)
intercept
[1] 4.949612

Y con esto solo hay que juntar ambos parámetros y… ¡ya tendríamos estimada nuestra regresión lineal! ¡Hagámoslo y veamos el residuos genera!

regres_lineal = function(x, slope, intercept){
  intercept + slope*x
}

data$predicted = sapply(data$loggdpPercap, regres_lineal,slope, intercept)

ggplot(data, aes(loggdpPercap, lifeExp)) +
  geom_point() + 
  geom_segment(aes(xend = loggdpPercap, yend = predicted)) +
  geom_point(aes(y = predicted), shape = 1) +
  theme_minimal()
Residuos del modelo de regresion lineal optimo simple hecho en R para predecir la esperanza de vida según el logaritmo del PIB Per Capita

Como vemos, parece que se ajusta mejor a todos los puntos, aunque sigue habiendo unos valores un tanto anómalos sobre los que está un poco lejos. En cualquier caso, este ajuste lo podemos comprobar tanto visualmente como con la suma de errores al cuadrado:


datos_predicho = sapply(c(5:11), regres_lineal, slope, intercept)

plot(data$lifeExp ~ data$loggdpPercap)
lines(c(5:11), datos_predicho, type = "l", col = "black")
Representación del modelo de regresion lineal optimo simple hecho en R para predecir la esperanza de vida según el logaritmo del PIB Per Capita

Como vemos, la recta ajusta muy bien, sobre todo para el cuadrante superior derecho, donde hay más datos.

Ya hemos aprendido a programar una regresión lineal en R. Ahora bien. Recuerda que esta predicción solo sirve si se cumplen ciertos supuestos… ¿Pero cuáles son esos supuestos?

Supuestos del Método de Mínimos Cuadrados Ordinales

Siempre que se vaya a realizar un modelo de regresión lineal, se debería comprobar que se cumplen estas condiciones. Ahora bien, el hecho de que alguna de las condiciones no se cumpla no significa que el modelo no sea válido, pero sí es algo que afectará al mismo y tendremos que tenerlo en cuenta.

Las condiciones de un modelo de regresión lineal son las siguientes:

No colinealidad o multicolinealidad

El principio de no colinealidad nos viene a decir que si utilizamos varios predictores, estos no deben estar linealmente relacionados entre ellos, ya que sino no se podrá identificar el efecto de cada variable, lo que derivará en una falta de significatividad estadística.

Para detectar si existe colinealidad entre las variables predictoras, podemos realizar distintas técnicas, aunque la más normal es calcular la matriz de correlaciones entre todas las variables predictoras. Si entre dos variables existe una correlación superior a 0.5, deberíamos descartar una de esas variables.

Voy a crear un ejemplo de una nueva variable que sea transformación lineal de otra pero con los datos un poco cambiados. Así veremos el efecto de incluir dos variable altamente correlacionadas en el modelo:

data_colinearity = data[, c("pop","gdpPercap","loggdpPercap","lifeExp")]
data_colinearity$gdpPercap2 = 2*data$gdpPercap
data_colinearity$gdpPercap2[1:20] = data_colinearity$gdpPercap2[1:20]*0.8
cor(data_colinearity)
                     pop  gdpPercap loggdpPercap    lifeExp  gdpPercap2
pop           1.00000000 -0.0556756  -0.02353029 0.04755312 -0.05284217
gdpPercap    -0.05567560  1.0000000   0.87510881 0.67866240  0.99591570
loggdpPercap -0.02353029  0.8751088   1.00000000 0.80898025  0.87020206
lifeExp       0.04755312  0.6786624   0.80898025 1.00000000  0.67732090
gdpPercap2   -0.05284217  0.9959157   0.87020206 0.67732090  1.00000000

Como vemos, existe una correlación importante entre las variables gdpPercap y gdpPercap2 algo normal, ya que, como he explicado, uno es una modificación de la otra variable.

Así pues, veremos como la capacidad predictiva de un modelo con las dos variables es menor que uno con solo una de las dos variables:

lm_colinearity = lm(lifeExp ~  gdpPercap + gdpPercap2 + pop , data = data_colinearity)


lm_non_colinearity = lm(lifeExp ~  gdpPercap + pop, data = data_colinearity)

cat(
  " R2 (adj) no colinealidad: ",summary(lm_non_colinearity)$adj.r.squared,
  "\n",
  "R2 (adj) colinealidad: ",summary(lm_colinearity)$adj.r.squared
) 
 R2 (adj) no colinealidad:  0.4602316 
 R2 (adj) colinealidad:  0.4565031

Efectivamente, el modelo no colineal tiene un mayor R2 ajustado que el modelo con colinealidad, esto es debido al problema de la colinealidad.

Nota: si no sabes qué es el R2 ajustado y por qué se usa, no te preocupes, lo veremos más adelante en este post;)

Relación lineal entre los parámetros

Otra cuestión imprescindible para que la regresión lineal funciones es que la relación entre el predictor y la variable a predecir debe ser lineal. De hecho, en los ejemplos he usado el logaritmo del PIB per Cápita en vez de el PIB per Cápita, precisamente, para evitar este problema. Pero, ¿cómo podemos saber si existe una relación lineal entre el predictor y la variable a predecir?

La forma más sencilla y fácil de conocer si existe una relación es creando un gráfico entre las dos variables. Veamos cómo es la relación entre la esperanza de vida y el PIB per Capita:

plot(data$lifeExp ~  data$gdpPercap)
Relación entre la esperanza de vida y el PIB per Cápita

Como vemos, sí existe una relación entre las dos variables, pero esta no es lineal. En estos casos debemos aplicar transformaciones lineales para conseguir una mayor relación lineal entre las variables.

La transformación a aplicar dependerá del tipo de dato que tengamos. En nuestro caso, la relación parece exponencial, por lo que habrá que aplicar su función inversa: la función logarítmica.

plot(data$lifeExp ~  data$loggdpPercap)
Relación entre la esperanza de vida el logaritmo del PIB per Capita

Como vemos, ahora sí que parece que exista una relación lineal entre ambas variables. De hecho, si comparáramos ambos modelos, aquel que aplica el logaritmo del PIB per Capita tendra un R2 ajustado mucho mayor. Lo comprobamos:


lm_linearity = lm(lifeExp ~  gdpPercap + loggdpPercap, data = data)


lm_non_linearity = lm(lifeExp ~  gdpPercap, data = data)

cat(
  " R2 (adj) sin linealidad: ",summary(lm_non_linearity)$adj.r.squared,
  "\n",
  "R2 (adj) con linealidad: ",summary(lm_linearity)$adj.r.squared
)
 R2 (adj) sin linealidad:  0.4567297 
 R2 (adj) con linealidad:  0.6531915

Principio de Homoscedasticidad

La homoscedasticidad se refiere a que los errores deben tener varianza constante, es decir, que para cada predictor, los errores tienen la misma varianza. Esto se puede ver claramente en este ejemplo del libro Modelos de regresión con R.

Diferencias entre un modelo Homoscedastico y un modelo Heteroscedastico

Para comprobar este supuesto se suele utilizar el test de Breusch-Pagan (aunque también podría usarse el test de White).

El test de Breusch-Pagan comprueba si hay formas lineales de heteroscedasticidad y se puede calcular usando la función bptest del paquete lmtest:

bptest(lm_non_linearity) 

    studentized Breusch-Pagan test

data:  lm_non_linearity
BP = 13.909, df = 1, p-value = 0.0001919

Como vemos, el p<0.05 por lo que hay evidencias de que no se cumple el principio de Homoscedasticidad. Pero… ¿qué implica esto?

Pues, básicamente, cuando no se cumple el principio de Homoscedasticidad, la estimación de la media realizada por el modelo seguirá siendo buena, pero sus intervalos de confianza no lo serán.

Ya conocemos cuáles son las hipótesis en las que se basa nuestra regresión lineal, cuáles se cumplen y cuáles no (y sus implicaciones). Ahora… ¡vamos a ver cómo de bueno es el modelo de regresión lineal que hemos aprendido a programar en R!

Medir la Bondad de Ajuste del modelo de una Regresión Lineal en R

Una vez tenemos el modelo, tenemos que medir cómo de útil es el modelo, ya que de nada nos serviría un modelo que tienen una capacidad predictiva muy baja. Para ello, hay dos métricas principales que se suelen calcular: el error estándar de los residuos y el coeficiente de determinación o R2.

Error Estándar de los Residuos

El error estándar de los residuos mide la desviación de cualquier punto respecto a la estimación del modelo en ese punto. Calcularlo es muy sencillo, ya que únicamente se debe dividir por los grados de libertad de los residuos, como muestra la siguiente fórmula:

 \(RSE = \sqrt{\frac{RSS}{n-p-1}}\)

Nota: p es igual a el número de predictores del modelo. En el caso de un modelo simple, con un predictor.

rse = function(rss, n, p){
  (rss/ (n-p-1))^(1/2)
}

model_rss =  rss(data$lifeExp, data$predicted) 

model_rse = rse(model_rss, nrow(data), 1)
model_rse
[1] 7.122255

El problema (en mi opinión) con el error estándar de los residuos es que interpretarlo no es del todo intuitivo, ya que requiere conocer la escala general de los datos.

Por suerte, existe otra forma muy común de conocer la bondad del modelo que en este caso sí es mucho más intuitiva: el R2.

Coeficiente de determinación R2

El coeficiente de determinación nos indica qué proporción de la varianza total de los datos es capaz de explicar nuestro modelo. A diferencia del RSE, el coeficiente de determinación es adimensional, es decir, es un ratio que va del 0 al 1 para todos los modelos. Esto hace que sea muy fácil de interpretar.

La forma de calcular el coeficiente de determinación es muy sencilla e intuitiva: podemos calcular la varianza total y la varianza de los residuos. Si dividimos la varianza de residuos entre la total, tendríamos el porcentaje de varianza no explicado por el modelo. Lo restante hasta 1, sería la varianza sí explicada por el modelo, es decir, le coeficiente de determinación:

\[R^2 = 1-\frac{RSS}{TSS} = 1 – \frac{\sum(\hat{y_i} – y_i)^2}{\sum(y_i – \bar{y}_i)^2}\]

¡Veámos cómo programar el R2 de nuestra regresión lineal con R!

tss = function(y){
  mean_y = mean(y)
  sum((y - mean_y)^2)
}

r_squared = function(y,y_pred){
  tss_temp = tss(y)
  rss_temp = rss(y, y_pred)
  1- (rss_temp/tss_temp)

}

r_squared = r_squared(data$lifeExp, data$predicted)
r_squared
[1] 0.654449

Como vemos, obtenemos que el modelo es capaz de explicar el 65% de la varianza total, lo cual parece que está bastante bien, teniendo en cuenta lo simple que es el modelo.

Sin embargo, usar el R2 genera un problema: cuantos más predictores incluyamos, más varianza se va a explicar, por lo que mayor será el R2. Esto nos podría llevar a incluir más variables al modelo, cuando esto no es necesariamente bueno. Otro problema es que si usamos el R2 no podríamos comparar modelos con distinto números de predictores, ya que aquellos con más predictores resultarán más beneficiados.

Coeficiente de determinación ajustado: R2ajustado

Para evitar estos problemas, se suele utilizar el coeficiente de determinación ajustado o R2 ajustado, el cual ajusta el R2 con el número de predictores del modelo. De esta forma, podremos comparar modelos con distintos números de predictores y evitaremos tender a incluir más predictores.

\[ R^2_{ajustado} = 1 – \frac{RSS}{TSS}x\frac{n-1}{n-p-1} = R^2 – (1-R^2)x\frac{n-1}{n-p-1} \]

Ahora que ya sabemos cómo ajustar el coeficiente de determinación… ¡vamos a programarlo para nuestra regresión lineal en R!

r_squared_adj = function(y, y_pred){

  n = ifelse(is.null(nrow(y)), length(y), nrow(y))
  p = ifelse(is.null(ncol(y)),1,ncol(y)) 

  tss_temp = tss(y)
  rss_temp = rss(y, y_pred)

  1-((rss_temp/tss_temp) * ((n-1)/(n-p-1)))
}

r_squared_adj(data$lifeExp, data$predicted)
[1] 0.6519808

¡Visto! Como vemos el R2 ajustado es algo inferior al R2 normal, pero tampoco mucho. Seguramente, en una regresión con muchos más predictores las diferencias serían mucho más grandes.

Perfecto, ahora tenemos formas de comparar modelos, pero… ¿cómo podemos saber si esos modelos son significativos? ¿Y los predictores, cómo sabemos si todos son significativos también o si, por el contrario, deberíamos prescindir de alguno de ellos? ¡Veámoslo!

Significancia estadística en la regresión lineal

Significancia del Modelo F-test

Lo primero de todo es comprobar si el modelo es significativo o no, es decir, de si el modelo es capaz de predecir mejor de lo que esperaríamos que acertara por azar. Y es que, las hipótesis del modelo son:

  • H0: β1 = … = βn = 0
  • H0: β1 = … = βn ≠ 0

Para calcular este contraste, se compara la suma de resíduos al cuadrado del modelo con predictores vs el modelo sin predictores. De esta forma, si la diferencia es significativa, es porque al menos uno de los predictores del modelo es significativo. ¡Veámos la fórmula!

\(F = \frac{\frac{TSS – RSS}{p – 1}}{\frac{RSS}{n – p}}\)

Veamos cómo programar esta parte tan importante de la regresión lineal en R!


f_test = function(y, y_pred){

  n = ifelse(is.null(nrow(y)), length(y), nrow(y))
  p = ifelse(is.null(ncol(y)),1,ncol(y)) + 1

  tss_temp = tss(y)
  rss_temp = rss(y, y_pred)

  ((tss_temp - rss_temp) * (n - p))/(rss_temp* (p - 1))  
}

f_value = f_test(data$lifeExp, data$predicted)
f_value
[1] 265.1501

Para saber si este valor es o no significativo, dado el valor de F y los grados de libertad (n-p) y (p-1), podemos obtener un valor de p. Comparando este valor con un treshold, generalmente 0.05, comprobaremos si el valor es o no significativo.

Aunque estemos programando la regresión en R, hacer este paso de forma manual no tiene mucho sentido. Por suerte, R ya cuenta con una función que hace este trabajo por nosotros, así que vamos a ver qué nos devuelve:

pf(f_value, 1, 140, lower.tail = F)
[1] 4.11537e-34

Como vemos, obtenemos un p-value inferior al 0.05. Esto significa que nuestro modelo es significativo, es decir, que nos sirve para predecir la esperanza de vida. Ahora bien, esto no significa que todos los parámetros de nuestro modelo sean significativos: quizás el logaritmo de la esperanza de vida es significativo, pero el intercerpt, no.

¡Veámos cómo analizar la significancia de los parámetros del modelo!

Significancia de los parámetros del modelo

Para comprobar si los parámetros del modelo son o no significativos, para cada uno de los parámetros se debe calcular su significatividad estadística. Para ello se suele utilizar el t-test, que se cacula dividiendo el parámetro por la desviación estándar del mismo:

\(t = \frac{\hatβ_j}{se(\hatβ_j)}\) \(SE(\hatβ_j) = \sqrt\frac{σ^2}{\sum^n_{i=1}(x_{ij}-\bar x)^2}\)

Una vez tenemos el valor de la t, esta se debe buscar en tablas para comprobar el valor de p. Como hemos hecho anteriormente, para ello usaremos la función pt de R. ¡Hagámoslo!

t_test_beta = function(x, estimator_variance, beta){
  n = length(x)
  se = sqrt(model_rse^2 / (var(x) * (n-1)))
  t = beta/se

  return(c(t,se) )

}

t_beta = t_test_beta(data$loggdpPercap, model_rse, slope)
t_beta
[1] 16.2834300  0.4423393

En el caso del intercept, el cálculo es algo diferente:

\(SE(_\hat\alpha) = \sqrt{\frac{S^2_ \hat\beta \; x \; \sum x_i^2}{n}}\)

t_test_alpha = function(x, t_beta, alpha){
  se = (t_beta^2*sum(x^2)/length(x))^0.5
  t = alpha/se
  return(c(t, se))
}

t_alpha = t_test_alpha(data$loggdpPercap, t_beta[2], intercept)
t_alpha
[1] 1.283053 3.857684

Ahora, habría para cada valor de t que hemos recibido, habría que revisarlos en tablas. De esta forma podremos saber si el valor es, o no, significativo. Esto es algo que podemos hacer fácilmente en R con la función pt (ventajas de programar una regresión lineal desde 0 en R;)

# p-value for alpha
2*pt(t_alpha[1],df = (length(data$loggdpPercap)-2), lower.tail = F)

# p-value for beta
2*pt(t_beta[1],df = (length(data$loggdpPercap)-2), lower.tail = F)
[1] 0.2015936
[1] 4.11537e-34

Como vemos, p_alpha>0.05, así que no es un estimador significativo, mientras que p_beta<0.05, por lo que sí es significativo. Lo normal suele ser eliminar del modelo aquellos parámetros que no son significativos, pero ¿ocurre lo mismo con el intercept?

Pues la respuesta es que no necesariamente, ya que puede ser que, a nivel conceptual, en nuestro modelo no tenga sentido que un cuando el predictor sea cero, el valor predicho sea cero.

Por ejemplo, si queremos predecir el precio de una casa según el número de metros cuadrados, tiene sentido que la casa valga 0 si tienen 0m2. Sin embargo, si queremos predecir los ingresos por ventas según el gasto en publicidad, es muy probable que aunque gastemos 0, seguiríamos vendiendo, es decir, el intercept no tiene sentido que sea 0.

En cualquier caso, si tuviéramos más variables en el modelo y estas no resultan significativas, deberíamos excluirlas del mismo hasta quedarnos con variables significativas.

Conclusiones de programar una regresión lineal en R

Y con esto… ya sabes todo lo que hay que saber sobre la regresión lineal: sabes en qué se basa, cómo se calculan sus parámetros, que supuestos siguen, cómo comprobar si el modelo y cada parámetro es significativo y cómo comparar distintos modelos. Todo ello, lo hemos aprendido programando la regresión lineal desde 0 en R.

Como ves, aunque la regresión lineal sea uno de los algoritmos más básicos dentro del mundo del machine-learning, sin duda es un algoritmo que tiene más miga de la que parece y que requiere de ciertos conocimientos estadísticos para su aplicación.

Y es que, aunque en la práctica no vayamos a programar la regresión lineal desde 0 en R, sin duda alguna, entender cómo funciona una regresión lineal, sus aspectos y detalles va a ser clave para poder aplicar la regresión lineal de forma correcta en R o cualquier otro lenguaje.

Y es que, como hemos visto, pequeñas diferencias (como conseguir la linealidad entre variables) pueden suponer grandes cambios en nuestros modelos.

Si te ha gustado conocer algortimos desde 0, te animo a leer también los posts de Cómo programar una red neuronal desde 0 en R, cómo programar el algoritmo K-means desde 0 en R y cómo programar una red neuronal desde 0 en Python.

Como siempre, espero que el post te haya gustado. Si es así, ¡suscríbete a la newsletter para no perderte mis próximos posts! ¡Nos vemos en el siguiente!

¿Te gusta el contenido? ¡Apoya el blog!

Si te gusta el contenido, si quieres puedes apoyar mi blog con una pequeña donación. Así, podré cubrir los costes que supone crear y mantener este blog y podré utilizar más herramientas Cloud con las que seguir creando contenido gratuito para que más personas mejoren como Data Scientist. 

¡No te pierdas ningún post!

Si te gusta lo que lees... suscríbete para estar al día de los contenidos que subo.
¡Serás el primero en enterarte!