La correlación lieal y la regresión lineal simple son métodos estadísticos que estudian la relación liean existente entre dos variables.
En la descripción de la ANOVA, vimos como los minimos cuadrados son implementaods dentro la función lm()
para ajustar un modelo con variables categoricas (factores). A diferencia de la ANOVA y la T de Student, el objetivo de la regresión no es probar una hipótesis sino estimar una relación que pueda ser usado para predicción.
La construcción de un modelo nos permite entender las relaciones entre las variables y hacer predicciónes usando futuras observaciones.
Un modelo lineal general se utiliza para predecir el valor de una variable continua (variable de respuesta) a partir de una o ma svariables explicativas. El modelo lieal general toma la forma:
\(Y = \beta_0 + \beta_1X_1 + \beta_2X_2+...+\beta_nX_n + \epsilon\)
donde \(y\) es la variable de respuesta, \(X_i\) son las variables explicativa, \(\beta_i\) son los coeficientes y \(\epsilon\) representa el error aleatorio.
Las variables explicativas pueden ser continuas o categorias.
Se asume que el error aleatorio es independiente, que sigue una distribución normal con una media igual a cero, y que tiene la misma varianza para todas las variables explicativas.
Como se ha mencionado, la regresión lieanl simple consiste en generar un modelo de regresión (ecuación de una recta) que permita explicar la relación lineal entre dos variables. A la variable dependiente o de respuesta se le identifica como y y a la variable predictora o independiente como x
El modelo de regresión lieal simple se describe de acuerdo a la ecuación:
\(y = a + bx\)
donde x es la variable explicativa, y es la variable de respuesta (predicción), a es el intercepto y b es la pendiente. Tambien hay una variación residual (error) que no se incluye en esta formula.
O alternativamente:
\(Y = \beta_0 + \beta_1X_1 + \epsilon\)
Para cada uno de los parámetros de la ecuación de regresión lineal simple se puede calcular su significancia y su intervalo de confianza. El test estadistio más empleado es el T-test.
EL test de significancia para la pendiente (\(b\)) del modelo lineal considera como hipótesis:
De esta misma forma tambien se aplica a \(a\).
Nuevamente utilizaremos el ejemplo dentro del libro de Andy Hector con datos de densidad y dureza en muestras de madera. Esto con el objetivo de predecir la dureza (desconocida) de muestras de madera basado en la medición de la densidad.
library(tidyverse)
janka <- read_csv("Data/janka.csv")
## Parsed with column specification:
## cols(
## Densidad = col_double(),
## Dureza = col_double(),
## Otra_variable = col_double()
## )
El primer paso antes de generar un modelo de regresión es representar los datos para poder intuir si existe una relació. Si en este paso no se detecta esta posible relacion, no tiene sentido hacer un modelo lieal.
ggplot(janka, aes(x = Densidad, y = Dureza))+
geom_point()+
labs(titile = "diagrama de dispersión",
y = "Dureza (lbf)", x = "Densidad (lbs pie-3)")
Ahora podemos construir el modelo lieanl con la funcion:
janka_lm <- lm(Dureza ~ Densidad, data = janka)
janka_lm
##
## Call:
## lm(formula = Dureza ~ Densidad, data = janka)
##
## Coefficients:
## (Intercept) Densidad
## -1160.50 57.51
En este modelo, no es necesario especificar el intercepto dentro del modelo ya que R lo incluye automaticamente. Para construir un modelo sin el intercepto se puede usar:
lm(Dureza ~ -1+Densidad, data = janka)
##
## Call:
## lm(formula = Dureza ~ -1 + Densidad, data = janka)
##
## Coefficients:
## Densidad
## 34.13
Si recordamos la formula \(y = a + bx\), la pendiente (\(b\)) corresponde al combio predicho en la dureza por cada unidad de cambio en la densidad.
El intercepto de la regresión es el valor de \(y\) (Dureza) predicho por la linea de regresión para un valor de \(x\) (densidad) de cero.
En este caso, el intercepto da un valor de dureza negativo. Esto se debe a que se esta extrapolando el valor a un valor de densidad de cero.
De esta manera, con los resultados del modelo podemos decir que:
Dureza = -1160.50 + 57.51 x Densidad
Para ver la validez del modelo podemos usar la función summary()
summary(janka_lm)
##
## Call:
## lm(formula = Dureza ~ Densidad, data = janka)
##
## Residuals:
## Min 1Q Median 3Q Max
## -338.40 -96.98 -15.71 92.71 625.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1160.500 108.580 -10.69 2.07e-12 ***
## Densidad 57.507 2.279 25.24 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 183.1 on 34 degrees of freedom
## Multiple R-squared: 0.9493, Adjusted R-squared: 0.9478
## F-statistic: 637 on 1 and 34 DF, p-value: < 2.2e-16
La primera columna (Estimate) devuelve el valor estimado para los dos parámetros de la ecuación del modelo lineal (\(a\) y \(b\)) que equivalen a la ordenada en el origen y la pendiente.
Se muestran los errores estándar, el valor del estadístico t y el p-value (dos colas) de cada uno de los dos parámetros. Esto permite determinar si los parámetros son significativamente distintos de 0.
Para el modelo generado, tanto la ordenada en el origen (intercepto) como la pendiente son significativas (p-values < 0.05).
El valor de \(R^2\) indica que el modelo calculado explica el 94% de la variabilidad presente en la variable respuesta mediante la variable independiente.
El p-value obtenido en el test F determina que sí es significativamente superior la varianza explicada por el modelo en comparación a la varianza total. Es el parámetro que determina si el modelo es significativo y por lo tanto se puede aceptar.
Podemos extraer los valores ajustados del modelo para cadavalor observado con la función fitted()
fit <- fitted(janka_lm)
fit
## 1 2 3 4 5 6 7 8
## 259.9152 265.6658 409.4325 472.6899 472.6899 507.1939 581.9525 719.9686
## 9 10 11 12 13 14 15 16
## 886.7379 1053.5073 1070.7593 1099.5126 1105.2633 1134.0166 1157.0193 1174.2713
## 17 18 19 20 21 22 23 24
## 1180.0220 1180.0220 1306.5366 1473.3060 1536.5633 1611.3220 1801.0940 1801.0940
## 25 26 27 28 29 30 31 32
## 1910.3567 2059.8741 2088.6274 2134.6328 2151.8848 2243.8954 2278.3994 2634.9408
## 33 34 35 36
## 2715.4502 2795.9595 2813.2115 2813.2115
y podemos obtener los residuales (diferencia entre el valor observado - el esperado) con la funcion residuals()
residuales <- residuals(janka_lm)
residuales
## 1 2 3 4 5 6
## 224.0848370 161.3341695 3.5674826 44.3101404 76.3101404 140.8061355
## 7 8 9 10 11 12
## 5.0474583 -15.9685611 92.2620821 -139.5072748 -0.7592772 -79.5126146
## 13 14 15 16 17 18
## 104.7367180 -145.0166194 2.9807107 -164.2712918 -80.0219592 -50.0219592
## 19 20 21 22 23 24
## -36.5366437 -293.3060005 -136.5633428 148.6779800 -91.0940467 208.9059533
## 25 26 27 28 29 30
## -30.3567287 -79.8740831 -268.6274205 -114.6327603 -171.8847628 66.1045576
## 31 32 33 34 35 36
## -338.3994472 625.0591692 -15.4501754 94.0404799 -73.2115225 326.7884775
Alternativamenete, Para poder representar el intervalo de confianza a lo largo de todo el modelo se recurre a la función predict()
. COn esto es posible generar una grafica de regresión con los limites superiores e inferiores del 95% IC.
Para esto, vamos a generar una secuencia con 100 valores
xseq <- seq( min(janka$Densidad), max(janka$Densidad), length= 100 )
Ahora podemos sustituir esta secuencia con los valores de densidad observados
prediccion <- predict(janka_lm, list(Densidad = xseq), se = TRUE, interval = "confidence")
plot(janka$Densidad, janka$Dureza, xlab= "Densidad", ylab= "Dureza" )
lines(prediccion$fit[,1]~xseq ) # regression line
lines(prediccion$fit[,2]~xseq, lty= 2) # 95% CI limite superior
lines(prediccion$fit[,3]~xseq, lty= 2) # 95% CI limite inferior
Alternativamente podemos usar la función geom_smooth
dentro de ggplot.
ggplot(janka, aes(x = Densidad, y = Dureza))+
geom_point()+
labs(titile = "diagrama de dispersión",
y = "Dureza (lbf)", x = "Densidad (lbs pie-3)")+
geom_smooth(method = "lm", se = TRUE)
## `geom_smooth()` using formula 'y ~ x'
A manera de ejemplo, veamos como se comporta la otra variable que tenemos en nuestra tabla
plot(janka)
Generamos un modelo para la otra variable
janka_lm2 <- lm(Dureza ~ Otra_variable, data = janka)
summary(janka_lm2)
##
## Call:
## lm(formula = Dureza ~ Otra_variable, data = janka)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.2 -496.6 -228.2 416.9 1807.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1856.176 484.790 3.829 0.000527 ***
## Otra_variable -7.270 8.757 -0.830 0.412271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 805.1 on 34 degrees of freedom
## Multiple R-squared: 0.01986, Adjusted R-squared: -0.008963
## F-statistic: 0.6891 on 1 and 34 DF, p-value: 0.4123
ggplot(janka, aes(x = Otra_variable, y = Dureza))+
geom_point()+
labs(titile = "diagrama de dispersión",
y = "Dureza (lbf)", x = "Otra variable")+
geom_smooth(method = "lm", se = TRUE)
## `geom_smooth()` using formula 'y ~ x'
Tal como hemos visto hasta ahora, la función lm()
realiza un análisis de minimos cuadrados el cual asume que los residuales tienen una distribución cercana a la normal. La varianza (o desviación estandar) debe ser aproximadamente constante a lo largo de la variable de respuesta.
Podemos usar la función plot()
en nuestro modelo para evaluar visualmente si los residuales del modelo cumplen con estas condiciones
plot(janka_lm, which = c(1:2))
La gráfica de la izquierda nos muestra dos patrones: la variabilidad no es constante sino que se incrementa de izquierda a derecha conforme los valores de dureza se incrementan y hay una curvatura inexplicada. Por otro lado, la gráfica de quantiles de la derecha nos indica que si bien, hay un buen ajuste a la distribución normal, hay algunos valores que pudieran considerarse extremos.
Para mejorar nuestro modelo podemos hacer una transformación usando la raiz cuadrada.
janka.rc <- lm(sqrt(Dureza) ~ Densidad, data = janka)
summary(janka.rc)
##
## Call:
## lm(formula = sqrt(Dureza) ~ Densidad, data = janka)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5387 -1.0912 -0.2802 1.0259 4.8130
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.25867 1.09397 2.065 0.0466 *
## Densidad 0.75795 0.02296 33.016 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.844 on 34 degrees of freedom
## Multiple R-squared: 0.9698, Adjusted R-squared: 0.9689
## F-statistic: 1090 on 1 and 34 DF, p-value: < 2.2e-16
library(patchwork)
modelo_1 <- ggplot(janka, aes(x = Densidad, y = Dureza))+
geom_point()+
geom_smooth(method="lm", se = TRUE)
modelo_2 <- ggplot(janka, aes(x = Densidad, y = sqrt(Dureza)))+
geom_point()+
geom_smooth(method = "lm", se = TRUE)
modelo_1 + modelo_2
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Evualuación del modelo:
plot(janka.rc, which = c(1:2))
Si bien el modelo logro una mejoria con la transformación, aun hay un patrón con IC angostos con los valores ajustados bajos y anchos a mayor valor.
En este ejercicio, vamos a identificar si existe correlación entre los parametros que se midieron, especificamente para buscar si existe relación entre la longitud del pico y la longitud de la aleta. Para esto, usaremos solamente la especie Gentoo
pin <- read_csv("Data/penguins_size.csv") %>%
filter(species == "Gentoo")
## Parsed with column specification:
## cols(
## species = col_character(),
## island = col_character(),
## culmen_length_mm = col_double(),
## culmen_depth_mm = col_double(),
## flipper_length_mm = col_double(),
## body_mass_g = col_double(),
## sex = col_character()
## )
#quitar casos con NA
pin <- pin[complete.cases(pin),]
Como primer paso, buscaremos si hay correlación entre alguna de las variables
plot(pin)
Ajuste del modelo para predecir la longitud de la aleta en función de la longitud del pico.
mod_pin <- lm(flipper_length_mm ~ culmen_length_mm, data = pin)
coef(mod_pin)
## (Intercept) culmen_length_mm
## 150.798118 1.397386
summary(mod_pin)
##
## Call:
## lm(formula = flipper_length_mm ~ culmen_length_mm, data = pin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.4316 -3.2787 0.0728 3.7300 11.2889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 150.7981 6.9428 21.720 <2e-16 ***
## culmen_length_mm 1.3974 0.1457 9.589 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.937 on 118 degrees of freedom
## Multiple R-squared: 0.438, Adjusted R-squared: 0.4332
## F-statistic: 91.95 on 1 and 118 DF, p-value: < 2.2e-16
Evaluación de los residuales
plot(mod_pin, which = c(1:2))
ggplot(pin, aes(x = culmen_length_mm, y = flipper_length_mm))+
geom_point()+
geom_smooth(method = "lm", se = TRUE)
## `geom_smooth()` using formula 'y ~ x'
Con base en estos resultados ¿crees que un modelo lineal es el que mas se ajusta a estos datos?
El paquete reports
proporciona una referencia sobre como reportar los resultados de un modelo lineal lm()
:
install.packages("devtools")
devtools::install_github("easystats/report", force = TRUE)
library(report)
report(mod_pin)
## We fitted a linear model (estimated using OLS) to predict flipper_length_mm with culmen_length_mm (formula = flipper_length_mm ~ culmen_length_mm). Standardized parameters were obtained by fitting the model on a standardized version of the dataset. Effect sizes were labelled following Funder's (2019) recommendations.
##
## The model explains a significant and substantial proportion of variance (R2 = 0.44, F(1, 118) = 91.95, p < .001, adj. R2 = 0.43). The model's intercept, corresponding to flipper_length_mm = 0 and culmen_length_mm = 0, is at 150.80 (SE = 6.94, 95% CI [137.05, 164.55], p < .001). Within this model:
##
## - The effect of culmen_length_mm is positive and can be considered as large and significant (beta = 1.40, SE = 0.15, std. beta = 0.66, p < .001).
La correlación lineal y la regresión lineal simple son métodos estadísticos que estudian la relación lineal existente entre dos variables.
La correlación cuantifica como estan relacionadas dos variables, meintras que la regresión lieal consiste en generar una ecuación (modelo) que, basándose en la relación existente entre ambas variables, permita predecir el valor de una a partir de la otra.
El cálculo de la correlación entre dos variables es independiente del orden o asignación de cada variable a X e Y, mide únicamente la relación entre ambas sin considerar dependencias. En el caso de la regresión lineal, el modelo varía según qué variable se considere dependiente de la otra (lo cual no implica causa-efecto).
A nivel experimental, la correlación se suele emplear cuando ninguna de las variables se ha controlado, simplemente se han medido ambas y se desea saber si están relacionadas. En el caso de estudios de regresión lineal, es más común que una de las variables se controle (tiempo, concentración de reactivo, temperatura…) y se mida la otra.
Por norma general, los estudios de correlación lineal preceden a la generación de modelos de regresión lineal. Primero se analiza si ambas variables están correlacionadas y, en caso de estarlo, se procede a generar el modelo de regresión.
Para estudiar la relación lineal existente entre dos variables continuas es necesario disponer de parámetros que permitan cuantificar dicha relación. Uno de estos parámetros es la covarianza, que indica el grado de variación conjunta de dos variables aleatorias.
Para poder hacer comparaciones se estandariza la covarianza, generando lo que se conoce como coeficientes de correlación. Existen diferentes tipos, de entre los que destacan el coeficiente de Pearson, Rho de Spearman y Tau de Kendall.
Se emplean como medida de fuerza de asociación (tamaño del efecto): * 0: asociación nula. * 0.1: asociación pequeña. * 0.3: asociación mediana. * 0.5: asociación moderada. * 0.7: asociación alta. * 0.9: asociación muy alta.
Para este ejercicio, usaremos la base de datos de metabolitos en Haliotis fulgens expuesto a un incremento en la temperatura
metabolitos <- read_csv("Data/Hfulgens_metabolitos.csv")
## Parsed with column specification:
## cols(
## Temperature = col_double(),
## Alanine = col_double(),
## Arginine = col_double(),
## Aspartate = col_double(),
## Glutamine = col_double(),
## Homarine = col_double(),
## Isoleucine = col_double(),
## Lactate = col_double(),
## Succinate = col_double(),
## Valine = col_double()
## )
Como una primera exploración de los datos, podemos usar la función plot()
plot(metabolitos)
Correlacion entre la valine y la isoleucina
ggplot(metabolitos, aes(x = Valine, y = Isoleucine))+
geom_point()
El diagrama de dispersión nos indica una posible relación lieal positiva entre ambas variables.
Calculo de la correlación
Para esto, podemos utilizar la funcion cor()
cor(x = metabolitos$Valine, y = metabolitos$Isoleucine, method="pearson")
## [1] 0.8758761
El test muestra una correlación alta (>0.8). Sin embargo para poder considerar que existe realmente una correlación entre las dos variables es necesario calcular su significancia, de lo contrario podría deberse al azar.
Para esto usamos la función cor.test()
cor.test(x = metabolitos$Valine, y = metabolitos$Isoleucine, method="pearson", conf.level = 0.95, alternative = "two.sided")
##
## Pearson's product-moment correlation
##
## data: metabolitos$Valine and metabolitos$Isoleucine
## t = 9.942, df = 30, p-value = 5.238e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7589861 0.9380712
## sample estimates:
## cor
## 0.8758761
El coeficiente de correlación calculado es significativo
**Calculo del coeficiente de determinación \(R^2\)
R2_pearson <- cor(x = metabolitos$Valine, y = metabolitos$Isoleucine, method="pearson")
R2_pearson <- R2_pearson^2
R2_pearson
## [1] 0.7671589
Por lo tanto, podemos concluir que existe una correlación significativa entre la concentración de la Alanina con la Isoleucina (r = 0.87; P < 0.0001), con un tamaño del efecto alto (\(R^2 = 0.76\))
Para esto necesitaremos instalar los siguientes paquetes:
install.packages("apaTables")
install.packages("corrr")
Uso de apaTables
para generar matrices de correlación
library(apaTables)
#no se muestra el resultado
apa.cor.table(metabolitos[,-1], filename = "Tabla_correlacion.doc", show.conf.interval = TRUE)
uso de corrr
para visualizar las correlaciones
library(corrr)
## Warning: package 'corrr' was built under R version 4.0.2
correlate(metabolitos[,-1], method = "pearson", quiet = TRUE) %>%
network_plot(legend = TRUE, min_cor = 0)
A continuación trabajaremos con datos de diversos parámetros en peces de agua dulce. Para este ejercicio usaremos solamente los datos de tilapia.
Queremos generar un modelo que nos permita predecir el peso del pez a partir del resto de los paramtros.
El primer paso para establecer un modelo lineal múltiple es estudiar la relación que existe entre variables. Esta relación en critica a la hora de identificar los mejores predictores para el modelo (p. ej, que variables presentan relación de tipo no lineal) y para identificar colinialidad entre predictores.
Para esto se puede hace representación gfráficas entre las variables o calcular el coeficiente de correlación:
peces <- read_csv("Data/peces_regresion.csv")
## Parsed with column specification:
## cols(
## Especie = col_character(),
## Peso = col_double(),
## Longitud1 = col_double(),
## Longitud2 = col_double(),
## Longitud3 = col_double(),
## Longitud4 = col_double(),
## Alto = col_double(),
## Ancho = col_double()
## )
tilapia <- peces %>%
filter(Especie == "Tilapia")
plot(tilapia)
Podemos usar la funcion pairs.panel()
del paquete psych
para visualziar los diagramas de disperción entre todas las variables asi como su valor de correlación.
install.packages("psych")
library(psych)
## Warning: package 'psych' was built under R version 4.0.2
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
pairs.panels(tilapia[,-1], stars = TRUE)
El paquete Psych tambien contien la función multi.hist()
la cual nos ayudaria a evaluar la distribución de los datos:
multi.hist(x = tilapia[,-1], dcol = c("blue", "red"))
Donde la linea roja corresponde a la distribución normal ajustada y la linea azul a la densidad observada.
Del análisis preliminar se pueden extraer las siguientes conclusiones:
tilapia_lm <- lm(Peso ~ Longitud1 + Longitud2 + Longitud3 + Longitud4 + Ancho + Alto, data = tilapia)
summary(tilapia_lm)
##
## Call:
## lm(formula = Peso ~ Longitud1 + Longitud2 + Longitud3 + Longitud4 +
## Ancho + Alto, data = tilapia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -164.544 -27.502 -4.932 23.936 115.699
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -939.0144 165.4496 -5.676 4.4e-06 ***
## Longitud1 16.4030 50.2746 0.326 0.74665
## Longitud2 18.9733 55.0268 0.345 0.73282
## Longitud3 -21.9432 43.2725 -0.507 0.61606
## Longitud4 0.8675 5.6808 0.153 0.87972
## Ancho 55.4062 44.2011 1.254 0.22039
## Alto 62.9731 19.8432 3.174 0.00364 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 54.89 on 28 degrees of freedom
## Multiple R-squared: 0.9433, Adjusted R-squared: 0.9312
## F-statistic: 77.64 on 6 and 28 DF, p-value: 3.817e-16
El modelo resultante usando todas las variables introducidas como predictores tiene una \(R^2\) alta (0.94), lo que implica que el modelo es capaz de explicar el 94% de la variabilidad observada en los valores de peso.
Es improtante notar que todas las variables menos Alto no son significativos, lo que es un indicativo de que podrían no contribuir al modelo.
Basado en los resultados anteriores, podemos concluir que solamente una variable no da la mejor predicción.
Manualmente, se puede comparar si un modelo mas sencillo es mejor a un modelo mas complejo usando la función anova()
tilapia_lm2 <- lm(Peso ~ Alto, data = tilapia)
summary(tilapia_lm2)
##
## Call:
## lm(formula = Peso ~ Alto, data = tilapia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -147.362 -35.831 -5.472 29.357 115.326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -941.559 74.896 -12.57 3.93e-14 ***
## Alto 102.705 4.893 20.99 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.06 on 33 degrees of freedom
## Multiple R-squared: 0.9303, Adjusted R-squared: 0.9282
## F-statistic: 440.5 on 1 and 33 DF, p-value: < 2.2e-16
anova(tilapia_lm, tilapia_lm2)
## Analysis of Variance Table
##
## Model 1: Peso ~ Longitud1 + Longitud2 + Longitud3 + Longitud4 + Ancho +
## Alto
## Model 2: Peso ~ Alto
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 84368
## 2 33 103699 -5 -19332 1.2831 0.299
Para encontrar el mejor modelo, es posible utilizar la función step()
el cual utiliza la strategia stepwise mixto y nos indica el mejor modelo basado en el AIC.
step(tilapia_lm, direction = "both", trace = 1)
## Start: AIC=286.57
## Peso ~ Longitud1 + Longitud2 + Longitud3 + Longitud4 + Ancho +
## Alto
##
## Df Sum of Sq RSS AIC
## - Longitud4 1 70.3 84438 284.60
## - Longitud1 1 320.8 84688 284.70
## - Longitud2 1 358.2 84726 284.71
## - Longitud3 1 774.8 85143 284.89
## - Ancho 1 4734.4 89102 286.48
## <none> 84368 286.57
## - Alto 1 30346.3 114714 295.32
##
## Step: AIC=284.59
## Peso ~ Longitud1 + Longitud2 + Longitud3 + Ancho + Alto
##
## Df Sum of Sq RSS AIC
## - Longitud1 1 313 84751 282.72
## - Longitud2 1 326 84764 282.73
## - Longitud3 1 732 85170 282.90
## <none> 84438 284.60
## - Ancho 1 5337 89775 284.74
## + Longitud4 1 70 84368 286.57
## - Alto 1 37700 122138 295.51
##
## Step: AIC=282.72
## Peso ~ Longitud2 + Longitud3 + Ancho + Alto
##
## Df Sum of Sq RSS AIC
## - Longitud3 1 494 85245 280.93
## - Longitud2 1 1555 86306 281.36
## <none> 84751 282.72
## - Ancho 1 5461 90211 282.91
## + Longitud1 1 313 84438 284.60
## + Longitud4 1 62 84688 284.70
## - Alto 1 39035 123786 293.98
##
## Step: AIC=280.93
## Peso ~ Longitud2 + Ancho + Alto
##
## Df Sum of Sq RSS AIC
## - Ancho 1 5000 90244 280.92
## <none> 85245 280.93
## - Longitud2 1 8012 93256 282.07
## + Longitud3 1 494 84751 282.72
## + Longitud1 1 75 85170 282.90
## + Longitud4 1 29 85216 282.92
## - Alto 1 38542 123787 291.98
##
## Step: AIC=280.92
## Peso ~ Longitud2 + Alto
##
## Df Sum of Sq RSS AIC
## <none> 90244 280.92
## + Ancho 1 5000 85245 280.93
## + Longitud4 1 762 89482 282.63
## + Longitud1 1 468 89776 282.74
## + Longitud3 1 33 90211 282.91
## - Longitud2 1 13455 103699 283.79
## - Alto 1 65408 155652 298.00
##
## Call:
## lm(formula = Peso ~ Longitud2 + Alto, data = tilapia)
##
## Coefficients:
## (Intercept) Longitud2 Alto
## -1013.12 16.35 71.77
Finalmente, generamos un modelo utilizando las variables Alto y Longitud2 y lo comparamos con el modelo anterior
tilapia_lm3 <- lm(Peso ~ Alto + Longitud2, data = tilapia)
summary(tilapia_lm3)
##
## Call:
## lm(formula = Peso ~ Alto + Longitud2, data = tilapia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -168.537 -24.862 -5.562 22.370 107.114
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1013.115 78.150 -12.964 2.79e-14 ***
## Alto 71.769 14.902 4.816 3.39e-05 ***
## Longitud2 16.348 7.485 2.184 0.0364 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.1 on 32 degrees of freedom
## Multiple R-squared: 0.9394, Adjusted R-squared: 0.9356
## F-statistic: 247.8 on 2 and 32 DF, p-value: < 2.2e-16
anova(tilapia_lm2, tilapia_lm3)
## Analysis of Variance Table
##
## Model 1: Peso ~ Alto
## Model 2: Peso ~ Alto + Longitud2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 33 103699
## 2 32 90244 1 13455 4.771 0.03638 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
To compare the fits of two models, you can use the anova() function with the regression objects as two separate arguments. The anova() function will take the model objects as arguments, and return an ANOVA testing whether the more complex model is significantly better at capturing the data than the simpler model. If the resulting p-value is sufficiently low (usually less than 0.05), we conclude that the more complex model is significantly better than the simpler model, and thus favor the more complex model. If the p-value is not sufficiently low (usually greater than 0.05), we should favor the simpler model.
Cada una de las pendientes de un modelo de regresión lineal múltiple (coeficientes parciales de regresión de los predictores) se define del siguiente modo: Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable (Y) varía en promedio tantas unidades como indica la pendiente. Para este ejemplo, por cada unidad que aumenta el predictor Alto, el peso aumenta en promedio 71.76 unidades, manteniéndose constantes el resto de predictores.
La relación lieal entre los predictores numéricos y la variable de respuesta se puede evaluar con diagramas de dispersión entre la variable dependiente y cada uno de los predictores o con diagramas de dispersión entre cada uno de los predictores y los residuos del modelo. Si la relación es lineal, los residuos deben de distribuirse aleatoriamente en torno a 0 con una variabilidad constante a lo largo del eje X.
library(patchwork)
Alto_res <- ggplot(tilapia, aes(x = Alto, y = tilapia_lm3$residuals))+
geom_point()+
geom_smooth(color = "firebrick")+
geom_hline(yintercept = 0)+
theme_bw()
Longitud2_res <- ggplot(tilapia, aes(x = Longitud2, y = tilapia_lm3$residuals))+
geom_point()+
geom_smooth(color = "firebrick")+
geom_hline(yintercept = 0)+
theme_bw()
Alto_res + Longitud2_res
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
qqnorm(tilapia_lm3$residuals)
qqline(tilapia_lm3$residuals)
shapiro.test(tilapia_lm3$residuals)
##
## Shapiro-Wilk normality test
##
## data: tilapia_lm3$residuals
## W = 0.93388, p-value = 0.03652
El modelo lineal múltiple:
\(Peso = -1013.115 + 71.76 Alto + 16.34 Longitud2\)
es capaz de explicar el 94% de la variabilidad observada en el peso. El test F muestra que es significatvio (\(P < 2.22^{-16}\)).
o
library(report)
report(tilapia_lm3)
## We fitted a linear model (estimated using OLS) to predict Peso with Alto and Longitud2 (formula = Peso ~ Alto + Longitud2). Standardized parameters were obtained by fitting the model on a standardized version of the dataset. Effect sizes were labelled following Funder's (2019) recommendations.
##
## The model explains a significant and substantial proportion of variance (R2 = 0.94, F(2, 32) = 247.83, p < .001, adj. R2 = 0.94). The model's intercept, corresponding to Peso = 0, Alto = 0 and Longitud2 = 0, is at -1013.12 (SE = 78.15, 95% CI [-1172.30, -853.93], p < .001). Within this model:
##
## - The effect of Alto is positive and can be considered as large and significant (beta = 71.77, SE = 14.90, std. beta = 0.67, p < .001).
## - The effect of Longitud2 is positive and can be considered as small and significant (beta = 16.35, SE = 7.48, std. beta = 0.31, p < .05).
NOTA: Puedes encontrar una explicación mas completa sobre el procedimeinto aqui
Los modelos lineales tienen la ventaja de ser fácilmente interpretables, sin embargo, pueden tener limitaciones importantes en capacidad predictiva. Esto se debe a que, la asunción de linealidad, es con frecuencia una aproximación demasiado simple para describir las relaciones reales entre variables. A continuación, se describen métodos que permiten relajar la condición de linealidad intentando mantener al mismo tiempo una interpretabilidad alta.
Los modelos lineales son relativamente simples de describir e implementar, y poseen la ventaja frente a otros modelos de ser fáciles de interpretar. Sim embargo, pueden estar limitados en cuanto a capacidad predictiva si la relación entre variables es más compleja.
Se reemplaza el modelo lineal:
\(y = \beta_0 + \beta_1x_i + \epsilon_i\)
con una función polinomial:
\(y = \beta_0 + \beta_ix_i + \beta_2x^{2}_i + \beta_3x^{3}_i\)
Para un polinomio de grado d lo suficientemente elevado, la regresión polinómica permite ajustar una curva extremadamente no lineal, aunque no es habitual utilizar un polinomio mayor de grado 3 o 4 para evitar una flexibilidad excesiva u oscilaciones indeseables, sobre todo en los extremos del predictor X.
Para estimar una ecuación de regresión polinómica:
visualización de los datos
Para este ejercicio, vamos a usar el resto de los peces en el archivo peces_regresion.csv
de manera que ahora tenemos que remover las tilapias
peces <- peces %>%
filter(Especie != "Tilapia")
plot(peces[,-1])
Se observa que la relación entre el peso y la Logitud1 no tiene una correlación lineal.
A. modelo lineal simple:
pez_lm1 <- lm(Peso ~ Longitud1, data = peces)
summary(pez_lm1)
##
## Call:
## lm(formula = Peso ~ Longitud1, data = peces)
##
## Residuals:
## Min 1Q Median 3Q Max
## -146.08 -51.14 -14.87 47.34 325.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -608.488 31.138 -19.54 <2e-16 ***
## Longitud1 38.481 1.237 31.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.53 on 91 degrees of freedom
## Multiple R-squared: 0.9141, Adjusted R-squared: 0.9132
## F-statistic: 968.5 on 1 and 91 DF, p-value: < 2.2e-16
ggplot(peces, aes(x = Longitud1, y = Peso))+
geom_point()+
geom_smooth(method = "lm", se = TRUE)+
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
B. Modelo polinomial, forma equivocada de escribir el modelo:
pez_lm2 <- lm(Peso ~ Longitud1 + Longitud1^2, data = peces)
summary(pez_lm2)
##
## Call:
## lm(formula = Peso ~ Longitud1 + Longitud1^2, data = peces)
##
## Residuals:
## Min 1Q Median 3Q Max
## -146.08 -51.14 -14.87 47.34 325.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -608.488 31.138 -19.54 <2e-16 ***
## Longitud1 38.481 1.237 31.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.53 on 91 degrees of freedom
## Multiple R-squared: 0.9141, Adjusted R-squared: 0.9132
## F-statistic: 968.5 on 1 and 91 DF, p-value: < 2.2e-16
C. Forma correcta: opción 1
pez_lm3 <- lm(Peso ~ Longitud1 + I(Longitud1^2), data = peces)
summary(pez_lm3)
##
## Call:
## lm(formula = Peso ~ Longitud1 + I(Longitud1^2), data = peces)
##
## Residuals:
## Min 1Q Median 3Q Max
## -185.752 -21.435 -3.306 17.927 231.189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 129.71739 63.14239 2.054 0.0428 *
## Longitud1 -22.65056 5.04028 -4.494 2.07e-05 ***
## I(Longitud1^2) 1.15052 0.09377 12.269 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 55.69 on 90 degrees of freedom
## Multiple R-squared: 0.9679, Adjusted R-squared: 0.9671
## F-statistic: 1355 on 2 and 90 DF, p-value: < 2.2e-16
D.Forma correcta: opción 2:
pez_lm4 <- lm(Peso ~ poly(Longitud1, degree = 2, raw = TRUE), data = peces)
summary(pez_lm4)
##
## Call:
## lm(formula = Peso ~ poly(Longitud1, degree = 2, raw = TRUE),
## data = peces)
##
## Residuals:
## Min 1Q Median 3Q Max
## -185.752 -21.435 -3.306 17.927 231.189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 129.71739 63.14239 2.054 0.0428
## poly(Longitud1, degree = 2, raw = TRUE)1 -22.65056 5.04028 -4.494 2.07e-05
## poly(Longitud1, degree = 2, raw = TRUE)2 1.15052 0.09377 12.269 < 2e-16
##
## (Intercept) *
## poly(Longitud1, degree = 2, raw = TRUE)1 ***
## poly(Longitud1, degree = 2, raw = TRUE)2 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 55.69 on 90 degrees of freedom
## Multiple R-squared: 0.9679, Adjusted R-squared: 0.9671
## F-statistic: 1355 on 2 and 90 DF, p-value: < 2.2e-16
Ahora podemos comparar el modelo lineal con el modelo que generamos con un polinomia de grado 2
anova(pez_lm1, pez_lm3)
## Analysis of Variance Table
##
## Model 1: Peso ~ Longitud1
## Model 2: Peso ~ Longitud1 + I(Longitud1^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 91 745873
## 2 90 279084 1 466789 150.53 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Recordemos que la hipotesis nula es que el modelo mas complejo no aporta mas información que el modelo mas sencillo, por lo tanto, los resultados nos indica una mejora significativa del ajuste polinomial
plot(pez_lm3, which = c(1:2))
Ahora probemos efecto de agregar un mayor grado a la ecuación:
pez_lm5 <- lm(Peso ~ poly(Longitud1, degree = 3, raw = TRUE), data = peces)
summary(pez_lm5)
##
## Call:
## lm(formula = Peso ~ poly(Longitud1, degree = 3, raw = TRUE),
## data = peces)
##
## Residuals:
## Min 1Q Median 3Q Max
## -189.655 -25.457 -0.578 18.112 220.024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 285.83743 140.61203 2.033 0.0451
## poly(Longitud1, degree = 3, raw = TRUE)1 -44.58656 18.36712 -2.428 0.0172
## poly(Longitud1, degree = 3, raw = TRUE)2 2.10251 0.77237 2.722 0.0078
## poly(Longitud1, degree = 3, raw = TRUE)3 -0.01275 0.01027 -1.242 0.2176
##
## (Intercept) *
## poly(Longitud1, degree = 3, raw = TRUE)1 *
## poly(Longitud1, degree = 3, raw = TRUE)2 **
## poly(Longitud1, degree = 3, raw = TRUE)3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 55.52 on 89 degrees of freedom
## Multiple R-squared: 0.9684, Adjusted R-squared: 0.9673
## F-statistic: 909.4 on 3 and 89 DF, p-value: < 2.2e-16
anova(pez_lm4, pez_lm5)
## Analysis of Variance Table
##
## Model 1: Peso ~ poly(Longitud1, degree = 2, raw = TRUE)
## Model 2: Peso ~ poly(Longitud1, degree = 3, raw = TRUE)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 90 279084
## 2 89 274332 1 4752.4 1.5418 0.2176
Grafica del modelo usando ggplot:
ggplot(peces, aes(x = Longitud1, y = Peso))+
geom_point()+
geom_smooth(method = "lm", formula = y~poly(x,2))+
geom_smooth(method = "lm", formula = y~x, col = "firebrick")+
theme_bw()
Agregando el modelo con 3 grado
ggplot(peces, aes(x = Longitud1, y = Peso))+
geom_point()+
geom_smooth(method = "lm",se = FALSE, formula = y~poly(x,2))+
geom_smooth(method = "lm", se = FALSE, formula = y~poly(x,3), col = "firebrick")+
theme_bw()
Puedes consultar una guia sobre modelos no lineales aqui.