Segundo ejemplo de construcción de modelos estadísticos utilizando regresión lineal simple
En este ejemplo vamos a entender la relación que existe entre diferentes medidas geométricas de unos árboles, utilizando la base de datos embebida en R
Medidas de árboles
En este ejemplo trabajaremos con el set de datos trees, que contiene datos sobre la circunferencia (en pulgadas), altura (en pies) y volumen (en pies cúbicos) del tronco de árboles de cerezos. Intentaremos ajustar un modelo de regresión lineal simple para predecir el volumen en función del diámetro.
Análisis exploratorio de datos (EDA)
• head(datos) -> Inspeccionar primeras observaciones del set de datos
• glimpse(datos) -> Estructura del set de datos (parecido a str())
• summary(datos) -> Resumen del set de datos (datos numéricos de las variables)
• pairs(datos) -> Gráficos dispersión del conjunto de variables
• cor.test() -> Test de correlación (“pearson”, “kendall”, “spearman”)
• ggpairs() -> Combina en un único gráfico diagramas de dispersión, distribución de las variables y los valores de correlación.
• ggplot()
• plot()
Consideraciones previas
Antes de generar el modelo, representaremos los datos para cuantificar la posible relación lineal entre variables y su magnitud. Si no detectáramos esta relación pasaríamos a plantearnos métodos de ajuste alternativo
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(xfun)##
## Attaching package: 'xfun'## The following objects are masked from 'package:base':
##
## attr, isFALSEhead(trees)## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7Conociendo los datos
glimpse(trees)## Rows: 31
## Columns: 3
## $ Girth <dbl> 8.3, 8.6, 8.8, 10.5, 10.7, 10.8, 11.0, 11.0, 11.1, 11.2, 11.3, …
## $ Height <dbl> 70, 65, 63, 72, 81, 83, 66, 75, 80, 75, 79, 76, 76, 69, 75, 74,…
## $ Volume <dbl> 10.3, 10.3, 10.2, 16.4, 18.8, 19.7, 15.6, 18.2, 22.6, 19.9, 24.…Resumen estadístico simple
summary(trees)## Girth Height Volume
## Min. : 8.30 Min. :63 Min. :10.20
## 1st Qu.:11.05 1st Qu.:72 1st Qu.:19.40
## Median :12.90 Median :76 Median :24.20
## Mean :13.25 Mean :76 Mean :30.17
## 3rd Qu.:15.25 3rd Qu.:80 3rd Qu.:37.30
## Max. :20.60 Max. :87 Max. :77.00Matriz de diagramas de dispersión
pairs(trees)
Visualizando un diagrama de dispersión
library(ggplot2)
ggplot(data = trees, mapping = aes(x = Girth, y = Volume)) +
geom_point(color = "firebrick", size = 2) +
labs(title = "Diagrama de dispersión", x = "Diámetro", y = "Volumen") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
## Coeficiente de correlación de pearson
cor.test(x = trees$Girth, y = trees$Volume, method = "pearson", digits = 3)##
## Pearson's product-moment correlation
##
## data: trees$Girth and trees$Volume
## t = 20.478, df = 29, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9322519 0.9841887
## sample estimates:
## cor
## 0.9671194Visualizando gráficamente la relación entre las variables
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(trees, lower = list(continuous = "smooth"), diag = list(continuous = "bar"), axisLabels = "none")## Warning in check_and_set_ggpairs_defaults("diag", diag, continuous =
## "densityDiag", : Changing diag$continuous from 'bar' to 'barDiag'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
De lo hasta ahora analizado podemos concluir que:
Observando los gráficos de dispersión podemos observar que la variable Girth (diámetro) está más linealmente asociada con la variable respuesta Volume, por lo que utilizaremos ésta para el modelo.
El coeficiente de correlación de Pearson es bastante alto (r = 0.967) y significativo (p-value = 2,2x10−16). Ello indica una correlación entre ambas variables bastante intensa.
Tiene sentido generar el modelo de regresión lineal cumpliéndose estos primeros requisitos.
2. Cálculo del modelo de regresión lineal simple
• lm(y ~ x, data) -> Modelo de regresión lineal simple
• summary(modelo) -> Resumen del ajuste del modelo
modelo.lineal <- lm(Volume ~ Girth, data = trees)
# Información del modelo
summary(modelo.lineal)##
## Call:
## lm(formula = Volume ~ Girth, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.065 -3.107 0.152 3.495 9.587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.9435 3.3651 -10.98 7.62e-12 ***
## Girth 5.0659 0.2474 20.48 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9331
## F-statistic: 419.4 on 1 and 29 DF, p-value: < 2.2e-16El summary del modelo contiene los errores estándar, el valor del estadístico t y el correspondiente p-value de ambos parámetros β0̂ y β1̂ . El p-value nos permite determinar si los estimadores de los parámetros son significativamente distintos de 0, es decir, que contribuyen al modelo. El parámetro que suele ser más útil en estos modelos es la pendiente.
Conclusiones del modelo:
Tanto la ordenada en el origen como la pendiente son significativas (p-value = 2,2x10−16 ).
El coeficiente de determinación R2indica que el modelo es capaz de explicar el 93% de la variabilidad presente en la variable respuesta (volumen) mediante la variable independiente (diámetro).
El p-value obtenido en el test F (2,2x10−16) determina que sí es significativamente superior la varianza explicada por el modelo en comparación con la varianza total, por lo que podemos aceptar nuestro modelo como válido y útil.
Ecuación del modelo: Volume = - 36.94 + 5.065 Girth -> por cada unidad que se incrementa el diámetro, el volumen aumenta en promedio 5,065 unidades.
3. Intervalos de confianza para los parámetros del modelo
• confint(modelo) -> Intervalos de confianza para los coeficientes del modelo
confint(modelo.lineal)## 2.5 % 97.5 %
## (Intercept) -43.825953 -30.060965
## Girth 4.559914 5.5717994. Representación gráfica del modelo
Representamos la línea de mínimos cuadrados, representaremos el intervalo de confianza (límites superior e inferior) para cada predicción con la función geom_smooth( ) del paquete ggplot2. Esto permite identificar la región en la que se encuentra el promedio de la variable respuesta según el modelo generado y para un determinado nivel de confianza:
ggplot(data = trees, mapping = aes(x = Girth, y = Volume)) +
geom_point(color = "firebrick", size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
labs(title = "Volumen ~ Diámetro", x = "Diámetro", y = "Volumen") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) ## `geom_smooth()` using formula 'y ~ x'
5. Verificar condiciones para aceptar el modelo
• plot(modelo) -> Análisis de los residuos (distribución, variabilidad…)
• shapiro.test(modelo$residuals) -> Test de hipótesis de Shapiro Wilk para el análisis de normalidad
• bptest(modelo) -> Test de contraste de homocedasticidad Breusch-Pagan
• influence.measures(modelo) -> Detección de observaciones influyentes
• influencePlot(modelo) -> Visualización de observaciones influyentes
• outlierTest(modelo) -> Test de detección de outliers
• rstudent(modelo) -> Cálculo de residuos estudentizados
Para evaluar las condiciones que permiten dar como válido el modelo lineal, haremos uso principalmente del análisis de los residuos:
par(mfrow=c(1,2))
plot(modelo.lineal)
## Análisis de normalidad de los residuos
# Contraste de hipótesis (normalidad de los residuos)
shapiro.test(modelo.lineal$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo.lineal$residuals
## W = 0.97889, p-value = 0.7811Autocorrelación de los residuos
# Análisis gráfico autocorrelación de los residuos
ggplot(data = trees, aes(x = seq_along(modelo.lineal$residuals),
y = modelo.lineal$residuals)) +
geom_point(aes(color = modelo.lineal$residuals)) +
scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
geom_line(size = 0.3) +
labs(title = "Distribución de los residuos", x = "index", y = "residuo")+
geom_hline(yintercept = 0) +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
En estadística se dice que un modelo predictivo presenta homocedasticidad cuando la varianza del error condicional a las variables explicativas es constante a lo largo de las observaciones
En este caso, la normalidad de los residuos podemos aceptarla, y tampoco parecen seguir una clara tendencia según el orden de registro de las observaciones, pero la condición de homocedasticidad parece no cumplirse. Al observar los gráficos podríamos sospechar que la observación 31 podría estar influyendo al modelo. Para analizar en qué medida pueda estar influyendo esta u otras observaciones, reajustaremos el modelo excluyendo posibles observaciones sospechosas. Dependiendo de la finalidad del modelo, la exclusión de posibles outliers debe analizarse con detalles, ya que estas observaciones podrían ser errores de medida, pero también podrían representar casos interesantes.
# Residuos estudentizados
studentized_residual <- rstudent(modelo.lineal)
which(abs(studentized_residual) > 3)## named integer(0)Que medidas tienen una mayor influencia sobre el modelo?
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodesummary(influence.measures(model = modelo.lineal))## Potentially influential observations of
## lm(formula = Volume ~ Girth, data = trees) :
##
## dfb.1_ dfb.Grth dffit cov.r cook.d hat
## 31 -1.20_* 1.37_* 1.49_* 0.82 0.89_* 0.22_*Ahora, vamos a ver los residuos estudentizados de forma gráfica
influencePlot(model = modelo.lineal)
## StudRes Hat CookD
## 1 1.315573 0.11514037 0.10983618
## 20 -2.030440 0.03328798 0.06408063
## 31 2.837732 0.21519431 0.88805814El análisis de los residuos estudentizados no detecta ninguna observación atípica. Sin embargo, tal y como sospechábamos, la observación 31 parece estar influenciando al modelo.
Procederemos a reajustar el modelo excluyendo la observación 31:
ggplot(data = trees, mapping = aes(x = Girth, y = Volume)) +
geom_point(color = "grey50", size = 2) +
#recta de regresión con todas las observaciones
geom_smooth(method = "lm", se = FALSE, color = "black") +
#se resalta el valor excluido
geom_point(data = trees[31, ], color = "red", size = 2) +
#se añade la nueva recta de regresión
geom_smooth(data = trees[-31, ], method="lm", se =FALSE,color = "blue") +
labs(x = "Diámetro", y = "Volumen") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5), legend.position = "none")## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Ahora, elaborando de nueva cuenta el modelo
modelo.lineal2 <- lm(Volume ~ Girth, data = trees[-31,])
summary(modelo.lineal2)##
## Call:
## lm(formula = Volume ~ Girth, data = trees[-31, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5036 -2.3834 -0.0489 2.3951 6.3726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.3104 3.2784 -10.16 6.76e-11 ***
## Girth 4.7619 0.2464 19.33 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.813 on 28 degrees of freedom
## Multiple R-squared: 0.9303, Adjusted R-squared: 0.9278
## F-statistic: 373.6 on 1 and 28 DF, p-value: < 2.2e-16La eliminación de la observación identificada como influyente a penas cambia la recta de mínimos cuadrados.
6. Uso del modelo para predecir nuevas observaciones
• predict() -> Predicciones de nuevas observaciones a partir del modelo
Dando por válido el modelo podemos utilizar la función genérica predict() para predecir el valor de nuevas observaciones. En nuestro ejemplo, podemos predecir el volumen a partir de mediciones del diámetro de nuevos árboles:
# Volumen PROMEDIO que esperaríamos de árboles de 15 pulgadas
predict(modelo.lineal, data.frame(Girth = 15), interval = "confidence")## fit lwr upr
## 1 39.04439 37.24858 40.84019Como es de esperar, el intervalo de predicción es mayor que el de confianza. Podemos interpretar la predicción de la siguiente manera: se espera que en promedio un árbol de diámetro 15 tenga un volumen de 39,04. Podemos afirmar con un 95% de confianza que el verdadero valor promedio se encuentra entre (37,24 – 40.84), mientras que el intervalo de predicción para un solo árbol de este diámetro será de (30,16 – 47, 92).