Correlación y regresión lineal

Gráfico de dispersión

Los diagramas de dispersión son las herramientas más comunes y efectivas para visualizar la relación entre dos variables numéricas.

El conjunto de datos ncbirths es una muestra aleatoria de 1,000 casos tomados de un conjunto de datos más grande recolectado en 2004. Cada caso describe el nacimiento de un solo niño nacido en Carolina del Norte, junto con varias características del niño (por ejemplo, peso al nacer, duración de la gestación, etc. .), la madre del niño (por ejemplo, la edad, el peso ganado durante el embarazo, los hábitos de fumar, etc.) y el padre del niño (por ejemplo, la edad). Puede ver el archivo de ayuda para estos datos ejecutando? Ncbirths en la consola.

Usando el conjunto de datos ncbirths, haga un diagrama de dispersión usando ggplot () para ilustrar cómo el peso al nacer de estos bebés varía según el número de semanas de gestación.

library(ggplot2)

library(openintro)

## Please visit openintro.org for free statistics materials

## 
## Attaching package: 'openintro'

## The following object is masked from 'package:ggplot2':
## 
##     diamonds

## The following objects are masked from 'package:datasets':
## 
##     cars, trees

# Scatterplot of weight vs. weeks
ggplot(data = ncbirths, aes(y = weight, x = weeks)) +
geom_point()

## Warning: Removed 2 rows containing missing values (geom_point).

Diagramas de caja como diagramas de dispersión discretos / condicionados

Si es útil, puede pensar en los diagramas de caja como diagramas de dispersión para los que se ha discretizado la variable en el eje x.

La función cut () toma dos argumentos: la variable continua que desea discretizar y el número de saltos que desea realizar en esa variable continua para discretizarla.

Usando el conjunto de datos ncbirths nuevamente, haga un diagrama de caja que ilustre cómo el peso al nacer de estos bebés varía según el número de semanas de gestación. Esta vez, use la función cut () para discretizar la variable x en seis intervalos (es decir, cinco descansos).

# Boxplot of weight vs. weeks
ggplot(data = ncbirths, 
       aes(x = cut(weeks, breaks = 5), y = weight)) + 
  geom_boxplot()

Crear diagramas de dispersión

Crear diagramas de dispersión es simple y son tan útiles que vale la pena exponerse a muchos ejemplos. Con el tiempo, se familiarizará con los tipos de patrones que ve. Comenzará a reconocer cómo los diagramas de dispersión pueden revelar la naturaleza de la relación entre dos variables.

En este ejercicio, y a lo largo de este capítulo, utilizaremos varios conjuntos de datos que se enumeran a continuación. Estos datos están disponibles a través del paquete openintro. Brevemente:

El conjunto de datos de mamíferos contiene información sobre 39 especies diferentes de mamíferos, incluido su peso corporal, peso cerebral, tiempo de gestación y algunas otras variables.
El conjunto de datos mlbBat10 contiene estadísticas de bateo para 1,199 jugadores de Grandes Ligas durante la temporada 2010.
El conjunto de datos bdims contiene mediciones de circunferencia corporal y diámetro esquelético para 507 individuos físicamente activos.
El conjunto de datos sobre fumar contiene información sobre los hábitos de fumar de 1.691 ciudadanos del Reino Unido.

Para ver documentación más exhaustiva, use el? o help() y dentro la función.

Utilizando el conjunto de datos de mamíferos, cree un diagrama de dispersión que ilustre cómo varía el peso del cerebro de un mamífero en función de su peso corporal.

# Mammals scatterplot
ggplot(data = mammals, aes(y = BrainWt, x = BodyWt)) +
  geom_point()

Usando el conjunto de datos mlbBat10, cree un diagrama de dispersión que ilustre cómo el porcentaje de slugging (SLG) de un jugador varía en función de su porcentaje en base (OBP).

# Baseball player scatterplot
ggplot(data = mlbBat10, aes(y = SLG, x = OBP)) +
  geom_point()

Usando el conjunto de datos bdims, cree un diagrama de dispersión que ilustre cómo varía el peso de una persona en función de su altura. Use el color para separar por sexo, que necesitará forzar a un factor con factor ().

# Body dimensions scatterplot
ggplot(data = bdims, aes(y = wgt, x = hgt, color = factor(sex))) +
  geom_point()

Utilizando el conjunto de datos sobre fumar, cree un diagrama de dispersión que ilustre cómo varía la cantidad que fuma una persona de lunes a viernes en función de su edad.

# Smoking scatterplot
ggplot(data = smoking, aes(y = amtWeekdays, x = age)) +
  geom_point()

## Warning: Removed 1270 rows containing missing values (geom_point).

Caracterización de diagramas de dispersión

Este diagrama de dispersión muestra la relación entre las tasas de pobreza y las tasas de graduación de secundaria de los condados en los Estados Unidos.

Transformaciones

La relación entre dos variables puede no ser lineal. En estos casos, a veces podemos ver patrones extraños e incluso inescrutables en un diagrama de dispersión de los datos. A veces realmente no hay una relación significativa entre las dos variables. Otras veces, una transformación cuidadosa de una o ambas variables puede revelar una relación clara.

Recuerde el extraño patrón que vio en el diagrama de dispersión entre el peso del cerebro y el peso corporal entre los mamíferos en un ejercicio anterior. ¿Podemos usar transformaciones para aclarar esta relación?

ggplot2 proporciona varios mecanismos diferentes para ver relaciones transformadas. La función coord_trans () transforma las coordenadas de la trama. Alternativamente, las funciones scale_x_log10 () y scale_y_log10 () realizan una transformación de registro de base 10 de cada eje. Tenga en cuenta las diferencias en la apariencia de los ejes.

El conjunto de datos de mamíferos está disponible en su espacio de trabajo.

Use coord_trans () para crear un diagrama de dispersión que muestre cómo varía el peso del cerebro de un mamífero en función de su peso corporal, donde los ejes x e y están en una escala "log10".

# Scatterplot with coord_trans()
ggplot(data = mammals, aes(x = BodyWt, y = BrainWt)) +
  geom_point() + 
  coord_trans(x = "log10", y = "log10")

Use scale_x_log10() and scale_y_log10() para ver el mismo efecto que en el gráfico anterior y ver fácilmente la asociaciòn existrente entre los pesos de los cerebros y del cuerpo al nacer

# Scatterplot with scale_x_log10() and scale_y_log10()
ggplot(data = mammals, aes(x = BodyWt, y = BrainWt)) +
  geom_point() +
  scale_x_log10() + 
  scale_y_log10()

dentificar valores atípicos

En el Capítulo 5, discutiremos cómo los valores atípicos pueden afectar los resultados de un modelo de regresión lineal y cómo podemos tratarlos. Por ahora, es suficiente simplemente identificarlos y observar cómo la relación entre dos variables puede cambiar como resultado de la eliminación de valores atípicos.

Recuerde que en el ejemplo del béisbol anterior en el capítulo, la mayoría de los puntos se agruparon en la esquina inferior izquierda de la trama, lo que dificulta ver el patrón general de la mayoría de los datos. Esta dificultad fue causada por unos pocos jugadores periféricos cuyos porcentajes en base (OBP) eran excepcionalmente altos. Estos valores están presentes en nuestro conjunto de datos solo porque estos jugadores tenían muy pocas oportunidades de bateo.

Tanto OBP como SLG se conocen como estadísticas de frecuencia, ya que miden la frecuencia de ciertos eventos (en oposición a su recuento). Para comparar estas tasas de manera sensata, tiene sentido incluir solo jugadores con un número razonable de oportunidades, de modo que estas tasas observadas tengan la oportunidad de acercarse a sus frecuencias a largo plazo.

En la Major League Baseball, los bateadores califican para el título de bateo solo si tienen 3.1 apariciones en el plato por juego. Esto se traduce en aproximadamente 502 apariciones en el plato en una temporada de 162 juegos. El conjunto de datos mlbBat10 no incluye las apariencias de placa como una variable, pero podemos usar at-bats (AB), que constituyen un subconjunto de apariencias de placa, como proxy.

Use filter () para mantener solo a los jugadores que tuvieron al menos 200 turnos al bate, asignando a ab_gt_200.
Usando ab_gt_200, cree un diagrama de dispersión para SLG en función de OBP.
Encuentra la fila de ab_gt_200 correspondiente al jugador (con al menos 200 turnos al bate) cuyo OBP estaba por debajo de 0.200.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Filter for AB greater than or equal to 200
ab_gt_200 <- mlbBat10 %>%
  filter(AB >= 200) 

# Scatterplot of SLG vs. OBP
ggplot(ab_gt_200, aes(x = OBP, y = SLG)) +
  geom_point()

# Identify the outlying player
ab_gt_200 %>%
  filter(OBP < 0.200)

Covarianza calculada

La función cor (x, y) calculará la correlación producto-momento de Pearson entre las variables, x e y. Como esta cantidad es simétrica con respecto a x e y, no importa en qué orden coloque las variables.

Al mismo tiempo, la función cor () es muy conservadora cuando encuentra datos faltantes (por ejemplo, NA). El argumento de uso le permite anular el comportamiento predeterminado de devolver NA siempre que cualquiera de los valores encontrados sea NA. Establecer el argumento de uso en “pairwise.complete.obs” permite que cor () calcule el coeficiente de correlación para aquellas observaciones en las que no faltan los valores de x e y.

\[-\infty<Cov(X,Y)=\frac{\sum_{i=1}^{N}(x_i-\mu_x)(y_i-\mu_y)}{N}<+\infty\]

Un par de variables asociadas positivamente

P <- c(70, 60)
E <- c(170, 160)
mean(P)

## [1] 65

mean(E)

## [1] 165

\[-\infty<Cov(X,Y)=\frac{(P_1-\mu_P)(E_1-\mu_E)+(P_2-\mu_P)(E_2-\mu_E)}{2}=\frac{(70-65)(170-165)+(60-65)(160-165)}{2}<+\infty\]

\[-\infty<Cov(X,Y)=\frac{(70-65)(170-165)+(60-65)(160-165)}{2}=\frac{5\cdot5+(-5)\cdot(-5)}{2}=\frac{25+25}{2}=25<+\infty\]

((P[1]-mean(P))*(E[1]-mean(E))+(P[2]-mean(P))*(E[2]-mean(E)))/2

## [1] 25

P <- c(70,60)
E <- c(1.70,1.60)
((P[1]-mean(P))*(E[1]-mean(E))+(P[2]-mean(P))*(E[2]-mean(E)))/2

## [1] 0.25

Un par de variables asociadas negativamente

CigarrillosDiarios <- c(5,17)
EsperanzadeVidaA <- c(85,63)
((CigarrillosDiarios[1]-mean(CigarrillosDiarios))*(EsperanzadeVidaA[1]-mean(EsperanzadeVidaA))+(CigarrillosDiarios[2]-mean(CigarrillosDiarios))*(EsperanzadeVidaA[2]-mean(EsperanzadeVidaA)))/2

## [1] -66

CigarrillosSemanales <- 7*c(5,17)
EsperanzadeVidaA <- 12*c(85,63)
((CigarrillosSemanales[1]-mean(CigarrillosSemanales))*(EsperanzadeVidaA[1]-mean(EsperanzadeVidaA))+(CigarrillosSemanales[2]-mean(CigarrillosSemanales))*(EsperanzadeVidaA[2]-mean(EsperanzadeVidaA)))/2

## [1] -5544

ncbirths %>%
  summarize(N = n(), covar = cov(weight, mage))

ncbirths %>%
  summarize(N = n(), covar = cov(weight, weeks, use = "pairwise.complete.obs"))

Correlación calculada

\[-1<\rho=\frac{Cov(X,Y)}{\sqrt{V(X)}\sqrt{V(Y)}}=\frac{Cov(X,Y)}{sd_{X}sd_{Y}}=\frac{\frac{\sum_{i=1}^{N}(x_i-\mu_x)(y_i-\mu_y)}{N}}{\sqrt{\frac{\sum_{i=1}^{N}(x_i-\mu_x)}{N}}\sqrt{\frac{\sum_{i=1}^{N}(y_i-\mu_y)}{N}}}<+1\]

Use cor() para calcular la correlación entre el peso al nacer de los bebés en el conjunto de datos ncbirths y la edad de su madre. No faltan datos en ninguna de las variables.

Calcule la correlación entre el peso al nacer y el número de semanas de gestación para todos los pares que no faltan.

ncbirths %>%
  summarize(N = n(), r = cor(weight, mage))

ncbirths %>%
  summarize(N = n(), r = cor(weight, weeks, use = "pairwise.complete.obs"))

Explorando Anscombe

En 1973, Francis Anscombe creó cuatro conjuntos de datos con propiedades numéricas notablemente similares, pero obviamente diferentes relaciones gráficas. El conjunto de datos Anscombe contiene las coordenadas xey para estos cuatro conjuntos de datos, junto con una variable de agrupación, conjunto, que distingue al cuarteto.

Puede ser útil recordar la relación gráfica al ver los cuatro diagramas de dispersión:

Para cada uno de los cuatro conjuntos de puntos de datos en el conjunto de datos Anscombe, calcule lo siguiente en el orden especificado. Los nombres se proporcionan en su llamada para resumir ().

Número de observaciones, N
Media de x
Desviación estándar de x
Media de y
Desviación estándar de y
Coeficiente de correlación entre x e y

library(datasets)
library(tidyr)

Anscombe <- anscombe %>% 
  pivot_longer(
    everything(),
    names_to = c(".value", "set"),
    names_pattern = "(.)(.)"
  ) %>% 
  as_tibble()

Anscombe %>%
  summarize(
    N = n(), 
    mean_of_x = mean(x), 
    std_dev_of_x = sd(x), 
    mean_of_y = mean(y), 
    std_dev_of_y = sd(y),
    correlation_between_x_and_y = cor(x, y),
    groups = set
  )

op <- par(no.readonly = TRUE)
xls <- c(2, 19)
yls <- c(2, 14)
par(oma = c(2, 2, 0, 2),
    mar = c(3, 3, 1, 0),
    mfrow = c(2, 2),
    pch = 16
    )

plot(y1 ~ x1, data = anscombe, xlim = xls, ylim = yls, xlab = "", ylab = "", col = "red")
plot(y2 ~ x2, data = anscombe, xlim = xls, ylim = yls, xlab = "", ylab = "", col = "blue")
plot(y3 ~ x3, data = anscombe, xlim = xls, ylim = yls, xlab = "", ylab = "", col = "green")
plot(y4 ~ x4, data = anscombe, xlim = xls, ylim = yls, xlab = "", ylab = "", col = "black")

mtext(text = "X axis", side = 1, line = 0, outer = TRUE)
mtext(text = "Y axis", side = 2, line = 0, outer = TRUE)

library(Tmisc)
data(quartet)
str(quartet)

## 'data.frame':    44 obs. of  3 variables:
##  $ set: Factor w/ 4 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ x  : int  10 8 13 9 11 14 6 4 12 7 ...
##  $ y  : num  8.04 6.95 7.58 8.81 8.33 ...

library(dplyr)
quartet %>% 
  group_by(set) %>% 
  summarize(mean(x), median(x), sd(x), mean(y), median(y), sd(y), cor(x, y), .groups = "drop")

library(ggplot2)
ggplot(quartet, aes(x, y)) + geom_point() + geom_smooth(method = lm, se = FALSE) + geom_rug() + 
    facet_wrap(~set)

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

Estimar el valor del coeficiente de correlación entre dos cantidades de su diagrama de dispersión puede ser complicado. Los estadísticos han demostrado que la percepción de las personas sobre la fuerza de estas relaciones puede verse influenciada por elecciones de diseño como las escalas x e y.

Sin embargo, con un poco de práctica mejorará su percepción de correlación. Alterne entre los cuatro diagramas de dispersión en la ventana de trazado, cada uno de los cuales ha visto en un ejercicio anterior. Anote su mejor estimación del valor del coeficiente de correlación entre cada par de variables. Luego, compare estos valores con los valores reales que calcula en este ejercicio.

Si tiene problemas para recordar nombres de variables, puede ser útil obtener una vista previa de un conjunto de datos en la consola con str () o glimpse ().

Dibuje la gráfica y luego calcule la correlación entre OBP y SLG para todos los jugadores en el conjunto de datos mlbBat10.

ggplot(data = mlbBat10, aes(x = OBP, y = SLG)) +
  geom_smooth(method = lm, se = FALSE) +
  geom_point()

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

mlbBat10 %>%
  summarize(N = n(), r = cor(OBP, SLG))

Dibuje la gráfica y luego calcule la correlación entre OBP y SLG para todos los jugadores en el conjunto de datos mlbBat10 con al menos 200 turnos al bate.

mlbBat10 %>% 
    filter(AB > 200) %>%
    ggplot(aes(x = OBP, y = SLG)) +
    geom_smooth(method = lm, se = FALSE) + 
    geom_point()

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

mlbBat10 %>%
  filter(AB >= 200) %>%
  summarize(N = n(), r = cor(OBP, SLG))

Dibuje la gráfica y luego calcule la correlación entre altura y peso para cada sexo en el conjunto de datos bdims.

ggplot(data = bdims, aes(x = hgt, y = wgt, color = factor(sex))) +
  geom_point() 

bdims %>%
  group_by(sex) %>%
  summarize(N = n(), r = cor(hgt, wgt))

## `summarise()` ungrouping output (override with `.groups` argument)

Dibuje la gráfica y luego calcule la correlación entre el peso corporal y el peso del cerebro para todas las especies de mamíferos. Junto con este cálculo, calcule la correlación entre las mismas dos cantidades después de tomar sus logaritmos naturales.

# Run this and look at the plot
ggplot(data = mammals, aes(x = BodyWt, y = BrainWt)) +
  geom_point() + scale_x_log10() + scale_y_log10()

# Correlation among mammals, with and without log
mammals %>%
  summarize(N = n(), 
            r = cor(BodyWt, BrainWt), 
            r_log = cor(log(BodyWt), log(BrainWt)))

Correlación espuria en datos aleatorios

Los estadísticos siempre deben ser escépticos de las correlaciones potencialmente espurias. Los seres humanos son muy buenos para ver patrones en los datos, a veces cuando los patrones en sí mismos son en realidad solo ruido aleatorio. Para ilustrar cuán fácil puede ser caer en esta trampa, buscaremos patrones en datos verdaderamente aleatorios.

El conjunto de datos de ruido contiene 20 conjuntos de variables x e y extraídas al azar de una distribución normal estándar. Cada conjunto, denotado como z, tiene 50 observaciones de pares x, y. ¿Ves algún par de variables que puedan estar significativamente correlacionadas? ¿Todos los coeficientes de correlación son cercanos a cero?

Cree un diagrama de dispersión facetado que muestre la relación entre cada uno de los 20 conjuntos de pares de variables aleatorias x e y. Necesitará la función facet_wrap () para esto.
Calcule la correlación real entre cada uno de los 20 conjuntos de pares de x e y.
Identifique los conjuntos de datos que muestran una correlación no trivial de más de 0.2 en valor absoluto.

noise <- data.frame(x = rnorm(1000), y = rnorm(1000), z = rep(seq(from = 1, to = 20, by = 1), len = 100))

# Create faceted scatterplot
ggplot(data = noise, aes(x = x, y = y)) +
  geom_point() +
  facet_wrap(~ z)

# Compute correlations for each dataset
noise_summary <- noise %>%
  group_by(z) %>%
  summarize(N = n(), spurious_cor = cor(x, y))

## `summarise()` ungrouping output (override with `.groups` argument)

# Isolate sets with correlations above 0.2 in absolute strength
noise_summary %>%
  filter(abs(spurious_cor) > 0.2)

La línea de “mejor ajuste”

El modelo de regresión lineal simple para una respuesta numérica en función de una variable explicativa numérica se puede visualizar en el diagrama de dispersión correspondiente mediante una línea recta. Esta es una línea de “mejor ajuste” que corta los datos de una manera que minimiza la distancia entre la línea y los puntos de datos.

Podríamos considerar la regresión lineal como un ejemplo específico de una clase más grande de modelos suaves. La función geom_smooth () le permite dibujar dichos modelos sobre un diagrama de dispersión de los datos en sí. Esta técnica se conoce como visualizar el modelo en el espacio de datos. El argumento del método para geom_smooth () le permite especificar qué clase de modelo liso desea ver. Como estamos explorando modelos lineales, estableceremos este argumento en el valor “lm”.

Tenga en cuenta que geom_smooth () también toma un argumento se que controla el error estándar, que ignoraremos por ahora.

Cree un diagrama de dispersión del peso corporal en función de la altura para todas las personas en el conjunto de datos bdims con un modelo lineal simple trazado sobre los datos.

# Scatterplot with regression line
ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = TRUE)

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

Singularidad de la línea de regresión de mínimos cuadrados

El criterio de mínimos cuadrados implica que la pendiente de la línea de regresión es única. En la práctica, la pendiente es calculada por R. En este ejercicio, experimentará tratando de encontrar el valor óptimo para la pendiente de regresión para el peso en función de la altura en el conjunto de datos bdims a través de prueba y error.

Para ayudar, hemos creado una función personalizada para usted llamada add_line (), que toma un solo argumento: el coeficiente de pendiente propuesto.

El conjunto de datos bdims está disponible en su espacio de trabajo. Experimente con diferentes valores (al entero más cercano) del parámetro my_slope hasta que encuentre uno que le parezca mejor.

add_line <- function (my_slope){

  bdims_summary <- bdims %>%
    summarize(N = n(), r = cor(hgt, wgt),
              mean_hgt = mean(hgt), mean_wgt = mean(wgt),
              sd_hgt = sd(hgt), sd_wgt = sd(wgt)) %>%
    mutate(true_slope = r * sd_wgt / sd_hgt, 
           true_intercept = mean_wgt - true_slope * mean_hgt)
  p <- ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
    geom_point() + 
    geom_point(data = bdims_summary, 
               aes(x = mean_hgt, y = mean_wgt), 
               color = "red", size = 3)
  
  my_data <- bdims_summary %>%
    mutate(my_slope = my_slope, 
           my_intercept = mean_wgt - my_slope * mean_hgt)
  p + geom_abline(data = my_data, 
                  aes(intercept = my_intercept, slope = my_slope), color = "dodgerblue")
}

# Estimate optimal value of my_slope
add_line(my_slope = 1.0)

Ajuste de un modelo lineal “a mano”

Recordemos el modelo de regresión lineal simple:

\[Y = b_0 + b_1{\cdot}X\]

Dos hechos le permiten calcular la pendiente b1 e interceptar b0

de un modelo de regresión lineal simple de algunas estadísticas de resumen básicas.

Primero, la pendiente se puede definir como:

\[b1 = r(X,Y)⋅sYsX\]

donde rX, Y representa la correlación (cor ()) de X e Y y sX y sY representan la desviación estándar (sd ()) de X e Y

, respectivamente.

Segundo, el punto (x¯, y¯) siempre está en la línea de regresión de mínimos cuadrados, donde x¯ e y¯ denotan el promedio de x e y

, respectivamente.

El marco de datos bdims_summary contiene toda la información que necesita para calcular la pendiente e intercepción de la línea de regresión de mínimos cuadrados para el peso corporal (Y ) en función de la altura (X). Es posible que debas hacer algo de álgebra para resolver b0

Imprima el marco de datos bdims_summary.
Use mutate () para agregar la pendiente e interceptar al marco de datos bdims_summary.

bdims_summary <- bdims %>%
  summarize(N = n(), r = cor(hgt, wgt),
            mean_hgt = mean(hgt), mean_wgt = mean(wgt),
            sd_hgt = sd(hgt), sd_wgt = sd(wgt))

# Print bdims_summary
bdims_summary

# Add slope and intercept
bdims_summary %>%
  mutate(slope = r * sd_wgt / sd_hgt, 
         intercept = mean_wgt - r * sd_wgt / sd_hgt * mean_hgt)

Regresión a la media

La regresión a la media es un concepto atribuido a Sir Francis Galton. La idea básica es que las observaciones aleatorias extremas tenderán a ser menos extremas en una segunda prueba. Esto se debe simplemente a la casualidad. Si bien “regresión a la media” y “regresión lineal” no son lo mismo, las examinaremos juntas en este ejercicio.

Una forma de ver los efectos de la regresión a la media es comparar las alturas de los padres con las de sus hijos. Si bien es cierto que las madres y los padres altos tienden a tener hijos altos, esos niños tienden a ser menos altos que sus padres, en relación con el promedio. Es decir, los padres que son 3 pulgadas más altos que el padre promedio tienden a tener hijos que pueden ser más altos que el promedio, pero de menos de 3 pulgadas.

Los conjuntos de datos Galton_men y Galton_women contienen datos recopilados originalmente por el mismo Galton en la década de 1880 sobre las alturas de hombres y mujeres, respectivamente, junto con las alturas de sus padres.

Compare la pendiente de la línea de regresión con la pendiente de la línea diagonal. ¿Qué te dice esto?

Cree un diagrama de dispersión de la altura de los hombres en función de la altura de su padre. Agregue la línea de regresión lineal simple y una línea diagonal (con pendiente igual a 1 e intersección igual a 0) a la gráfica.
Cree un diagrama de dispersión de la altura de las mujeres en función de la altura de su madre. Agregue la línea de regresión lineal simple y una línea diagonal a la gráfica.

library(HistData)

GaltonFamilies

Galton_men <- GaltonFamilies %>%
  filter(gender == "male") %>%
  mutate(sex = case_when(gender == "male" ~ "M"), nkids = children, height = childHeight) %>%
  select(family, father, mother, sex, height, nkids)

Galton_women <- GaltonFamilies %>%
  filter(gender == "female") %>%
  mutate(sex = case_when(gender == "female" ~ "W"), nkids = children, height = childHeight) %>%
  select(family, father, mother, sex, height, nkids)

# Height of children vs. height of father
ggplot(data = Galton_men, aes(x = father, y = height)) +
  geom_point() + 
  geom_abline(slope = 1, intercept = 0) + 
  geom_smooth(method = "lm", se = FALSE)

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

# Height of children vs. height of mother
ggplot(data = Galton_women, aes(x = mother, y = height)) +
  geom_point() + 
  geom_abline(slope = 1, intercept = 0) + 
  geom_smooth(method = "lm", se = FALSE)

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

¡Excelente! Debido a la pendiente de la línea de regresión es menor que 1 (pendiente de la línea diagonal) tanto para hombres como para mujeres, ¡podemos verificar la regresión de Sir Francis Galton al concepto medio!

Ajuste de modelos lineales simples

Si bien la función geom_smooth (method = “lm”) es útil para dibujar modelos lineales en un diagrama de dispersión, en realidad no devuelve las características del modelo. Sin embargo, como sugiere esa sintaxis, la función que crea modelos lineales es lm (). Esta función generalmente toma dos argumentos:

Una fórmula que especifica el modelo.
Un argumento de datos para el marco de datos que contiene los datos que desea utilizar para ajustarse al modelo.

La función lm () devuelve un objeto modelo que tiene la clase “lm”. Este objeto contiene mucha información sobre su modelo de regresión, incluidos los datos utilizados para ajustar el modelo, la especificación del modelo, los valores ajustados y los residuos, etc.

Usando el conjunto de datos bdims, cree un modelo lineal para el peso de las personas en función de su altura.
Usando el conjunto de datos mlbBat10, cree un modelo lineal para SLG en función de OBP.
Usando el conjunto de datos de mamíferos, cree un modelo lineal para el peso corporal de los mamíferos en función de su peso cerebral, después de tomar el registro natural de ambas variables.

# Linear model for weight as a function of height
lm(wgt ~ hgt, data = bdims)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Coefficients:
## (Intercept)          hgt  
##    -105.011        1.018

# Linear model for SLG as a function of OBP
lm(SLG ~ OBP, data = mlbBat10)

## 
## Call:
## lm(formula = SLG ~ OBP, data = mlbBat10)
## 
## Coefficients:
## (Intercept)          OBP  
##    0.009407     1.110323

# Log-linear model for body weight as a function of brain weight
lm(log(BodyWt) ~ log(BrainWt), data = mammals)

## 
## Call:
## lm(formula = log(BodyWt) ~ log(BrainWt), data = mammals)
## 
## Coefficients:
##  (Intercept)  log(BrainWt)  
##       -2.509         1.225

mod <- lm(wgt ~ hgt, data = bdims)

# Show the coefficients
coef(mod)

## (Intercept)         hgt 
## -105.011254    1.017617

# Show the full output
summary(mod)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.743  -6.402  -1.231   5.059  41.103 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -105.01125    7.53941  -13.93   <2e-16 ***
## hgt            1.01762    0.04399   23.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared:  0.5145, Adjusted R-squared:  0.5136 
## F-statistic: 535.2 on 1 and 505 DF,  p-value: < 2.2e-16

Español

Valores ajustados y residuales

Una vez que ha ajustado un modelo de regresión, a menudo le interesan los valores ajustados (y ^ i ) y los residuos (ei), donde i

indexa las observaciones. Recordar que:

ei = yi − y ^ i

El procedimiento de ajuste de mínimos cuadrados garantiza que la media de los residuos sea cero (n.b., la inestabilidad numérica puede provocar que los valores calculados no sean exactamente cero). Al mismo tiempo, la media de los valores ajustados debe ser igual a la media de la variable de respuesta.

En este ejercicio, confirmaremos estos dos hechos matemáticos accediendo a los valores ajustados y los residuales con las funciones ajustado.values () y residuals (), respectivamente, para el siguiente modelo:

mod <- lm (wgt ~ hgt, data = bdims)

Instrucciones 100 XP

mod (definido anteriormente) está disponible en su espacio de trabajo.

Confirme que la media de los pesos corporales es igual a la media de los valores ajustados de mod.
Calcule la media de los residuos de mod.

# Mean of weights equal to mean of fitted values?
mean(bdims$wgt) == mean(fitted.values(mod))

## [1] TRUE

# Mean of the residuals
mean(residuals(mod))

## [1] -1.538036e-15

Poner en orden su modelo lineal

A medida que se ajusta a un modelo de regresión, hay algunas cantidades (por ejemplo, R2) que se aplican al modelo como un todo, mientras que otros se aplican a cada observación (por ejemplo, y ^ i) Si hay varias de estas cantidades por observación, a veces es conveniente adjuntarlas a los datos originales como nuevas variables.

La función augment () del paquete de la escoba hace exactamente esto. Toma un objeto modelo como argumento y devuelve un marco de datos que contiene los datos en los que se ajustó el modelo, junto con varias cantidades específicas del modelo de regresión, incluidos los valores ajustados, los residuos, las puntuaciones de apalancamiento y los residuos estandarizados.

El mismo modelo lineal del último ejercicio, mod, está disponible en su espacio de trabajo.

Cargue el paquete de la escoba.
Cree un nuevo marco de datos llamado bdims_tidy que es el aumento del modelo lineal mod.
Vea el marco de datos bdims_tidy usando glimpse ().

# Load broom
library(broom)

# Create bdims_tidy
bdims_tidy <- augment(mod)

# Glimpse the resulting data frame
glimpse(bdims_tidy)

## Rows: 507
## Columns: 9
## $ wgt        <dbl> 65.6, 71.8, 80.7, 72.6, 78.8, 74.8, 86.4, 78.4, 62.0, 81.6…
## $ hgt        <dbl> 174.0, 175.3, 193.5, 186.5, 187.2, 181.5, 184.0, 184.5, 17…
## $ .fitted    <dbl> 72.05406, 73.37697, 91.89759, 84.77427, 85.48661, 79.68619…
## $ .se.fit    <dbl> 0.4320546, 0.4520060, 1.0667332, 0.7919264, 0.8183471, 0.6…
## $ .resid     <dbl> -6.4540648, -1.5769666, -11.1975919, -12.1742745, -6.68660…
## $ .hat       <dbl> 0.002154570, 0.002358152, 0.013133942, 0.007238576, 0.0077…
## $ .sigma     <dbl> 9.312824, 9.317005, 9.303732, 9.301360, 9.312471, 9.314716…
## $ .cooksd    <dbl> 5.201807e-04, 3.400330e-05, 9.758463e-03, 6.282074e-03, 2.…
## $ .std.resid <dbl> -0.69413418, -0.16961994, -1.21098084, -1.31269063, -0.721…

Haciendo predicciones

La función fits.values () o el marco de datos augment () ed nos proporcionan los valores ajustados para las observaciones que estaban en los datos originales. Sin embargo, una vez que hayamos ajustado el modelo, es posible que queramos calcular los valores esperados para las observaciones que no estaban presentes en los datos en los que se ajustó el modelo. Este tipo de predicciones se llaman fuera de muestra.

El marco de datos de ben contiene una observación de altura y peso para una persona. El objeto mod contiene el modelo ajustado para el peso en función de la altura para las observaciones en el conjunto de datos bdims. Podemos usar la función predic () para generar valores esperados para el peso de nuevos individuos. Debemos pasar el marco de datos de nuevas observaciones a través del argumento newdata.

El mismo modelo lineal, mod, se define en su espacio de trabajo.

Imprimir ben a la consola.
Use predic () con el argumento newdata para calcular el peso esperado del individuo en el marco de datos de ben.

ben <- data.frame(wgt = 74.8, hgt = 182.8)

# Print ben
ben

# Predict the weight of ben
predict(mod, ben)

##        1 
## 81.00909

Agregar una línea de regresión a un gráfico manualmente

La función geom_smooth () facilita la adición de una línea de regresión lineal simple a un diagrama de dispersión de las variables correspondientes. Y, de hecho, hay modelos de regresión más complicados que se pueden visualizar en el espacio de datos con geom_smooth (). Sin embargo, todavía puede haber momentos en los que deseemos agregar líneas de regresión a nuestro diagrama de dispersión manualmente. Para hacer esto, usaremos la función geom_abline (), que toma argumentos de pendiente e intercepción. Naturalmente, tenemos que calcular esos valores con anticipación, pero ya vimos cómo hacerlo (por ejemplo, usando coef ()).

El marco de datos de coefs contiene las estimaciones del modelo recuperadas de coef (). Pasar esto a geom_abline () como argumento de datos le permitirá dibujar una línea recta en su diagrama de dispersión.

Use geom_abline () para agregar una línea definida en el marco de datos de coeficientes a un diagrama de dispersión de peso versus altura para individuos en el conjunto de datos bdims.

tittles <- c(names(coef(mod)))
coefs <- data.frame(coef(mod)[1],coef(mod)[2])
colnames(coefs) <- c(tittles)
rownames(coefs) <- c(1)

# Add the line to the scatterplot
ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
  geom_point() + 
  geom_abline(data = coefs, 
              aes(intercept = `(Intercept)`, slope = `hgt`),  
              color = "dodgerblue")

Error estándar de residuos

Una forma de evaluar la fuerza de ajuste es considerar qué tan lejos está el modelo para un caso típico. Es decir, para algunas observaciones, el valor ajustado estará muy cerca del valor real, mientras que para otros no. La magnitud de un residuo típico puede darnos una idea general de cuán cercanas están nuestras estimaciones.

Sin embargo, recuerde que algunos de los residuos son positivos, mientras que otros son negativos. De hecho, el procedimiento de ajuste de mínimos cuadrados garantiza que la media de los residuos es cero. Por lo tanto, tiene más sentido calcular la raíz cuadrada del residual cuadrático medio o el error cuadrático medio (RMSE)

) R llama a esta cantidad el error estándar residual.

Para que esta estimación sea imparcial, debe dividir la suma de los residuos al cuadrado por los grados de libertad en el modelo. Así,

RMSE = ∑ie2id.f .−−−−− √ = SSEd.f .−−−−− √

Puede recuperar los residuos de mod con residuales () y los grados de libertad con df.residual ().

Ver un resumen () de mod. Calcule la media de los residuos () y verifique que sea aproximadamente cero. Utilice los residuales () y df.residual () para calcular el error cuadrático medio (RMSE), también conocido como error estándar residual.

# View summary of model
summary(mod)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.743  -6.402  -1.231   5.059  41.103 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -105.01125    7.53941  -13.93   <2e-16 ***
## hgt            1.01762    0.04399   23.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared:  0.5145, Adjusted R-squared:  0.5136 
## F-statistic: 535.2 on 1 and 505 DF,  p-value: < 2.2e-16

# Compute the mean of the residuals
mean(residuals(mod))

## [1] -1.538036e-15

# Compute RMSE
sqrt(sum(residuals(mod)^2) / df.residual(mod))

## [1] 9.30804

Evaluar el ajuste del modelo lineal simple

Recordemos que el coeficiente de determinación (R2 ), se puede calcular como R2 = 1 − SSE / SST = 1 − Var (e) / Var (y), donde e es el vector de los residuos e y es la variable de respuesta. Esto nos da la interpretación de R2

como el porcentaje de la variabilidad en la respuesta que explica el modelo, ya que los residuales son la parte de esa variabilidad que permanece sin explicar por el modelo.

El marco de datos bdims_tidy es el resultado de aumentar () - el marco de datos bdims con el mod para wgt en función de hgt.

Use la función summary () para ver los resultados completos de mod.
Use el marco de datos bdims_tidy para calcular el R2

de mod manualmente usando la fórmula anterior, calculando la relación de la varianza de los residuos a la varianza de la variable de respuesta.

# View model summary
summary(mod)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.743  -6.402  -1.231   5.059  41.103 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -105.01125    7.53941  -13.93   <2e-16 ***
## hgt            1.01762    0.04399   23.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared:  0.5145, Adjusted R-squared:  0.5136 
## F-statistic: 535.2 on 1 and 505 DF,  p-value: < 2.2e-16

# Compute R-squared
bdims_tidy %>%
  summarize(var_y = var(wgt), var_e = var(.resid)) %>%
  mutate(R_squared = 1 - var_e / var_y)

Interpretación de R ^ 2

El R2

reportado para el modelo de regresión para la tasa de pobreza de los condados de EE. UU. en términos de tasa de graduación de secundaria es 0.464.

lm (fórmula = pobreza ~ hs_grad, data = countyComplete)%>% resumen()

¿Cómo debe interpretarse este resultado?

El 46.4% de la variabilidad en la tasa de graduados de secundaria entre los condados de EE. UU. Puede explicarse por la tasa de pobreza.
El 46.4% de la variabilidad en la tasa de pobreza entre los condados de EE. UU. Puede explicarse por la tasa de graduación de la escuela secundaria. ¡Derecha!
Este modelo es 46.4% efectivo.
La correlación entre la tasa de pobreza y la tasa de graduación de la escuela secundaria es de 0.464.

al versus promedio

El R2 nos da una medición numérica de la fuerza de ajuste en relación con un modelo nulo basado en el promedio de la variable de respuesta: y ^ nulo = y¯

Este modelo tiene un R2 de cero porque SSE = SST. Es decir, dado que los valores ajustados (y ^ nulo) son todos iguales al promedio (y¯

), el residuo para cada observación es la distancia entre esa observación y la media de la respuesta. Como siempre podemos ajustar el modelo nulo, sirve como una línea de base contra la cual se compararán todos los demás modelos.

En el gráfico, visualizamos los residuos para el modelo nulo (mod_null a la izquierda) frente al modelo de regresión lineal simple (mod_hgt a la derecha) con la altura como una variable explicativa única. Intente convencerse de que, si cuadra las longitudes de las flechas grises a la izquierda y las suma, obtendrá un valor mayor que si realizara la misma operación en las flechas grises a la derecha.

Puede ser útil obtener una vista previa de estos marcos de datos aumentados () ed con glimpse ():

vislumbre (mod_null) vislumbre (mod_hgt)

Calcule la suma de los residuos cuadrados (SSE) para el modelo nulo mod_null.
Calcule la suma de los residuos cuadrados (SSE) para el modelo de regresión mod_hgt.

mod_hgt <- augment(lm(wgt ~ hgt, data = bdims))
mod_null <- augment(lm(wgt ~ 1, data = bdims))

# Compute SSE for null model
mod_null %>%
  summarize(SSE = sum(.resid^2))

# Compute SSE for regression model
mod_hgt %>%
  summarize(SSE = sum(.resid^2))

Apalancamiento

El apalancamiento de una observación en un modelo de regresión se define completamente en términos de la distancia de esa observación desde la media de la variable explicativa. Es decir, las observaciones cercanas a la media de la variable explicativa tienen un bajo apalancamiento, mientras que las observaciones lejos de la media de la variable explicativa tienen un alto apalancamiento. Los puntos de alto apalancamiento pueden o no ser influyentes.

La función augment () del paquete de escoba agregará los puntajes de apalancamiento (.hat) a un marco de datos modelo. Instrucciones 100 XP

Use augment () para enumerar las 6 observaciones principales por sus puntajes de apalancamiento, en orden descendente.

# Rank points of high leverage
mod %>%
  augment() %>%
  arrange(desc(.hat)) %>%
  head()

Influencia

Como se señaló anteriormente, las observaciones de alto apalancamiento pueden o no ser influyentes. La influencia de una observación depende no solo de su influencia, sino también de la magnitud de su residuo. Recuerde que si bien el apalancamiento solo tiene en cuenta la variable explicativa (x ), el residual depende de la variable de respuesta (y) y del valor ajustado (y ^

)

Es probable que los puntos de influencia tengan un alto apalancamiento y se desvíen de la relación general entre las dos variables. Medimos la influencia usando la distancia de Cook, que incorpora tanto el apalancamiento como el residuo de cada observación. Instrucciones 100 XP

Use augment () para enumerar las 6 observaciones principales por la distancia de su Cook (.cooksd), en orden descendente.

# Rank influential points
mod %>%
  augment() %>%
  arrange(desc(.cooksd)) %>%
  head()

Eliminar valores atípicos

Las observaciones pueden ser atípicas por varias razones diferentes. Los estadísticos siempre deben ser cuidadosos, y lo que es más importante, transparentes, cuando se trata de valores atípicos. A veces, se puede lograr un mejor ajuste del modelo simplemente eliminando los valores atípicos y volviendo a ajustar el modelo. Sin embargo, uno debe tener una fuerte justificación para hacer esto. Un deseo de tener un R2 más alto

¡No es una razón suficientemente buena!

En los datos de mlbBat10, el valor atípico con un OBP de 0.550 es Bobby Scales, un jugador de cuadro que tuvo cuatro hits en 13 turnos al bate para los Cachorros de Chicago. Scales también caminó siete veces, lo que resultó en su OBP inusualmente alto. La justificación para eliminar Escalas aquí es débil. Si bien su desempeño fue inusual, no hay nada que sugiera que no sea un punto de datos válido, ni hay una buena razón para pensar que de alguna manera aprenderemos más sobre los jugadores de las Grandes Ligas al excluirlo.

Sin embargo, podemos demostrar cómo eliminarlo afectará nuestro modelo. Instrucciones 100 XP

Use filter () para crear un subconjunto de mlbBat10 llamado nontrivial_players que consta solo de aquellos jugadores con al menos 10 turnos al bate y OBP de menos de 0.500.
Ajuste el modelo lineal para SLG en función de OBP para los jugadores no triviales. Guarde el resultado como mod_cleaner.
Vea el resumen () del nuevo modelo y compare la pendiente y R2 a los de mod, el modelo original se ajusta a los datos de todos los jugadores.
Visualice el nuevo modelo con ggplot () y las funciones geom _ * () apropiadas.

# Create nontrivial_players
nontrivial_players <- mlbBat10 %>%
  filter(10 <= AB, OBP < 0.500)

# Fit model to new data
mod_cleaner <- lm(SLG ~ OBP, data = nontrivial_players)

# View model summary
summary(mod_cleaner)

## 
## Call:
## lm(formula = SLG ~ OBP, data = nontrivial_players)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31383 -0.04165 -0.00261  0.03992  0.35819 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.043326   0.009823  -4.411 1.18e-05 ***
## OBP          1.345816   0.033012  40.768  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07011 on 734 degrees of freedom
## Multiple R-squared:  0.6937, Adjusted R-squared:  0.6932 
## F-statistic:  1662 on 1 and 734 DF,  p-value: < 2.2e-16

# Visualize new model
ggplot(data = nontrivial_players, aes(x = OBP, y = SLG)) +
  geom_point() + 
  geom_smooth(method = "lm")

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

Elevados puntos de apalancamiento

No todos los puntos de alto apalancamiento son influyentes. Si bien la observación de alto apalancamiento correspondiente a las escalas de Bobby en el ejercicio anterior es influyente, las tres observaciones para jugadores con valores OBP y SLG de 0 no son influyentes.

Esto se debe a que de todos modos se encuentran cerca de la regresión. Por lo tanto, mientras que su OBP extremadamente bajo les da el poder de ejercer influencia sobre la pendiente de la línea de regresión, su bajo SLG les impide usarlo.

El modelo lineal, mod, está disponible en su espacio de trabajo. Use una combinación de aumentar (), organizar () con dos argumentos y encabezar () para encontrar las 6 observaciones principales con el mayor apalancamiento pero la menor distancia de Cook.

mod <- lm(formula = SLG ~ OBP, data = filter(mlbBat10, AB >= 10))

# Rank high leverage points
mod %>%
  augment() %>%
  arrange(desc(.hat)) %>%
  head()