1: Generalidades

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.

1.1: Modelos lineales generales

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.

1.2: Regresión lineal simple

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:

  • \(H_0\): No hay relación lieal entre ambas variables por lo que la pendiente del modelo lineal es cero. \(b = 0\)
  • \(H_a\): Sí hay relación lieneal entre ambas variables por lo que la pendiente del modelo lieal es distinta a cero \(b \neq 0\)

De esta misma forma tambien se aplica a \(a\).

2: Densidad y dureza de madera

2.1: El modelo

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.

2.2: Representación gráfica

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'

2.3: Otra variable

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'

2.4: Validez del modelo

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.

Ejercicio: Pinguinos

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).

3: Correlacion

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.

3.1: Metabolitos

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\))

3.2: Visualización de correlaciones entre multiples variables

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)

4: Correlación lineal multiple

4.1: Analizar la relación entre variables

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:

  • Las variables que tienen mayor relación lineal con el peso son: ancho (r = 0.96), longitud3 y longitud2 (r = 0.95).
  • Las longitudes 1,2 y 3 estan altamente correlacionadas entre ellas, por lo que posiblemente no sea útil introducir ambos predictores en el modelo.

4.2: Generar el modelo

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.

4.3: Selección de los mejores predictores

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.

4.4: Validación de condiciones para la regresión lineal multiple

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

4.5: Conlcusión

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

5: Correlación no lineal

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.

5.1: Regresión polinómica

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.