library(ggplot2)
library(kableExtra)

Regresión Lineal Simple

La relación entre dos variables puede ser una de dependencia funcional; en este caso la magnitud de una de las variables (la variable dependiente o respuesta) se asume que está determinada - es función de - por la magnitud de la segunda variable (la variable independiente o predictora); mientras que lo contrario no es verdad.

Por ejemplo, la presión sanguínea en los humanos está relacionada con la edad, y podemos considerar a la presión como la variable dependiente, y a la edad como la independiente; lo contrario no tiene mucho sentido, pero tampoco podemos decir que la edad es el único factor que determina la presión sanguínea, o que exista una directa relación causa-efecto entre la presión y la edad.

Esta relación de dependencia se denomina regresión, y es simple porque se consideran solo dos variables. Como veremos más adelante, el término lineal se refiere al carácter aditivo del modelo a usar.

Es muchos casos, sin embargo, la relación entre dos variables no es de dependencia, aunque una cambie en magnitud y de manera equivalente la otra; en este caso tenemos una correlación entre las dos variables.

Ejemplo 1
Vamos a usar los datos de la longitud (cm) del ala derecha de gorriones de diferenctes edades (días). En este caso, podemos asumir que la variable dependiente es la longitud del ala y la independiente, la edad de las aves.

edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
sparrow <- data.frame(edad,l.ala)
#tabla de datos
knitr::kable(sparrow, col.names = c("Edad (días)", "Longitud ala (cm)"), caption = 'Ejemplo 1.  Longitud del ala derecha de gorriones de diferente edad')
Ejemplo 1. Longitud del ala derecha de gorriones de diferente edad
Edad (días) Longitud ala (cm)
3 1.4
4 1.5
5 2.2
6 2.4
8 3.1
9 3.2
10 3.2
11 3.9
12 4.1
14 4.7
15 4.5
16 5.2
17 5.0

Gráfica de los datos

#gráfica x-y
plot(edad,l.ala,xlab = "Edad, días", ylab = "Longitud, cm")

Figura 1. Gráfica de puntos de la edad (días), y la correspondiente longitud del ala (cm), en gorriones.

Modelo de regresión lineal simple

El modelo básico de regresión lineal simple, tomando en cuenta la variación no explicada por el modelo (residual o error), es el siguiente: \[Y_i = \alpha + \beta X_i + \epsilon_i\] Podemos observar que este modelo corresponde a la ecuación de una línea recta. La línea recta con un intercepto \(\alpha\), pendiente \(\beta\), y una suma de \(\epsilon_i\) igual a cero, se considera como el ‘mejor ajuste’ a los puntos con coordenadas \((X_i,Y_i)\).

Si todos los puntos cayeran en la recta (es decir: \(\epsilon_i = 0\), para todo i), el cálculo de \(\alpha\) y \(\beta\) sería el de una línea recta. Dado que este no es el caso, se utiliza el método de los cuadrados mínimos (‘least squares’), el cual consiste en encontrar el mínimo valor de esta expresión:

\[\sum_{i=1}^n (Y_i - \hat{Y_i})^2\] es decir el mínimo valor de la suma de las diferencias verticales de la Y de cada punto y la Y sobre la línea de regresión, elevada al cuadrado (ver la Figura 2).

Figura 2. Desviaciones de \(Y_i\) con respecto a \(\bar Y\) y \(\hat Y\). (tomada de Zar, 2014)

Cálculos del Coeficiente de Regresión y el Intercepto en Y

Para el cálculo de los estimadores de \(\alpha\) y \(\beta\), según el método de los cuadrados mínimos, se necesitan las siguientes cantidades:

Suma de Cuadrados \(\sum x^2\)
(todas las sumatorias son de i = 1 hasta n, número de puntos)
\[\sum(X_i - \bar X)^2 = \sum X_i^{2} - \frac{(\sum X_i)^2}{n} \]

Suma de Productos Cruzados \(\sum xy\) \[\sum(X_i - \bar X)(Y_i - \bar Y) = \sum X_iY_i - \frac{(\sum X_i)(\sum Y_i)}{n}\]

El coeficiente de regresión (parámetro \(\beta\)) es la pendiente de la línea de regresión y su estimador (b) es el siguiente:
\[b = \frac{\sum xy}{\sum x^2}\]

Cálculos manuales (a y b)

Ejercicio en clase: Calcular b con los datos del Ejemplo 1 (edad y longitud del ala).

#tamaño muestra
n <- length(l.ala)
#medias
X <- mean(sparrow$edad)
Y <- mean(sparrow$l.ala)
#sumatorias
sumX <- sum(sparrow$edad)
sumY <- sum(sparrow$l.ala)
sumX2 <- sum(sparrow$edad^2)
sumXY <- sum(sparrow$edad*sparrow$l.ala)
#suma de cuadrados
SSx <- sumX2 - (sumX^2)/n
Sxy <- sumXY - (sumX*sumY)/n
#cálculo de beta
est.b <- Sxy/SSx
#tablas
df1 <- data.frame(n,X,Y)
kable(df1, format = "markdown", col.names = c("n", "media Edad", "media Longitud"))
n media Edad media Longitud
13 10 3.415385 
df2 <- data.frame(sumX,sumY,sumX2,sumXY)
kable(df2, format = "markdown", col.names = c("sumatoria-X", "sumatoria-Y", "sumatoria-X^2", "sumatoria-XY"))
sumatoria-X sumatoria-Y sumatoria-X^2 sumatoria-XY
130 44.4 1562 514.8 
df3 <- data.frame(SSx,Sxy,est.b)
kable(df3, format = "markdown", col.names = c("SumaCuadradosX","SumaProductoCruzado","Pendiente"))
SumaCuadradosX SumaProductoCruzado Pendiente
262 70.8 0.270229

Para calcular el estimador de \(\alpha\), a, utilizamos el estimador b y las medias de X y Y: \[a = \bar Y - b\bar X\]

#estimador del intercepto, a
est.a <- Y - est.b*X
#ecuación de la recta de regresión
print('Ecuación de la recta de regresión:')
## [1] "Ecuación de la recta de regresión:"
sprintf("Y = %.3g + %.3fX", est.a, est.b)
## [1] "Y = 0.713 + 0.270X"
Cálculo de los coeficientes usando R

La función lm (‘linear model’) produce una lista (en este caso llamada reg1). Esta lista contiene los resultados del análisis de regresión, para el modelo: l.ala = a + b*edad + residual(error). Por ahora estamos interesados en el valor de los coeficientes (función coefficients) de la recta de regresión: intercepto (a) y el coeficiente de regresión (b) para la variable independiente (edad).

edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
reg1 <- lm(l.ala ~ edad)
coefficients(reg1)
## (Intercept)        edad 
##   0.7130945   0.2702290

Gráfica de la recta de regresión

Usando las funciones plot y abline podemos obtener una gráfica de los puntos y la recta de regresión, usando los datos de la lista reg1.

plot(edad,l.ala, xlab = 'Edad, días', ylab = 'Longitud del ala, cm', asp = 1)
abline(reg1, col = 'blue')

Figura 3. Línea de regresión para la relación entre la longitud del ala (cm) de gorriones y la edad (días).

¿Qué asumimos para una regresión paramétrica?

El análisis de regresión paramétrico incluye varias suposiciones, pero en general no interfieren con el cálculo de los coeficientes de la regresión. Sin embargo, para obtener un modelo confiable podemos examinar cómo se distribuyen los residuales (errores) del estimado de Y (\(\hat Y\)):

plot(edad, residuals(reg1), xlab = "Edad, días", ylab = "Residuales")

Figura 4. Gráfica de los residuales de la regresión lineal de la Figura 3.

Si los residuales no muestran una tendencia definida, podemos pensar que tienen una distribución aleatoria, lo cual es una de las suposiciones principales del análisis de regresión.

Prueba de hipótesis para la regresión.

Luego de obtener los coeficientes de la regresión, podemos realizar una prueba de hipótesis sobre el coeficiente b (estimado de \(\beta\)): \[H_0:\beta = 0\qquad H_A:\beta \neq 0\] Necesitamos tres nuevas cantidades para realizar una prueba usando el estadístico F:

Suma de Cuadrados de Y o Total \(\sum y^2\):
\[total\ SS = \sum Y_i{^2} - \frac{(\sum Y_i)^2}{n}\] que tiene n - 1 grados de libertad, total GL = 12.

Suma de Cuadrados de la Regresión: \[regresión\ SS = \frac{(\sum xy)^2}{\sum x^2} = b\sum xy\]

Suma de Cuadrados de Residuales: \[residual\ SS = total\ SS - regresión\ SS\]

Tabla 1: Resumen de cálculos para la prueba de hipótesis del coeficiente de regresión b.
Fuente.de.Variacion SS GL MS
Total total SS n - 1
Regresión regresión SS 1 reg.SS/reg.GL
Residual residual SS n - 2 res.SS/res.GL

Ahora podemos calcular el estadístico F: \[F = \frac{regresión\ MS}{residual\ MS}\] con un valor crítico: \[F_{\alpha(1), \nu_1, \nu_2}\] dónde \(\nu_1\) son los grados de libertad de la regresión (1 para una regresión lineal simple), y \(\nu_2\) son los grados de libertad de los residuales (n - 2).

Cálculos manuales (suma de cuadrados totales, de la regresión y de los residuales, F calculado)
edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
n <- length(edad)
reg1 <- lm(l.ala ~ edad)
coe <- as.matrix(coefficients(reg1))
b <- coe[2,1]
#Total SS
TotalSS <- sum(l.ala^2) - (sum(l.ala)^2)/n
#Regresión SS
RegresSS <- b*(sum(l.ala*edad)-(sum(edad)*sum(l.ala))/n)
#Residual SS
ResidualSS <- TotalSS - RegresSS
#MS
RegresMS <- RegresSS/1
ResidualMS <- ResidualSS/(n-2)
#tabla fuentes de variaciones
SS <- c(TotalSS, RegresSS, ResidualSS)
GL <- c(n-1, 1, n-2)
MS <- c(" ", RegresMS, ResidualMS)
tabla2 <- data.frame(
  Fuente.de.Variacion = c("Total", "Regresión", "Residual"),
  SS, GL, MS)
kable(tabla2) %>%
  kable_styling(full_width = F) %>%
  column_spec(1, bold = T, border_right = T) %>%
  column_spec(2:4, width = "10em", background = "orange")
Fuente.de.Variacion SS GL MS
Total 19.6569231 12
Regresión 19.1322137 1 19.132213740458
Residual 0.5247093 11 0.0477008487695559 
#Fs
Fcalc <- RegresMS/ResidualMS
Ftabla <- qf(.95,1,n-2)
tabla3 <- data.frame(Fcalc, Ftabla, row.names = " ")
kable(tabla3, format = "markdown")
Fcalc Ftabla
401.0875 4.844336

Conclusión sobre la prueba: El valor del F calculado es mucho mayor que el F de la tabla (\(F_{0.05(1),1,11}\)), por lo tanto podemos rechazar con seguridad la hipótesis nula, y proponer que la pendiente de la línea de regresión es diferente de cero.

Coeficiente de Determinación, \(r^2\)

La proporción (o porcentaje, si se multiplica por 100) de la variación total en Y que es explicada por el modelo de regresión (lineal, en este caso), se denomina coeficiente de determinación y se representa como \(r^2\). Se calcula de la siguiente manera: \[r^2 = \frac{regresión\ SS}{total\ SS}\] A \(r^2\) también se considera que expresa la bondad de ajuste de la línea a los datos, o la precisión de la regresión.

Análisis de regresión simple con R

La función lm provee todos los resultados para un análisis de regresión, y estos pueden obtenerse usando la función summary(objeto con lm).

edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
reg1 <- lm(l.ala ~ edad)
summary(reg1)
## 
## Call:
## lm(formula = l.ala ~ edad)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30699 -0.21538  0.06553  0.16324  0.22507 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.71309    0.14790   4.821 0.000535 ***
## edad         0.27023    0.01349  20.027 5.27e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2184 on 11 degrees of freedom
## Multiple R-squared:  0.9733, Adjusted R-squared:  0.9709 
## F-statistic: 401.1 on 1 and 11 DF,  p-value: 5.267e-10
  • Call: variables sobre las que se aplica la función lm, variable dependiente ~ variable predictora.
  • Residuals: cuartiles de la distribución de los valores del error (residuales); ayuda a determinar si los residuales tienen una distribución no-sesgada con respecto a la variable predictora.
  • Coefficients: estimados de los coeficientes de la regresión, sus errores estándares, el valor t para su diferencia con 0, y la probabilidad (error tipo I) para este valor t.
  • Residual standard error: el error estándar de los residuales, con n - 2 grados de libertad.
  • Multiple & Adjusted R-squared: los coeficientes de determinación (\(r^2\)), múltiple (total) y ajustado al número de variables.
  • F-statistic: el valor del estadístico F calculado, sus grados de libertad, y el valor p correspondiente.

Estimado de Y y sus limitaciones

Una vez que obtenemos los valores de a y b para la regresión lineal, podemos predecir el valor de la variable dependiente (\(\hat Y\)) dado un valor de \(X_i\) (Ejemplo 1): \[\hat Y = a + bX_i\] \[\hat Y = 0.715 + 0.270X_i\]

¿cuánto mide el ala de un gorrión de 13 días?

longitud del ala = 0.715 cm + (0.270 cm/día)(13 días) = 4.225 cm

Es importante tener en cuenta que la línea de regresión es aplicable para el ámbito de edades a partir del cuál se calculó, - en este caso, entre 3 y 17 días - fuera del mismo no podemos esperar que la regresión sea la misma. Por ejemplo, usando la ecuación de regresión obtenida: ¿cuál sería la longitud del ala de un gorrión de seis meses (180 días)?

Error Estándar e Intervalo de Confianza para el estimado de Y

El estimado de Y proviene de una ‘población’ de valores con un error estándar dado por la ecuación: \[s_{\hat Y_i} = \sqrt{s_{Y·X}^2[\frac{1}{n} + \frac{(X_i - \bar X)^2}{\sum x^2}]}\] donde \(s_{Y·X}^2\) es el cuadrado medio de residuales: \[ residual\ MS = \frac{residual\ SS}{n - 2}\]

Continuando con el Ejemplo 1, vamos a calcular el error estándar para el estimado de Y (4.225 cm) con X = 13 días.

#Y estimado = 4.225
#SumaCuadradosX = 262.0
#residual MS = 0.047701
#mediaX = 10
#n = 13
#error estándar del estimado de Y
syi <- sqrt(0.047701*(1/13 + ((13 - 10)^2)/262))
#intervalo de confianza del estimado
linf <- 4.225 - qt(.975, 11)*syi
lsup <- 4.225 + qt(.975, 11)*syi
result <- data.frame(syi, "4.225", linf, lsup)
kable(result, format = "markdown", col.names = c("Error Est.", "Y estimada", "Límite inferior", "Límite superior"))
Error Est. Y estimada Límite inferior Límite superior
0.0728553 4.225 4.064647 4.385353
Cálculos con R
edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
reg1 <- lm(l.ala ~ edad)
intconf <- predict.lm(reg1, newdata=list(edad=13.0), se.fit = TRUE, 
           interval="confidence", level=.95)
table4 <- data.frame(intconf)
kable(table4[,c(4,1:3)], format = "markdown")
se.fit fit.fit fit.lwr fit.upr
0.0728552 4.226072 4.065719 4.386425

El error estándar e intervalo de confianza del estimado están basados en la población de valores de nuestra muestra; pero también se pueden calcular como error estándar e intervalo de confianza de una predicción en particular.

intconf <- predict.lm(reg1, newdata=list(edad=13.0), se.fit = TRUE, 
           interval="prediction", level=.95)
table4 <- data.frame(intconf)
kable(table4[,c(4,1:3)], format = "markdown")
se.fit fit.fit fit.lwr fit.upr
0.0728552 4.226072 3.719325 4.732818
Gráfica del intervalo de confianza del estimado de Y

La línea de regresión puede ser acompañada con la franja del intervalo de confianza para los estimados de Y. Los puntos que caigan dentro de la franja, tienen una probabilidad dada (usualmente 95 %) de estar representados por la línea de regresión; esto es equivalente a decir que no son significativamente diferentes del valor estimado por la regresión, \(\hat Y\), para un \(\alpha\) dado.

gorrion <- data.frame(edad, l.ala)
ggplot(data=gorrion, aes(x=edad, y=l.ala)) +
  geom_point(pch=19, color="purple", size=2) +
  geom_smooth(method="lm", color="violet", linetype=1, level = 0.95) +
  labs(x="Edad, días", y="Long. ala, cm")

Figura 5. Línea de regresión con intervalo de confianza de 95 %, para la relación entre la longitud del ala de gorriones (cm) y la edad (días),

Preguntas
¿Cómo cambia la franja de intervalo de confianza si usamos una probabilidad de 99.9 %?

¿Dónde quedan los puntos de los datos?

¿Qué significa que los puntos estén dentro de la franja?

¿Qué tipo de error puede ocurrir al incluir más puntos en la franja?

Regresión Múltiple

Vamos a considerar el caso dónde tenemos una variable dependiente y más de una variable independiente. Usaremos datos de las etiquetas de diversos alimentos.

library(readxl)
datafood <- read_excel("datafood.xlsx", sheet = "fooddata")
head(datafood)

Estamos interesados en encontrar un modelo de regresión que nos indique cuáles variables determinan el contenido calórico de los alimentos.

modelo

supuestos

análisis de regresión múltiple

modRM <- lm(cal ~ serving + fat + chol + sodium + carbs + protein, 
              data = datafood)
summary(modRM)
## 
## Call:
## lm(formula = cal ~ serving + fat + chol + sodium + carbs + protein, 
##     data = datafood)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.1547  -0.7625   0.8633   1.9267   5.2134 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.022947   2.930160  -0.008    0.994    
## serving     -0.020711   0.028064  -0.738    0.471    
## fat          8.644047   0.280778  30.786 1.14e-15 ***
## chol         0.001097   0.121105   0.009    0.993    
## sodium       0.008753   0.006057   1.445    0.168    
## carbs        3.956071   0.117028  33.805 2.61e-16 ***
## protein      3.817924   0.409050   9.334 7.11e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.873 on 16 degrees of freedom
##   (7 observations deleted due to missingness)
## Multiple R-squared:  0.9948, Adjusted R-squared:  0.9929 
## F-statistic: 514.3 on 6 and 16 DF,  p-value: < 2.2e-16

Regresión Logística

Cuando la variable dependiente es una de tipo binomial: 0 o 1, positivo o negativo, falso o verdadero, podemos construir un modelo de regresión logística, para tratar de determinar cuáles son las variables independientes que determinan el resultado binario.

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
head(mydata)

modelo logit

mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
summary(mylogit)
## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
##     data = mydata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5802  -0.8848  -0.6382   1.1575   2.1732  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.449548   1.132846  -3.045  0.00233 ** 
## gre          0.002294   0.001092   2.101  0.03564 *  
## gpa          0.777014   0.327484   2.373  0.01766 *  
## rank        -0.560031   0.127137  -4.405 1.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 459.44  on 396  degrees of freedom
## AIC: 467.44
## 
## Number of Fisher Scoring iterations: 4

Referencias

Kabacoff, R. 2015. R in Action. Data analysis and graphics with R, Second Ed. ed. Manning, Shelter Island.

Rosner, B. 2016. Fundamentals of biostatistics, 8th edition. ed. Cengage Learning, Boston, MA.

Zar, J.H. 2014. Biostatistical analysis. Pearson India, Noida.