Code
library(tidyverse)
library(broom)
library(purrr)
library(plotly)Unidad 3: Determinación de Sistemas de Aprendizaje automático
Regresión Lineal
Unidad 3 (RLin): Determinación de sistemas de aprendizaje automático
Algoritmos de Machine Learning
Regresión Lineal
Regresión Logística
Árboles de Decisión
k Nearest Neighbor
SVM
Naive Bayes
Clustering jerárquico
K-Means
PCA
Redes Neuronales
Aprendizaje profundo
¿Qué es la regresión lineal?
Bagnato (2022)
La regresión lineal es un algoritmo de aprendizaje supervisado que se utiliza en Machine Learning y en estadística. En su versión más sencilla, lo que haremos es “dibujar una recta” que nos indicará la tendencia de un conjunto de datos continuos. En estadística, regresión lineal es una aproximación para modelar la relación entre una variable escalar dependiente “y” y una o más variables explicativas nombradas con “X”.
Conceptos Básicos de regresión
Recordemos rápidamente la fórmula de la recta: Y = mX + b
Donde Y es el resultado, X es la variable, m la pendiente (o coeficiente) de la recta y b la constante o también conocida como el “punto de corte con el eje Y” en la gráfica (cuando x=0).
¿Cómo funciona el algoritmo de regresión lineal en Machine Learning?
Recordemos que los algoritmos de Machine Learning Supervisados, aprenden por sí mismos y -en este caso- a obtener automáticamente esa “recta” que buscamos con la tendencia de predicción. Para hacerlo se mide el error con respecto a los puntos de entrada y el valor “Y” de salida real. El algoritmo deberá minimizar el coste de una función de error cuadrático y esos coeficientes corresponderán con la recta óptima. Hay diversos métodos para conseguir minimizar el coste.
NOTA: cuando hablo de “recta” es en el caso particular de regresión lineal simple. Si hubiera más variables, hay que generalizar el término. Con 2 variables predictoras sería la ecuación de un plano en el espacio.
Métricas:
- Error cuadrático Medio ECM o MSE (en inglés):
\[\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\] - Para pasarlo a unidades (sin estar elevadas al cuadrado), se suele usar RMSE:
\[\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}\]
Ejemplo Regresión lineal simple: Temperatura ambiente en función del canto de los grillos
Se va a trabajar con datos estudiados por biólogos en los que se analiza la relación entre la frecuencia del canto de los grillos y la temperatura del aire. Usaremos el canto de los grillos como “termómetro”. Por lo tanto, nuestra variable objetivo va a ser la temperatura.
# Leemos el fichero: viene en formato .rds (formato de r comprimido)
grillos <- read_rds("https://raw.githubusercontent.com/jesusturpin/curintel2324/main/data/grillos.rds")# Examinamos el dataset:
glimpse(grillos)Rows: 170
Columns: 3
$ freq <dbl> 32, 85, 142, 110, 92, 36, 28, 135, 78, 18, 98, 24, 147, 100, 2…
$ temp <dbl> 9.399821, 17.160850, 26.720099, 19.462094, 16.282250, 10.96216…
$ especie <chr> "Gryllus campestris", "Gryllus campestris", "Gryllus campestri…# Visualizamos el resumen de los datos
summary(grillos) freq temp especie
Min. : -1.00 Min. :-0.8629 Length:170
1st Qu.: 35.25 1st Qu.: 9.1866 Class :character
Median : 63.50 Median :13.5386 Mode :character
Mean : 67.25 Mean :14.3767
3rd Qu.: 93.75 3rd Qu.:18.1872
Max. :151.00 Max. :41.4029 # La variable especie es categórica, pero está en formato chr. La convertimos a factor, de este modo, averiguamos cuántas especies hay y cuáles son:
grillos$especie <- factor(grillos$especie)# Visualizamos de nuevo el resumen de los datos
summary(grillos) freq temp especie
Min. : -1.00 Min. :-0.8629 Ensifera :62
1st Qu.: 35.25 1st Qu.: 9.1866 Gryllus campestris:59
Median : 63.50 Median :13.5386 Oecanthus fultoni :49
Mean : 67.25 Mean :14.3767
3rd Qu.: 93.75 3rd Qu.:18.1872
Max. :151.00 Max. :41.4029 # Buscamos posibles errores en los datos: hay una frecuencia negativa!!
grillos %>%
filter(freq <= 0) freq temp especie
1 -1 -0.8629387 Oecanthus fultoni# Eliminamos el error
grillos <- grillos %>%
filter(freq > 0)Visualizamos la relación temperatura vs frecuencia:
Code
grillos %>%
ggplot(aes(freq, temp))+
geom_point()+
geom_smooth(method = "lm", se = FALSE)+
theme_bw()
# calculamos la correlación lineal:
cor(grillos$freq, grillos$temp)[1] 0.6812606#Creamos un modelo de regresión lineal simple que sirva para predecir la temperatura en función de la frecuencia.
mdl_temp_vs_freq <- lm(data = grillos, formula = temp ~ freq)
mdl_temp_vs_freq
Call:
lm(formula = temp ~ freq, data = grillos)
Coefficients:
(Intercept) freq
4.4639 0.1478 \(R^2 = 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2}\)
# Evaluamos el modelo: Coeficiente de determinación (R^2) (calculado)
1 - (sum((residuals(mdl_temp_vs_freq))**2) / (sum((grillos$temp - mean(grillos$temp))**2)))[1] 0.464116# Resumen del modelo (fórmula, coeficientes y métricas más importantes)
summary(mdl_temp_vs_freq)
Call:
lm(formula = temp ~ freq, data = grillos)
Residuals:
Min 1Q Median 3Q Max
-13.7663 -5.7758 -0.3755 5.1364 14.9096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.46392 0.94991 4.699 5.43e-06 ***
freq 0.14785 0.01229 12.026 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.965 on 167 degrees of freedom
Multiple R-squared: 0.4641, Adjusted R-squared: 0.4609
F-statistic: 144.6 on 1 and 167 DF, p-value: < 2.2e-16# Resumen "ordenado": métricas
broom::glance(mdl_temp_vs_freq)# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.464 0.461 5.96 145. 2.14e-24 1 -541. 1087. 1097.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int># Resumen "ordenado": coeficientes
broom::tidy(mdl_temp_vs_freq)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.46 0.950 4.70 5.43e- 6
2 freq 0.148 0.0123 12.0 2.14e-24Predicción 100 cantos/s
Según el modelo, si tenemos una frecuencia de 100 cantos por segundo las 9 de la mañana. ¿Cuál es la temperatura?
# Obtenemos los coeficientes del modelo:
coeffs <- coefficients(mdl_temp_vs_freq)
coeffs(Intercept) freq
4.4639245 0.1478483 b <- coeffs[1] # ordenada en el origen (Intercept)
m <- coeffs[2] # pendiente (freq)
freq_9am <- 100t_9am <- b + m*freq_9am # Y = b + mx
t_9am(Intercept)
19.24876 Predicción 200 cantos/s
Con una temperatura de 19.25ºC a las 9am, si la frecuencia se duplica a las 3 de la tarde, ¿Cuántos grados habría aumentado la temperatura? Calcula el error cuadrático medio. ¿Según el error cuadrático medio RMSE, con qué error afirmas los resultados?
# A las 3 de la tarde: 200 cantos por segundo
freq_3pm <- 200
t_3pm <- b + m*freq_3pm # Y = b + mx
t_3pm(Intercept)
34.03359 # Tenemos 34.03 º a las 3 pm. Calculamos el RMSE
# Error cuadrático medio: MSE --> Raiz MSE RMSE
MSE <- sum(residuals(mdl_temp_vs_freq)**2)/(nrow(grillos))
RMSE <- sqrt(MSE)
RMSE # grados[1] 5.929384# Aumento de la temperatura en grados
t_3pm - t_9am(Intercept)
14.78483 Predicciones automáticas:
test_grillos <- expand.grid(
freq = c(100,200)
)test_grillos %>%
mutate(temp = predict(mdl_temp_vs_freq,
newdata = select(., freq))) #select(test_grillos, freq) freq temp
1 100 19.24876
2 200 34.03359Regresión Lineal Múltiple
Modelo una variable explicativa categórica
mdl_temp_vs_esp <- lm(temp ~ especie, grillos)
mdl_temp_vs_esp
Call:
lm(formula = temp ~ especie, data = grillos)
Coefficients:
(Intercept) especieGryllus campestris
20.598 -6.566
especieOecanthus fultoni
-13.517 # Añadir + 0, modifica el nombre del coeficiente, pero no cambia el modelo
mdl_temp_vs_freq_esp_dummy <- lm(temp ~ especie + 0, grillos)
mdl_temp_vs_freq_esp_dummy
Call:
lm(formula = temp ~ especie + 0, data = grillos)
Coefficients:
especieEnsifera especieGryllus campestris
20.598 14.032
especieOecanthus fultoni
7.081 Code
# Crear un nuevo data frame con las predicciones
predicciones <- data.frame(especie = unique(grillos$especie),
Predicted = predict(mdl_temp_vs_esp, newdata = data.frame(especie = unique(grillos$especie))))
# Datos originales y las predicciones
ggplot(grillos, aes(x = especie, y = temp)) +
geom_jitter(width = 0.2, alpha = 0.5) + # Usamos geom_jitter para evitar superposiciones
geom_point(data = predicciones, aes(x = especie, y = Predicted), color = "red", size = 4, shape = "square") +
labs(title = "Predicción de la temperatura ambiente por Especie de grillo", x = "Especie", y = "Temperatura") +
theme_bw()
summary(mdl_temp_vs_freq_esp_dummy)
Call:
lm(formula = temp ~ especie + 0, data = grillos)
Residuals:
Min 1Q Median 3Q Max
-11.4131 -4.7005 -0.8878 3.9985 20.8045
Coefficients:
Estimate Std. Error t value Pr(>|t|)
especieEnsifera 20.5984 0.7716 26.696 < 2e-16 ***
especieGryllus campestris 14.0321 0.7910 17.741 < 2e-16 ***
especieOecanthus fultoni 7.0814 0.8769 8.075 1.31e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.075 on 166 degrees of freedom
Multiple R-squared: 0.8681, Adjusted R-squared: 0.8657
F-statistic: 364.2 on 3 and 166 DF, p-value: < 2.2e-16grillos %>%
group_by(especie) %>%
summarise("coeficientes(media)" = mean(temp))# A tibble: 3 × 2
especie `coeficientes(media)`
<fct> <dbl>
1 Ensifera 20.6
2 Gryllus campestris 14.0
3 Oecanthus fultoni 7.08Rectas paralelas: Variable categórica + numérica
2 variables predictoras(numérica+categórica): Pendiente común + 3 ordenadas en el origen
mdl_temp_vs_freq_esp <- lm(temp ~ especie + freq, grillos)
coefficients(mdl_temp_vs_freq_esp) (Intercept) especieGryllus campestris especieOecanthus fultoni
10.744776 -7.511183 -13.719990
freq
0.151369 Simplifica los coeficientes:
mdl_temp_vs_freq_esp <- lm(temp ~ especie + freq + 0, grillos)
coefficients(mdl_temp_vs_freq_esp) especieEnsifera especieGryllus campestris especieOecanthus fultoni
10.744776 3.233593 -2.975214
freq
0.151369 Code
grillos %>%
ggplot(aes(freq, temp, color = especie)) +
geom_point() +
moderndive::geom_parallel_slopes(se = FALSE)+
theme_bw()
summary(mdl_temp_vs_freq_esp)
Call:
lm(formula = temp ~ especie + freq + 0, data = grillos)
Residuals:
Min 1Q Median 3Q Max
-7.7070 -1.1981 0.1876 1.2444 8.1042
Coefficients:
Estimate Std. Error t value Pr(>|t|)
especieEnsifera 10.744776 0.397746 27.014 < 2e-16 ***
especieGryllus campestris 3.233593 0.422853 7.647 1.61e-12 ***
especieOecanthus fultoni -2.975214 0.428283 -6.947 8.21e-11 ***
freq 0.151369 0.004443 34.069 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.15 on 165 degrees of freedom
Multiple R-squared: 0.9836, Adjusted R-squared: 0.9832
F-statistic: 2472 on 4 and 165 DF, p-value: < 2.2e-16Ejemplo predicciones con Regresión lineal múltiple: 100 y 200 cantos/s para las 3 especies
Ahora la especie también es predictora:
test_grillos <- expand.grid(
freq = c(100,200),
especie = unique(grillos$especie)
)
test_grillos %>%
mutate(temp = predict(mdl_temp_vs_freq_esp,
newdata = .)) # el . indica todas las predictoras freq especie temp
1 100 Gryllus campestris 18.37049
2 200 Gryllus campestris 33.50739
3 100 Oecanthus fultoni 12.16169
4 200 Oecanthus fultoni 27.29859
5 100 Ensifera 25.88168
6 200 Ensifera 41.01858Code
glance(mdl_temp_vs_freq_esp)# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.984 0.983 2.15 2472. 4.71e-146 4 -367. 744. 760.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Precisión del modelo múltiple
¿Cuál es el error cuadrático medio en ºC del nuevo modelo?
Regresión Lineal Múltiple: Interacciones
El modelo anterior, la variable especie no condiciona la pendiente de la recta. No hay interacciones:
\[Y = b_1 D_1 + b_2 D_2 + \ldots + b_k D_k + \beta \times \text{freq} + \epsilon\]
\[\hat{Y} = b_1 D_1 + b_2 D_2 + \ldots + b_k D_k + \beta \times \text{freq}\]
Ejercicio 1:
Identifica todos los elementos en la ecuación anterior. Vamos a usar el modelo para predecir con un grillo Ensifera y otro campestris la temperatura a 50 cantos/s
Code
test_grillos <- expand.grid(
freq = c(50),
especie = c("Ensifera", "Gryllus campestris")
)
test_grillos %>%
mutate(temp = predict(mdl_temp_vs_freq_esp,
newdata = .)) # el . indica todas las predictoras freq especie temp
1 50 Ensifera 18.31323
2 50 Gryllus campestris 10.80204Code
# Ensifera
b1 + beta*freqespecieEnsifera
18.31323 Code
# Campestris
b2 + beta*freqespecieGryllus campestris
10.80204 - Y –> Es el valor real, no el del modelo
- \(\hat{Y}\) –> Es el valor que calcula el modelo
- b1
- D1
- b2
- D2
- b3
- D3
- \(\beta\)
- freq
- \(\epsilon\) –> El error residual (variable aleatoria)
Modelo con interacciones: mdl_temp_vs_freq_esp_inter
En un modelo de regresión lineal en R, las interacciones entre variables se especifican utilizando el operador : o * en la fórmula del modelo. Estos operadores permiten modelar cómo el efecto de una variable predictora sobre la variable respuesta puede cambiar en función de los niveles de otra variable predictora.
mdl_temp_vs_freq_esp_inter <- lm(temp ~ especie + freq + especie:freq + 0, grillos)
# Equivalente: lm(temp ~ especie * freq + 0)coeffs <- coefficients(mdl_temp_vs_freq_esp_inter)
coeffs especieEnsifera especieGryllus campestris
7.91713697 5.04735869
especieOecanthus fultoni freq
-0.66246478 0.19480648
especieGryllus campestris:freq especieOecanthus fultoni:freq
-0.06886208 -0.07824837 La fórmula con interacciones se complica, apareciendo nuevos coeficientes:
\[\hat{Y} = b_1 D_1 + b_2 D_2 + \ldots + b_k D_k + \beta \times \text{freq} + \gamma_1 \times D_1 \times \text{freq} + \gamma_2 \times D_2 \times \text{freq} + \ldots + \gamma_k \times D_k \times \text{freq}\] Para esto se traduce visualmente en que las pendientes dependen de la especie
Code
grillos %>%
ggplot(aes(freq, temp, color = especie)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme_bw()
Ejercicio 2:
Identifica todos los elementos en la ecuación anterior (con interacciones). Vamos a usar el modelo para predecir con un grillo de cada especie la temperatura a 50 cantos/s. Verifica con predict y usando la ecuación que las predicciones son iguales.
freq especie temp
1 50 Gryllus campestris 11.344579
2 50 Oecanthus fultoni 5.165441
3 50 Ensifera 17.657461especieEnsifera
17.65746 especieGryllus campestris
11.34458 especieOecanthus fultoni
5.165441 Ejercicio 3:
Sustituye el modelo mdl_temp_vs_freq_esp_inter por tres modelos simples, uno por cada especie. Para ello, deberás separar el dataset grillos en 3 dataset.
Indica los coeficientes de cada modelo y realiza las predicciones con predict y operando con los coeficientes.
Utilizando un solo panel, representa los 3 modelos con todos sus datos y predicciones.
2 variables numéricas ejemplo RL con y sin interacciones
Añadimos datos sintéticos: variable medida, teniendo en cuenta las medidas según las fuentes consultadas sobre tamaño medio de las especies. Para simplificar, no se tiene en cuenta el sexo.
set.seed(0)
grillos <- grillos %>%
mutate(medida = map_dbl(especie, ~case_when(
.x == "Gryllus campestris" ~ rnorm(1, mean = 22.5, sd = 6),
.x == "Oecanthus fultoni" ~ rnorm(1, mean = 16.5, sd = 6),
.x == "Ensifera" ~ rnorm(1, mean = 45, sd = 15),
TRUE ~ NA_real_)))glimpse(grillos)Rows: 169
Columns: 4
$ freq <dbl> 32, 85, 142, 110, 92, 36, 28, 135, 78, 18, 98, 24, 147, 100, 2…
$ temp <dbl> 9.399821, 17.160850, 26.720099, 19.462094, 16.282250, 10.96216…
$ especie <fct> Gryllus campestris, Gryllus campestris, Gryllus campestris, Gr…
$ medida <dbl> 30.07773, 30.13458, 16.92860, 36.92792, 15.61406, 20.03094, 25…Code
grillos %>%
ggplot() +
geom_point(aes(medida, temp, color = especie)) +
geom_smooth(aes(medida, temp, color = especie), method = "lm", se = FALSE)+
theme_bw()
Modelo con 2 predictoras numéricas en 3D (sin interacciones)
mdl_temp_vs_freq_med <- lm(temp ~ freq + medida, data = grillos)
coefficients(mdl_temp_vs_freq_med)(Intercept) freq medida
-3.0126707 0.1494446 0.2646969 Ecuación y coeficientes:
\[ z = \beta_0 + \beta_1x + \beta_2y + \varepsilon\]
Modelo con 2 predictoras numéricas en 3D (con interacciones)
mdl_temp_vs_freq_med_inter <- lm(temp ~ freq * medida, data = grillos)
coefficients(mdl_temp_vs_freq_med_inter) (Intercept) freq medida freq:medida
-1.1189767577 0.1232495070 0.1986752264 0.0009154477 Ecuación y coeficientes:
\[ z = \beta_0 + \beta_1x + \beta_2y + \beta_3xy + \varepsilon\]
References
Bagnato, Juan Ignacio. 2022. Aprende Machine Learning En Español: Teoría + Práctica Python. Leanpub. http://leanpub.com/aprendeml.