1 Introducción

En esta cuarta clase, profundizaremos en la aplicación de análisis estadísticos básicos utilizando R. Hasta ahora, hemos aprendido a visualizar y manipular datos usando ggplot2 y dplyr; hoy, daremos un paso más allá al integrar estas habilidades en un contexto de análisis estadístico formal.

La estadística es fundamental para comprender y extraer conclusiones de los datos, y en esta sesión, nos enfocaremos en cómo comparar variables, identificar diferencias significativas entre grupos y evaluar relaciones entre variables cuantitativas. Utilizaremos el dataset “Palmer Penguins,” que contiene información sobre tres especies de pingüinos en la Antártida, como nuestro caso de estudio.

La clase se dividirá en tres secciones principales:

Análisis Exploratorio de Datos (EDA)
Comparación de Medias entre Dos Grupos (t-test)
Evaluación de Relaciones entre Variables Cuantitativas (Regresión Lineal)

A lo largo de cada sección, aplicaremos análisis de supuestos para asegurarnos de que los métodos estadísticos utilizados sean apropiados para nuestros datos. Al final de la clase, serás capaz de aplicar estos conceptos a tus propios conjuntos de datos y llevar a cabo análisis estadísticos de manera efectiva en R.

2 Librerias y Datos

Para esta clase, utilizaremos los paquetes palmerpenguins, tidyverse, corrgram, ggfortify y GGally. Asegúrate de que estén instalados y luego cárgalos en tu entorno de trabajo.

library(tidyverse)
library(palmerpenguins)
library(GGally) #extensión de ggplot
library(corrgram) # para graficar correlaciones
library(ggfortify) # para graficos de supuestos modelos

A continuación, cargaremos el dataset “Palmer Penguins” y le daremos el nombre peng. También haremos una revisión rápida para asegurarnos de que los datos estén listos para el análisis.

# Cargar el dataset y asignarlo a un nuevo objeto llamado 'peng'
peng <- penguins

# Un vistazo rápido a los datos
glimpse(peng)

## Rows: 344
## Columns: 8
## $ species           <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island            <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm     <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g       <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex               <fct> male, female, female, NA, female, male, female, male…
## $ year              <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…

Este dataset contiene información sobre tres especies de pingüinos (Adelie, Chinstrap y Gentoo), así como medidas relacionadas con la longitud y profundidad del pico, la longitud de la aleta, la masa corporal y el sexo de los pingüinos. Ahora que tenemos los datos y librerías cargados, estamos listos para comenzar con el Análisis Exploratorio de Datos (EDA).

3 Análisis Exploratorio de Datos (EDA)

El Análisis Exploratorio de Datos (EDA) es un paso crucial antes de aplicar cualquier modelo estadístico o de realizar pruebas de hipótesis. Nos permite comprender la estructura y las características fundamentales del conjunto de datos, identificar patrones, detectar anomalías, y formular hipótesis iniciales. A través del EDA, podemos descubrir relaciones interesantes entre las variables y asegurarnos de que nuestros datos estén listos para el análisis formal.

3.1 Análisis Descriptivo de las Variables

Comencemos revisando las estadísticas descriptivas de las variables numéricas para conocer sus medias, medianas, mínimos, máximos, y rangos intercuartílicos.

# Resumen estadístico de las variables numéricas
summary(peng)

##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
##                                  Mean   :43.92   Mean   :17.15  
##                                  3rd Qu.:48.50   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##                                  NA's   :2       NA's   :2      
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197.0     Median :4050   NA's  : 11   Median :2008  
##  Mean   :200.9     Mean   :4202                Mean   :2008  
##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231.0     Max.   :6300                Max.   :2009  
##  NA's   :2         NA's   :2

Tambien podemos crear nuestra propia tabla resumen:

peng %>% 
  group_by(species) %>% 
  summarise(across(c(body_mass_g, ends_with("_mm")), 
                   # Aquí se define una lista de funciones que se aplicarán a cada columna seleccionada por across()
                   list(mean = ~mean(.x, na.rm = TRUE), 
                        sd = ~sd(.x, na.rm = TRUE))))

3.2 Limpieza de datos faltantes

Antes de realizar cualquier análisis, es importante asegurarnos de que no haya valores faltantes en nuestro conjunto de datos, ya que estos podrían afectar los resultados de nuestro análisis estadístico.

# Eliminar filas con datos faltantes
peng <- drop_na(peng)

# Verificar que no queden datos faltantes
sum(is.na(peng))  # Esto debería retornar 0 si todos los valores faltantes fueron eliminados

## [1] 0

3.3 Exploración visual de las variables

El paquete GGally es una extensión de ggplot2 que facilita la creación de visualizaciones más complejas y combinadas. Uno de los principales atractivos de GGally es la función ggpairs(), que nos permite generar de manera sencilla una matriz de gráficos de pares para visualizar relaciones entre múltiples variables.

Utilizaremos ggpairs() para explorar las relaciones entre las variables continuas (body_mass_g, bill_length_mm, bill_depth_mm, etc.) de nuestro conjunto de datos peng, distinguiendo a los pingüinos según su especie (species). Esto nos brindará una visión general de cómo se relacionan las diferentes mediciones y si existen diferencias notables entre las especies.

peng %>%
  select(species, body_mass_g, ends_with("_mm")) %>% 
  GGally::ggpairs(aes(color = species)) +
  scale_colour_manual(values = c("darkorange","purple","cyan4")) +
  scale_fill_manual(values = c("darkorange","purple","cyan4"))+
  theme_bw()

Para profundizar en el análisis, exploraremos cómo se relacionan las variables dentro de la especie “Adelie” en función del sexo (male y female). Esta exploración nos permitirá ver si existen diferencias en las características morfológicas entre machos y hembras de la misma especie, lo que puede ser relevante para entender posibles dimorfismos sexuales.

peng %>%
  filter(species == "Adelie") %>% 
  select(sex, body_mass_g, ends_with("_mm")) %>% 
  GGally::ggpairs(aes(color = sex)) +
  scale_colour_manual(values = c("skyblue2","pink2")) +
  scale_fill_manual(values = c("male"= "skyblue2",
                               "female" = "pink2"))+
  theme_bw()

3.4 Distribuciones Univariadas de las Variables Continuas

las distribuciones univariadas de nuestras variables continuas (body_mass_g, bill_length_mm, bill_depth_mm, etc.) para cada especie de pingüino. Esto nos ayudará a identificar la variabilidad, la presencia de outliers, y las diferencias entre las especies.

3.4.1 Transformación del Dataset a Formato Largo

La función pivot_longer() del paquete tidyverse se utiliza para transformar un dataset de formato “ancho” (wide) a formato “largo” (long). Esta transformación es especialmente útil cuando queremos visualizar o analizar múltiples variables de manera simultánea.

# Transformar el dataset a formato largo
peng_long <- peng %>% 
  select(species, body_mass_g, ends_with("_mm")) %>% 
  pivot_longer(cols = -species, # Columnas a convertir
               names_to = "variable",  # Nombre de la nueva columna que almacenará los nombres de las variables originales
               values_to = "valor")   # Nombre de la nueva columna que almacenará los valores

peng_long

cols = -species: Indica que queremos transformar todas las columnas excepto species a formato largo.
names_to = “variable”: Crea una nueva columna llamada variable que contiene los nombres originales de las variables (body_mass_g, bill_length_mm, etc.).
values_to = “valor”: Crea una nueva columna llamada valor que contiene los valores correspondientes a cada observación.

3.4.2 Visualización del Data Set Transformado

Una vez que nuestro dataset ha sido transformado a formato largo, podemos crear visualizaciones que nos permitan comparar y analizar las variables continuas según la especie de los pingüinos. Esto nos ayudará a identificar diferencias y patrones importantes en las mediciones.

Boxplot de las Variables Continuas por Especie

El boxplot es una herramienta eficaz para visualizar la distribución, la mediana, y la presencia de valores atípicos (outliers) en cada variable continua según la especie de pingüino.

# Boxplot de las variables continuas por especie
peng_long %>% 
  ggplot(aes(y = valor, x = species, color = species)) +
  geom_boxplot(outlier.colour = "red") +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Distribución de Variables Continuas por Especie",
       x = "Especie", y = "Valor de la Variable") +
  theme_minimal() +
  theme(legend.position = "none") # Eliminamos la leyenda redundante

facet_wrap(~variable, scales = “free”): Permite que cada variable tenga su propia escala, haciendo que la visualización sea más clara y comprensible.
outlier.colour = “red”: Resalta los valores atípicos en rojo, facilitando su identificación.

Análisis de Outliers
Los outliers, o valores atípicos, son observaciones que se alejan significativamente del resto de los datos y suelen destacarse en los boxplots. Es importante revisar estos valores porque pueden influir de manera considerable en los resultados de los análisis estadísticos y la interpretación de los datos.

Los outliers pueden surgir por errores de medición, variabilidad natural, o representar casos especiales que merecen atención. Antes de eliminarlos, es fundamental analizar su origen y considerar si son parte inherente del fenómeno que estamos estudiando. En algunos casos, puede ser más apropiado tratarlos por separado o aplicar transformaciones que reduzcan su impacto en los análisis posteriores.

La decisión de cómo manejar los outliers debe basarse en el contexto y los objetivos del estudio.

Histogramas de las Variables Continuas por Especie

Los histogramas nos permiten observar la distribución de cada variable continua para cada especie y evaluar la forma, la dispersión y la concentración de los datos.

# Histogramas de las variables continuas por especie
peng_long %>% 
  ggplot(aes(x = valor, fill = species)) +
  geom_histogram(color = "gray30", bins = 30, alpha = 0.7) +
  facet_grid(species ~ variable, scales = "free") +
  labs(title = "Distribuciones de las Variables Continuas por Especie",
       x = "Valor de la Variable", y = "Frecuencia") +
  theme_minimal() +
  scale_fill_manual(values = c("darkorange", "purple", "cyan4"))

facet_grid(species ~ variable, scales = “free”): Crea una cuadrícula de gráficos con las especies en las filas y las variables en las columnas, permitiendo ver la distribución de cada variable dentro de cada especie.
bins = 30: Ajusta el número de barras en los histogramas para que la distribución sea más detallada.

Gráfico de Densidad de Probabilidad para las Variables Continuas por Especie

El gráfico de densidad de probabilidad es otra excelente herramienta para visualizar la distribución de las variables continuas. A diferencia de los histogramas, los gráficos de densidad muestran una estimación suave de la distribución, lo que facilita la identificación de patrones, asimetrías y diferencias entre los grupos.

# Gráficos de densidad de probabilidad por especie
peng_long %>% 
  ggplot(aes(x = valor, fill = species, color = species)) +
  geom_density(alpha = 0.4) +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Gráfico de Densidad de Probabilidad por Especie",
       x = "Valor de la Variable", y = "Densidad") +
  theme_minimal() +
  scale_fill_manual(values = c("darkorange", "purple", "cyan4")) +
  scale_color_manual(values = c("darkorange", "purple", "cyan4"))

geom_density(alpha = 0.4): Crea el gráfico de densidad y establece la transparencia con alpha = 0.4 para que las áreas superpuestas entre especies sean visibles.
facet_wrap(~variable, scales = “free”): Permite ver la densidad de cada variable por separado y con su propia escala.

3.5 Correlación entre variables cuantitativas

La correlación es una medida que nos indica la fuerza y la dirección de la relación entre dos variables cuantitativas. Calcular y visualizar la matriz de correlación nos permite identificar qué variables están asociadas entre sí, lo que es fundamental para comprender cómo las diferentes características de los pingüinos se relacionan entre sí.

3.5.1 Calculo de la matriz de correlacion

Primero, vamos a calcular la matriz de correlación utilizando las variables cuantitativas del dataset peng.

# Seleccionar las variables cuantitativas y calcular la matriz de correlación
cor_matrix <- peng %>% 
  select(body_mass_g, ends_with("_mm")) %>% 
  cor(use = "complete.obs", # Calcular la correlación ignorando los valores NA
      method = "pearson")   # Método a utilizar

# Visualizar la matriz de correlación
cor_matrix

##                   body_mass_g bill_length_mm bill_depth_mm flipper_length_mm
## body_mass_g         1.0000000      0.5894511    -0.4720157         0.8729789
## bill_length_mm      0.5894511      1.0000000    -0.2286256         0.6530956
## bill_depth_mm      -0.4720157     -0.2286256     1.0000000        -0.5777917
## flipper_length_mm   0.8729789      0.6530956    -0.5777917         1.0000000

3.5.2 Visualización de la Matriz de Correlación

Para facilitar la interpretación, podemos utilizar el paquete corrplot para crear un gráfico visual de la matriz de correlación:

# Crear el gráfico de la matriz de correlación
corrgram(cor_matrix,
         lower.panel = panel.shade,
         upper.panel = panel.cor)

Argumentos corrgram()

panel: Especifica la función de panel que se usa para mostrar los coeficientes de correlación. Algunos ejemplos comunes son:

panel.shade: Sombreado de las correlaciones.
panel.pie: Representa las correlaciones como gráficos de sectores.
panel.ellipse: Usa elipses para mostrar la relación entre variables.
panel.minmax: Muestra minimos y maximos.

lower.panel, upper.panel y diag.panel: Estos argumentos permiten especificar diferentes funciones de panel para la parte inferior y superior del correlograma, respectivamente. Esto te permite personalizar el tipo de visualización que deseas en cada mitad.

text.panel: Define la función para el panel de texto, que usualmente muestra los nombres de las variables. El valor por defecto es panel.txt.

4 Comparación de Medias entre Dos Grupos: Prueba t de Student

En esta sección, aprenderemos a comparar una variable cuantitativa entre dos grupos definidos por una variable cualitativa utilizando la prueba t de Student. Esta prueba nos permite determinar si existe una diferencia significativa en las medias de la variable cuantitativa entre los dos grupos.

¿Qué es una prueba de hipótesis estadística?

Una prueba de hipótesis estadística es un procedimiento utilizado para tomar decisiones o realizar inferencias sobre una población basándose en datos de una muestra. El proceso implica plantear dos hipótesis:

Hipótesis nula (H₀): Representa la afirmación de “no efecto” o “sin diferencia”. Es el supuesto que se considera verdadero hasta que se demuestre lo contrario.
Hipótesis alternativa (H₁): Representa la afirmación de “efecto” o “diferencia” que queremos probar.

La prueba de hipótesis evalúa si hay suficiente evidencia en los datos de la muestra para rechazar la hipótesis nula en favor de la hipótesis alternativa.

Nivel de significancia (α)

El nivel de significancia (α) es el umbral que define el criterio para rechazar la hipótesis nula. Es la probabilidad máxima de cometer un error tipo I, es decir, rechazar H₀ cuando en realidad es verdadera. Los niveles de significancia más comunes son:

α = 0.05 (5%): Existe un 5% de probabilidad de rechazar la hipótesis nula cuando es verdadera. Es el nivel más utilizado en muchas disciplinas.
α = 0.01 (1%): Utilizado cuando se requiere mayor rigurosidad y se tolera un error menor. Por ejemplo, en investigaciones médicas.
α = 0.10 (10%): Menos estricto y más flexible. A veces se usa en investigaciones exploratorias.

Si el valor p (p-value) resultante de la prueba es menor que el nivel de significancia α, se rechaza la hipótesis nula. Por el contrario, si el valor p es mayor que α, no hay suficiente evidencia para rechazar H₀.

Conclusiones de la Prueba de Hipótesis

Rechazar H₀: Hay suficiente evidencia para concluir que existe un efecto o diferencia significativa.
No rechazar H₀: No hay suficiente evidencia para concluir que existe un efecto o diferencia significativa.

4.1 ¿Qué es la Prueba t de Student?

La prueba t de Student es un método estadístico que se utiliza para comparar las medias de dos grupos y determinar si son significativamente diferentes entre sí. Hay dos tipos de pruebas t:

Prueba t para muestras independientes: Se utiliza cuando los dos grupos son independientes.
Prueba t para muestras pareadas: Se utiliza cuando las observaciones están relacionadas (por ejemplo, antes y después de un tratamiento en el mismo grupo).

En este caso, utilizaremos la prueba t para muestras independientes.

4.2 Aplicación de la Prueba t en el Dataset “Palmer Penguins”

Vamos a comparar la masa corporal (body_mass_g) entre los sexos (male y female) de la especie “Adelie” para ver si existe una diferencia significativa.

Primero, filtraremos los datos para quedarnos solo con la especie “Adelie” en un nuevo set de datos:

# Filtrar datos para la especie "Adelie"
peng_adelie <- peng %>% 
  filter(species == "Adelie")

4.2.1 Verificar Supuestos de la Prueba t

Antes de realizar la prueba, es importante verificar que se cumplen los supuestos:

Independencia de las Observaciones: Las observaciones deben ser independientes entre sí. Esto significa que el valor de una observación no debe influir o estar relacionado con el valor de otra. En el caso de la prueba t para muestras independientes, esto implica que los dos grupos deben ser completamente independientes. En este caso asumiremos que el supuesto se cumple.
Normalidad: La variable dependiente debe estar distribuida normalmente dentro de cada grupo que se está comparando. Comprobamos si body_mass_g sigue una distribución normal dentro de cada grupo (male y female). Para comprobar utilizamos la prueba de normalidad de Shapiro-Wilk o un Gráfico Q-Q (quantile-quantile plot):

# Creamos dos set, uno para cada prueba

# Machos
peng_adelie_m <- peng_adelie %>% 
  filter(sex == "male") %>% select(body_mass_g) %>% 
  as_vector()
  
# Hembras
peng_adelie_f <- peng_adelie %>% 
  filter(sex == "female") %>% select(body_mass_g) %>% 
  as_vector()

# Prueba machos
shapiro.test(peng_adelie_m)

## 
##  Shapiro-Wilk normality test
## 
## data:  peng_adelie_m
## W = 0.98269, p-value = 0.416

# Prueba hembras
shapiro.test(peng_adelie_f)

## 
##  Shapiro-Wilk normality test
## 
## data:  peng_adelie_f
## W = 0.97684, p-value = 0.1985

Interpretación de la Prueba de Shapiro-Wilk

La prueba de Shapiro-Wilk es un test estadístico que se utiliza para evaluar si una muestra sigue una distribución normal. Es uno de los métodos más comunes para verificar la normalidad de los datos, especialmente cuando trabajamos con análisis estadísticos que asumen que las variables siguen una distribución normal (como la prueba t de Student o ANOVA).

Cómo se realiza la prueba Cuando aplicas la función shapiro.test() en R, obtienes dos resultados principales:

W: El estadístico de la prueba. Indica la proximidad de la distribución de los datos a una distribución normal.
p-value: El valor p indica la significancia estadística del test.

Cómo interpretar los resultados

Hipótesis nula (H0): Los datos siguen una distribución normal.

Hipótesis alternativa (H1): Los datos no siguen una distribución normal.

Interpretación del valor p

Si el valor p (p-value) es mayor que 0.05: No hay suficiente evidencia para rechazar la hipótesis nula, por lo que podemos asumir que los datos siguen una distribución normal.

Si el valor p (p-value) es menor o igual a 0.05: Rechazamos la hipótesis nula, lo que indica que los datos no siguen una distribución normal.

Homocedasticidad: Las varianzas de la variable dependiente deben ser aproximadamente iguales en los dos grupos que se están comparando. La prueba de Levene es una prueba estadística que se utiliza para verificar el supuesto de homogeneidad de varianzas entre dos o más grupos. Verificamos si las varianzas de body_mass_g son similares entre los grupos_ con la prueba de levene:

# Cargar el paquete 'car' si no está instalado
library(car)

# Prueba de Levene para verificar la homocedasticidad
leveneTest(body_mass_g ~ sex, data = peng_adelie)

Interpretación de la Prueba de Levene

La prueba de Levene es un test estadístico que se utiliza para evaluar si las varianzas de dos o más grupos son iguales. Es especialmente útil para verificar el supuesto de homocedasticidad antes de realizar pruebas que asumen igualdad de varianzas, como la prueba t de Student o ANOVA.

Cómo se realiza la prueba
Cuando aplicas la función leveneTest() en R, obtienes tres resultados principales:

F-value: El estadístico F de la prueba, que mide la relación entre la variabilidad entre grupos y dentro de los grupos.
Df: Los grados de libertad asociados con la prueba.
p-value: El valor p indica la significancia estadística del test.

Cómo interpretar los resultados

Hipótesis nula (H0): Las varianzas de los grupos son iguales (homocedasticidad).

Hipótesis alternativa (H1): Las varianzas de los grupos son diferentes (heterocedasticidad).

Interpretación del valor p

Si el valor p (p-value) es mayor que 0.05: No hay suficiente evidencia para rechazar la hipótesis nula, por lo que podemos asumir que las varianzas de los grupos son iguales.

Si el valor p (p-value) es menor o igual a 0.05: Rechazamos la hipótesis nula, lo que indica que las varianzas de los grupos son significativamente diferentes.

4.2.2 Realizar la prueba t de Student

Ya que los supuestos se cumplen:

# Realizar la prueba t para comparar la masa corporal entre sexos
t_test_result <- t.test(body_mass_g ~ sex, 
                        data = peng_adelie, 
                        var.equal = TRUE)
t_test_result

## 
##  Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -13.126, df = 144, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -776.2484 -573.0666
## sample estimates:
## mean in group female   mean in group male 
##             3368.836             4043.493

Interpretación de los Resultados El valor p (p-value) indica si existe una diferencia significativa entre las medias de male y female.

Dado que p < 0.05, rechazamos la hipótesis nula y concluimos que hay una diferencia significativa en la masa corporal entre los sexos.

4.3 Visualización de los Resultados

Finalmente, visualizamos la comparación de las medias utilizando ggplot2:

# Visualizar la comparación de masa corporal por sexo
ggplot(peng_adelie, aes(x = sex, y = body_mass_g, fill = sex)) +
  geom_boxplot(outlier.colour = "red", alpha= 0.5) +
  labs(title = "Comparación de Masa Corporal entre Sexos",
       subtitle = "Pingüinos de Adelia",
       x = "Sexo", y = "Masa Corporal (g)") +
  scale_fill_manual(values = c("skyblue2", "pink2")) +
  theme_minimal()

El gráfico muestra claramente que hay una diferencia significativa en la masa corporal entre machos y hembras en la especie Adelie, con los machos siendo, en promedio, más pesados. Esto se alinea con los resultados de la prueba t de Student para comparar la masa corporal entre sexos.

Tambien podemos agregar los resultados del t-test al gráfico:

# Extraer el valor p
p_value <- t_test_result$p.value
p_value_text <- ifelse(p_value < 0.001, "p < 0.001", paste("p =", round(p_value, 3)))

# Gráfico del boxplot con anotación del valor p
ggplot(peng_adelie, aes(x = sex, y = body_mass_g, fill = sex)) +
  geom_boxplot(outlier.colour = "red", alpha= 0.5) +
  labs(title = "Comparación de Masa Corporal entre Sexos",
       subtitle = "Pingüinos de Adelia",
       x = "Sexo", y = "Masa Corporal (g)") +
  scale_fill_manual(values = c("skyblue2", "pink2")) +
  theme_minimal()+
  annotate("text", x = 1.5, 
           y = max(peng_adelie$body_mass_g, na.rm = TRUE)+50, 
           label = paste("t-test:", p_value_text), size = 4, color = "black")

5 Evaluación de Relaciones entre Variables Cuantitativas: Regresión Lineal

La regresión lineal es un método estadístico que nos permite modelar la relación entre dos variables cuantitativas. En este caso, utilizaremos la regresión lineal para evaluar cómo la longitud de la aleta flipper_length_mm influye en la masa corporal body_mass_g de las tres especies de pingüinos del archipielago de Palmer.

5.1 ¿Qué es la Regresión Lineal?

La regresión lineal simple es un modelo que permite predecir el valor de una variable dependiente (Y) a partir de una variable independiente (X). En nuestro caso:

Variable independiente (X): Longitud de la aleta (flipper_length_mm)
Variable dependiente (Y): Masa corporal (body_mass_g)

El modelo de regresión lineal se expresa como:

\[ Y = \beta_0 + \beta_1X + \epsilon \]

Donde:

\(\beta_0\): Intercepto del modelo (valor de Y cuando X = 0).
\(\beta_1\): Pendiente del modelo (cambio en Y por cada unidad de cambio en X).
\(\epsilon\): Término de error.

5.2 Visualización de la Relación entre las Variables

Primero, visualicemos la relación entre la longitud de la aleta y la masa corporal utilizando ggplot2:

# Visualizar la relación entre la longitud de la aleta y la masa corporal
peng %>% 
  select(body_mass_g, flipper_length_mm) %>% 
  ggplot(aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point() +
  geom_smooth(method = lm, se = TRUE, color = "blue", fill = "lightblue") +
  labs(title = "Relación entre Longitud de la Aleta y Masa Corporal",
       x = "Longitud de la Aleta (mm)", 
       y = "Masa Corporal (g)") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

geom_point(): Muestra los puntos de datos.
geom_smooth(method = lm, se = TRUE): Añade la línea de regresión y el intervalo de confianza alrededor de la línea.

5.3 Modelo de Regresión Lineal

5.3.1 Ajuste del Modelo de Regresión Lineal

Ajustemos el modelo de regresión lineal en R usando lm():

# Ajustar el modelo de regresión lineal
modelo_regresion <- lm(body_mass_g ~ flipper_length_mm, data = peng)

# Resumen del modelo
resultados_lm <- summary(modelo_regresion)

5.3.2 Interpretación de Resultados

resultados_lm

## 
## Call:
## lm(formula = body_mass_g ~ flipper_length_mm, data = peng)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1057.33  -259.79   -12.24   242.97  1293.89 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5872.09     310.29  -18.93   <2e-16 ***
## flipper_length_mm    50.15       1.54   32.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 393.3 on 331 degrees of freedom
## Multiple R-squared:  0.7621, Adjusted R-squared:  0.7614 
## F-statistic:  1060 on 1 and 331 DF,  p-value: < 2.2e-16

Interpretación de los Resultados Regresión Lineal
El resumen del modelo incluye información importante:

Intercepto (β₀): Representa el valor esperado de la masa corporal cuando la longitud de la aleta es 0 (aunque este valor puede no ser relevante en un contexto realista).
Pendiente (β₁): Indica cuánto cambia la masa corporal por cada unidad de cambio en la longitud de la aleta. Si β₁ es positivo, hay una relación positiva; si es negativo, hay una relación negativa.
Valor p: Indica si la relación entre las variables es estadísticamente significativa.
R²: Mide qué proporción de la variación en la masa corporal puede explicarse por la longitud de la aleta.

Tambien podemos ingresar a la lista de resultados indexando:

# Coeficientes
resultados_lm$coefficients

##                      Estimate Std. Error   t value      Pr(>|t|)
## (Intercept)       -5872.09268 310.285155 -18.92483  1.183941e-54
## flipper_length_mm    50.15327   1.540231  32.56217 3.132836e-105

# R2
resultados_lm$r.squared

## [1] 0.7620922

5.3.3 Verificación de Supuestos del Modelo

Linealidad: existe una relación lineal entre la variable independiente (X) y la variable dependiente (Y). Ya observamos en el punto 5.2 la relación lineal entre body_mass_g y flipper_length_mm.
Independencia de los Residuos: Los residuos (las diferencias entre los valores observados y los valores predichos) deben ser independientes entre sí. No deben mostrar ningún patrón o dependencia en sus valores. La falta de independencia indica que los errores en las predicciones no son aleatorios.
Normalidad de los Residuos: Los residuos deben seguir una distribución normal. Esto significa que los errores se distribuyen de manera simétrica alrededor de la línea de regresión. Si los residuos no son normales, los intervalos y pruebas pueden no ser precisos.
Homocedasticidad de los Residuos: La varianza de los residuos debe ser constante a lo largo de todos los valores de X. Si hay heterocedasticidad (es decir, la varianza de los residuos no es constante), las estimaciones del modelo pueden ser ineficientes y los intervalos de confianza y pruebas de hipótesis podrían no ser válidos.

Mediante la funcion autoplot() del paquete ggfortify podemos explorar los últimos 3 supuestos:

autoplot(modelo_regresion)

Gráfico	Qué muestra	Si es correcto	Si hay problemas
1. Residuals vs Fitted	Muestra los residuos frente a los valores ajustados.	Los residuos se dispersan aleatoriamente alrededor de la línea horizontal en 0.	Si hay un patrón claro (curva o abanico), indica problemas de linealidad o heterocedasticidad.
2. Normal Q-Q	Compara la distribución de los residuos con una distribución normal.	Los puntos se alinean a lo largo de la línea diagonal.	Los puntos se desvían notablemente de la línea, indicando que los residuos no son normales.
3. Scale-Location	Muestra la raíz cuadrada de los residuos estandarizados frente a los valores ajustados.	Los puntos se distribuyen uniformemente sin mostrar un patrón claro.	Un patrón de abanico indica que la varianza de los residuos no es constante (heterocedasticidad).
4. Residuals vs Leverage	Muestra la relación entre residuos estandarizados y la influencia de cada observación (leverage).	Identifica valores influyentes; pocos puntos fuera de las líneas punteadas de “Cook’s distance”.	Puntos fuera de las líneas de “Cook’s distance” indican observaciones influyentes.

Otra manera de identificar observaciones que influyen en mi modelo es con un gráfico de distancias de Cook. Para ello crearemos un data frame con la funcion fortify().

# Fortificar el modelo de regresión
fortified_model <- fortify(modelo_regresion) %>% 
  mutate(Observation = row_number())  # Agregar el número de observación

# Crear el gráfico de distancias de Cook
  ggplot(fortified_model, aes(x = Observation, y = .cooksd)) +
  geom_bar(stat = "identity", fill = "skyblue", color = "black") +
  geom_hline(yintercept = 4/(nrow(fortified_model)-3), 
             linetype = "dashed", color = "red") +
  labs(title = "Distancias de Cook",
       x = "Observación",
       y = "Distancia de Cook") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Las distancias de Cook nos permiten identificar observaciones que influyen significativamente en el modelo de regresión. Un valor alto indica que una observación tiene un gran impacto en los resultados y podría distorsionar el modelo.

5.4 Predicciones con el modelo y gráfico final

Realizaremos predicciones con la función predict():

# Crear nuevos datos para predicción
nuevos_datos <- data.frame(flipper_length_mm = seq(min(peng$flipper_length_mm), 
    max(peng$flipper_length_mm), 
    length.out = 100))

# Realizar la predicción con intervalos de confianza y predicción
predicciones <- predict(modelo_regresion, newdata = nuevos_datos, 
                        interval = "prediction")

# Agregar las predicciones al data frame
nuevos_datos <- cbind(nuevos_datos, bm = predicciones)

nuevos_datos

Utilizando la función predict(), hemos generado predicciones para la masa corporal (body_mass_g) basadas en diferentes valores de la longitud de la aleta (flipper_length_mm).

El gráfico a continuación muestra los datos observados, la línea de regresión ajustada (en color azul) y los intervalos de predicción (área sombreada en azul claro) que indican el rango donde se espera que se encuentren futuras observaciones para un valor dado de flipper_length_mm.

# Obtener la ecuación del modelo
intercept <- coef(modelo_regresion)[1]
slope <- coef(modelo_regresion)[2]
equation <- paste("y = ", round(intercept, 2), " + ", round(slope, 2),
                    "x", sep = "")


# Gráfico del modelo de regresión con intervalos de confianza y predicción
ggplot(peng, aes(x = flipper_length_mm)) +
  geom_point(alpha = 0.6, aes(y = body_mass_g)) +
  geom_line(data = nuevos_datos, aes(y = bm.fit), 
            color = "navy") +
  geom_ribbon(data = nuevos_datos, aes(ymin = bm.lwr, ymax = bm.upr), 
              alpha = 0.2, fill = "lightblue") +
  labs(title = "Modelo de Regresión Lineal",
       subtitle = equation,
       x = "Longitud de la Aleta (mm)", 
       y = "Masa Corporal (g)") +
  theme_minimal()

Esto nos permite evaluar cómo el modelo ajustado describe la relación entre estas dos variables y la precisión de las predicciones que realiza.

Programación y Análisis de Datos en Dendrocronología y Ciencias Ambientales

Clase 4: Aplicaciones estadísticas de R

Nicolás Carrasco Ilabaca - Laboratorio de Dendrocronología y Cambio Global