1 Regresión lineal simple y múltiple: teoría y práctica

Este tema es parte del Curso Fundamental de Posgrado: Análisis Estadístico en Ciencias Biológicas Utilizando R, semestre 2017-1, de la Universidad Nacional Autónoma de México, impartido entre Agosto-Diciembre 2016 en el Centro de Ciencias Genómicas. Para más información consultar la página del curso en: http://www.ccg.unam.mx/~vinuesa/R4biosciences/.

La parte teórica está basada en (Crawley 2015), (Davies 2016) y (A. P. Field, Miles, and Field 2012).

Este documento está aún en construcción y es generado con R (R Core Team 2016), rstudio (RStudio Team 2016), knitr (Xie 2016), rmarkdown (Allaire et al. 2016), pandoc (MacFerlane 2016) y \(LaTeX\).

1.1 Introducción: el concepto de regresión

El análisis de regresión engloba a un conjunto de métodos estadísticos que usamos cuando tanto la variable de respuesta como la la(s) variable(s) predictiva(s) son contínuas y queremos predecir valores de la primera en función de valores observados de las segundas. En esencia, el análisis de regresión consiste en ajustar un modelo a los datos, estimando coeficientes a partir de las observaciones, con el fin de predecir valores de la variable de respuesa a partir de una (regresión simple) o más variables (regresión múltiple) predictivas o explicativas.

El análisis de regresión juega un papel central en la estadística moderna y se usa para:

  • identificar a las variables predictivas relacionadas con una variable de respuesta
  • describir la forma de la relación entre estas variables y para derivar una función matemática óptima que modele esta relación
  • predecir la variable de respuesta a partir de la(s) explicativas o predictoras

1.1.1 Tipos de regresión

El término regresión puede ser confuso porque existen muchas variantes especializadas de regresión. Además, R tiene muchas funciones para ajustar una gran gama de modelos de regresión.

Tabla 1. Algunos tipos básicos de regresión (hay muchos más)

Tipo de regresión Uso típico
lineal simple Predicción de una variable de respuesta cuantitativa a partir de una variable predictora cuantitativa
polinomial Predicción de una variable de respuesta cuantitativa a partir de una variable predictora cuantitativa, donde la relación se modela como una función polinomial de orden \(n\)
lineal múltiple Predicción de una variable de respuesta cuantitativa a partir de dos o más variables predictoras cuantitativas
multivariada Predicción de más de una variable de respuesta cuantitativa a partir de una o más variables predictoras cuantitativas
logística Predicción de una variable categórica a partir de una o más predictoras
de Poisson Predicción de una variable de respuesta que representa un conteo a partir de una o más predictoras
no lineal Predicción de una variable de respuesta cuantitativa a partir de una o más predictoras, donde el modelo no es lineal
robusta Predicción de una variable de respuesta cuantitativa a partir de una o más predictoras, usando una aproximación resistente al efecto de observaciones influyentes

1.1.2 Modelado estadístico: selección de modelos, ajuste a los datos y varianza en la estima

Dadas dos variables cuantitativas existen virtualmente cientos de modelos que podrían describir la relación matemática entre ellas. El reto está en elegir el modelo que mejor se ajuste a estos datos para minimizar el error en la estima que se haga a partir del modelo.

Usamos los modelos para estimar el valor promedio de la variable de respuesta en función de parámetros estimados de los datos. De manera general, podemos predecir valores de la variable de respuesta usando esta fórmula:

\(estima_i = (modelo) + error_i\).

Para obtener la máxima precisión (mínimo \(error_i\)) en nuestra estima o predicción, tendremos que:

  1. elegir una familia de modelos adecuados a nuestros datos (modelos lineales, polinomiales, exponenciales, no lineales …)
  2. determinar el grado de parametrización adecuado del modelo
  3. obtener estimas de máxima verosimilitud de dichos parámetros

Sólo así podremos llegar a un compromiso óptimo entre realismo del modelo, grado de ajuste del modelo a los datos, y mínima varianza de la estima.
Veamos como ejemplo el ajuste de modelos lineal y cuadrático (polinomial) para predecir el peso de mujeres de 30-39 años en función de su altura (peso ~ altura) ¿Cuál de ellos se ajusta mejor a los datos?

  • Estadísticas asociadas al modelo lineal: \(\widehat{peso} = -87.52 +3.45 * altura\)
  • Residual standard error: 1.525 on 13 degrees of freedom
  • Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
  • F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14

  • Estadísticas asociadas al modelo polinomial: \(\widehat{peso} = 261.88 -7.35 * altura +0.083 * altura^2\)
  • Residual standard error: 0.3841 on 12 degrees of freedom
  • Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
  • F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16

Visualmente es claro que \(el modelo polinomial (cuadrático) se ajusta mejor a los datos\), como indican los valores de \(R^2\) (coef. de determinación) y los \(valores-p\) asociados al estadístico\(-F\).

¿Pero es este ajuste significativamente mejor? Comparando el ajuste de ambos modelos mediane el criterio de información de Akaike (AIC), confirmamos que el modelo polinomial es mucho mejor: aunque tiene un parámetro libre más a estimar de los datos (el coeficiente \(b2*altura^2\)) que el modelo lineal, el ajuste es ~ 40 unidades de AIC mejor, por lo que rechazamos el modelo lineal en favor el cuadrático.

##      df      AIC
## fit   3 59.08158
## fit2  4 18.51270

Aprenderemos más adelante los detalles relativos a la selección del mejor modelo de regresión, que es un aspecto esencial del modelado matemático.

En este capítulo nos enfocaremos en métodos que encajan en el rubro de modelos de regresión de mínimos cuadrados ordinarios (OLS). El énfasis estará en los modelos lineales simple y polinomiales, es decir, modelos en los que la analizamos el comportamiento de la variable de respuesta cuantitativa en función de una variable predictora o regresora cuantitativa, asumiendo una relación lineal o curvilínea entre ellas.

1.2 Regresión lineal simple

Se basa en modelos lineales con la fórmula genaral: \[\displaystyle Y_i = (a + bX_i) + \epsilon_i\]
donde:

  • a = punto de corte en el eje de ordenadas
  • b = pendiente o gradiente de la recta, que son los coeficientes de regresión
  • \(\epsilon_i\) corresponde al término de resíduos, que representa la diferencia entre el valor observado y el estimado para el individuo \(i\).

Los coeficientes de regresión los tenemos que estimar de los datos, usando el método de mínimos cuadrados, basado en las siguientes fórmulas, y el criterio de optimización de máxima verosimilitud. \[\displaystyle b = \frac{\sum\limits_{i=1}^n(x_i-\bar{x})(yi-\bar{y})} {\sum\limits_{i=1}^n(x_i-\bar{x})^2} = \frac{\sum{xy} - \frac{\sum{x}\sum{y}}{n}} {\sum{x^2} - \frac{\left(\sum{x}\right)^2}{n}} = \frac{SSXY}{SSX} = \frac{suma \, corregida \, de \, productos}{suma\,de \, cuadrados \, de \, x}\] \[\displaystyle a = \bar{y} - b\bar{x} = \frac{\sum{y}}{n} -b \frac{\sum{x}}{n}\]

Una vez estimados estos coeficientes los sustituímos en la fórmula general de la recta, definiendo así el modelo lineal que ajustamos en el ejemplo gráfico anterior: \(\widehat{peso} = -87.52 +3.45 * altura\). Nótese que en este caso el punto de corte no tiene mucho sentido ya que no encontraremos una persona de 0 pulgadas de altura.

1.2.1 Cálculo manual de los coeficientes de regresión de un modelo lineal

Es fácil estimar “a ojo” los coeficientes de regresión, razonablemente precisos, de un modelo lineal.
Ajustemos un modelo lineal a los siguientes datos:

En pseudocódigo los pasos a seguir son:

  1. Evaluar el decremento de la variable de respuesta: ha decrecido de 12 a 2, es decir un decremento de -10
  2. Determinar el rango de variación de la variable regresora: \(x = +8\)
  3. Calcular la pendiente o gradiente de la variación de la variable de respuesta: \(b\approx-10/8=-1.25\)
  4. Determinar el punto de corte: ¿Cuánto vale \(y\) cuando \(x = 0\)?, \(a\approx12\)

Por tanto la recta de regresión simple ajustada manualmente a los datos tiene la fórmula: \(y=12 - 1.25x\)

Estos coeficientes no difieren mucho de los obtenidos al ajustar un modelo lineal con R:

# Noten que la var. de respuesta va a la izda de la tilde ' ~' y se lee como 
# modelado de Y en función de X. Esta notación corresponde a lo que se conoce como 
# una fórmula en R. Las fórmulas son centrales en el modelado estadístico y tienen una 
# sintaxis particular y compleja, que hay que aprender para manejar R eficientemente
lm(dfr1$Y ~ dfr1$X)
## 
## Call:
## lm(formula = dfr1$Y ~ dfr1$X)
## 
## Coefficients:
## (Intercept)       dfr1$X  
##      11.756       -1.217

1.2.2 ¿Cómo calcula R los coeficientes de un modelo lineal de regresión?

R busca las estimas de máxima verosimilitud de los parámetros: \(L\approx Pr(D|M)\). Es decir, habiendo seleccionado un modelo lineal, necesitamos encontrar los valores de los parámetros del Modelo (\(a\) y \(b\)) que maximizan la probabilidad de observar los Datos.

1.2.2.1 El método de mínimos cuadrados

Este método garantiza encontrar el valor de máxima verosimilitud de los parámetros del modelo lineal para un conjunto de datos, al minimizar la magnitud de la diferencia entre los puntos y la recta. Veamos ésto gráficamente, al añadir la recta más ajustada a los datos medienta la función \(abline()\) y dibujando líneas que representan los resíduos, es decir, la distancia de cada punto a la recta.

# 1. Graficamos los puntos
plot(dfr1$X, dfr1$Y, pch=21, bg="blue", ylab= "Var. de respuesta", 
     xlab="Var. explicativa", main = "Recta de regresión y resíduos")

# 2. Añadimos a la gráfica de dispersión la recta de regresión 
abline(lm(dfr1$Y ~ dfr1$X), col="green") 

# Podemos obtener el valor predicho de Y para cada valor de X con predict()
fitted <- predict(lm(dfr1$Y ~ dfr1$X))

# Con un bucle for y la función lines() podemos representar los resíduos
for(i in 1:9)
{
  lines( c(dfr1$X[i],dfr1$X[i]), c(dfr1$Y[i], fitted[i]), col="red")
}

# Veamos los valores Y ajustados que calculamos previamente:
fitted
##         1         2         3         4         5         6         7 
## 11.755556 10.538889  9.322222  8.105556  6.888889  5.672222  4.455556 
##         8         9 
##  3.238889  2.022222

Los resíduos describen la bondad de ajuste de la recta de regresión. Nuestro modelo de máxima verosimilitud se define como el modelo que minimiza la suma de los cuadrados de estos resíduos.

1.2.2.2 Simulación del algoritmo que usa R para minimizar la suma de cuadrados

R implementa un algoritmo para minimizar esta suma de cuadrados. Concretamente el algoritmo, mediante cáluculo diferencial (tomando derivadas) calcula el mínimo de la función de verosimilitud para cada coeficiente a partir de la fórmula: \[\displaystyle d = y - \hat{y}\] Dado \(\hat{y}\) está sobre la recta \(a+bx\), podemos reescribir la ecuación anterior como: \[\displaystyle d = y - (a+bx)=y-a-bx\]

La estima por mínimos cuadrados de la pendiente de regresión (\(b\)) se calcula por tanto rotando la línea hasta que el error de la suma de cuadrados, \(SSE\), sea minimizada. El algoritmo por tanto tiene que encontrar el mínimo de \(\sum{(y-a-bx)^2}\). Tomando la derivada de SSE con respecto a \(b\): \[\displaystyle \frac{dSSE}{db} = -2 \sum{x(y-a-bx)}\]

Operando llegamos a que:

\[\displaystyle b = \frac{\sum{xy} - \frac{\sum{x}\sum{y}}{n}} {\sum{x^2} - \frac{\left(\sum{x}\right)^2}{n}}\]

Que se abrevia como:

\[\displaystyle b = \frac{SSXY}{SSX} = \frac{suma \, corregida \, de \, productos}{suma\,de \, cuadrados \, de \, x}\]

Para que tengan una impresión de los cálculos involucrados, es útil graficar la suma de cuadrados de los resíduos \(SSE\) contra el valor del parámetro que estamos tratando de estimar. Veamos el caso de la pendiente \(b\).

  • Sabemos que la recta de regresión pasa por un punto en el centro de la nube de puntos, cuyas coordenadas son \((\bar{x},\bar{y})\). La línea de mejor ajuste tendrá que pivotarse en torno a las medias de \(x\) e \(y\), seleccionando aquella línea cuya pendiente \(b\) minimiza la suma de cuadrados de las distancias marcadas en rojo en la gráfica anterior.
  • El valor de \(b\) que minimiza la \(SSE\) representa la estima de máxima verosimilitud de este parámetro. Veamos una simulación del algoritmo. Sabemos de nuestras estimas que \(b\approx -1.25\).
  • Por tanto vamos a localizar el valor de máxima verosimilitud en el rango de \(-1.4 < b < -1.0\), estimando la \(SSE\) para cada valor de \(b\) en este rango.

  • En pseudocódigo los pasos a realizar son:
  • recorrer en pequeños incrementos el valor de \(b\) dentro del rango establecido
  • calcular el punto de corte correspondiente (\(a=\bar{y}-b\bar{x}\) )
  • estimar los valores ajustados de \(y\) para cada nivel de x (\(y = a+bx\))
  • calcular los resíduos asociados (\(d = y-a-bx\)) a cada estima puntual
  • elevar los resíduos al cuadrado y sumarlos (\(\sum(y-a-bx)^2\))
  • asociar este valor de \(sse[i]\) con la estima actual de la pendient \(b[i]\)

Una vez ejecutado este algoritmo, produciremos una gráfica en forma de “U”, con la suma de resíduos cuadrados \(SSE\) en el eje de ordenadas, expresada en función de la estima de la pendiente \(b\). El mínimo de esta gráfica corresponderá a la estima de máxima verosimilitud de \(b\) para la cual se minimiza la \(SSE\).

# estima de maxima verosimilitud de la pendiente

b <-  seq(-1.43,-1,0.002)
sse <- numeric(length(b))
for (i in 1:length(b)) {
a <- mean(dfr1$Y)-b[i]*mean(dfr1$X)
residual <- dfr1$Y - a - b[i]*dfr1$X
sse[i] <- sum(residual^2)
}

# generemos la gŕafica de sse ~ b
plot(b,sse,type="l",ylim=c(19,24), 
     main = "Estima de máxima verosimilitud del parámetro b")

# añadamos una flecha en el mínimo
arrows(-1.216,20.07225,-1.216,19,col="red")

# dibujemos una recta verde que pase por el mínimo
abline(h=20.07225,col="green",lty=2)
lines(b,sse)

# el valor mínimo de sse lo obtenemos con la fución min()
min(sse)
## [1] 20.07225
# así obtenemos el valor de "b" asociado al valor mínimo de sse
b4minSSE <- b[which(sse==min(sse))]
b4minSSE
## [1] -1.216

1.2.2.3 Fórmulas se suma de cuadrados corregidas y cálculo de los coeficientes de regresión

En resumen, en base a las fórmulas antes enunciadas para los coeficientes de regresión:

\[\displaystyle b = \frac{\sum\limits_{i=1}^n(x_i-\bar{x})(yi-\bar{y})} {\sum\limits_{i=1}^n(x_i-\bar{x})^2} = \frac{\sum{xy} - \frac{\sum{x}\sum{y}}{n}} {\sum{x^2} - \frac{\left(\sum{x}\right)^2}{n}} = \frac{SSXY}{SSX} = \frac{suma \, corregida \, de \, productos}{suma\,de \, cuadrados \, de \, x}\] \[\displaystyle a = \bar{y} - b\bar{x} = \frac{\sum{y}}{n} -b \frac{\sum{x}}{n}\]

y recordando las fórmulas de las sumas de cuadrados y productos corregidas, que son absolutamente centrales en todo lo que seguirá sobre regresión y análisis de la varianza:

  • suma total de cuadrados \(SSY\): \[\displaystyle SSY = {\sum{y^2} - \frac{\left(\sum{y}\right)^2}{n}}\]
  • suma de cuadrados de x \(SSX\): \[\displaystyle SSX = {\sum{x^2} - \frac{\left(\sum{x}\right)^2}{n}}\]
  • suma de productos corregida: \[\displaystyle SSXY= \frac{\sum{xy} - \frac{\sum{x}\sum{y}}{n}} {\sum{x^2} - \frac{\left(\sum{x}\right)^2}{n}}\]

podemos calcular fácilmente los coeficientes de regresión, como muestra el siguiente bloque de código:

SSX <- sum(dfr1$X^2)-sum(dfr1$X)^2/length(dfr1$X)
SSY <- sum(dfr1$Y^2)-sum(dfr1$Y)^2/length(dfr1$Y)
SSXY <- sum(dfr1$X*dfr1$Y)-sum(dfr1$X)*sum(dfr1$Y)/length(dfr1$X)

cat("SSX =", SSX, "| SSY =", SSY, "| SSXY =", SSXY, "\n")
## SSX = 60 | SSY = 108.8889 | SSXY = -73
# por tanto:
b <- SSXY/SSX
a <- (sum(dfr1$Y)/length(dfr1$Y)) - (b*sum(dfr1$X)/length(dfr1$X))

cat("los coeficientes de regresión: a =", a, "| b =", b, "\n")
## los coeficientes de regresión: a = 11.75556 | b = -1.216667

1.2.3 Determinación de la bondad del ajuste: \(suma\ de\ cuadrados\), \(R\) y \(R^2\)

Lo que hemos visto hasta ahora es sólo la “mitad de la historia”. Además de los parámetros \(a\) y \(b\), necesitamos calcular la incertidumbre asociada a cada una de estas estimas. Es decir, necesitamos calcular el error estándar de cada parámetro de nuestro modelo de regresión.

El error estándar (\(SE\)) es la desviación estándar de la distribución de muestreo de un estadístico. Para un estadístico (p.ej. la \(media\)), nos indica cuánta variación existe entre muestras de la misma población. Por tanto valores grandes de \(SE\) indican que el estadístico no representa adecuadamente a la población de la que procede la muestra.

Para determinar la bondad del ajuste de un modelo necesitamos compararlo con un modelo nulo, el más simple de todos, la \(media\), que es un “modelo de no relación entre variables”. Cuantitativamente la determinamos usamos la siguiente ecuación general:

\[\displaystyle desviación = \sum{(observación - modelo)^2}\]

1.2.3.1 Explicación gráfica de los conceptos subyacentes a la bondad de ajuste: \(SSY\), \(SSE\) y \(SSR\)

La mejor manera de explicar este concepto es hacerlo gráficamente:

opar <- par(no.readonly = TRUE)

par(mfrow = c(1,2))

## I. Gráfica de desviaciones de los datos de la media (para cálculo de SSY)
# 1. Graficamos los puntos (gráfico de dispersión)
plot(dfr1$X, dfr1$Y, pch=21, bg="blue", ylab= "Y", xlab="X",
     main = "H0 = Media; +dev. (SSY)")

# 2. Añadimos una recta, a la altura de la media de Y, como modelo nulo
abline(h=mean(dfr1$Y), col="black") 

# 3. Con un bucle for y la función lines() podemos representar las desviaciones 
#    de la media
for(i in 1:9)
{
  lines( c(dfr1$X[i],dfr1$X[i]), c(mean(dfr1$Y), dfr1$Y[i]), col="blue", lty=3)
}

## II: Gráfica de resíduos (para cálculo de SSE)
# 2. Añadimos a la gráfica de dispersión la recta de regresión 
plot(dfr1$X, dfr1$Y, pch=21, bg="blue", ylab= "Y", xlab="X",
     main = "H1 = lm; + resíduos (SSE)")

abline(lm(dfr1$Y ~ dfr1$X), col="black") 

# Podemos obtener el valor predicho de Y para cada valor de X con predict()
fitted <- predict(lm(dfr1$Y ~ dfr1$X))


for(i in 1:9)
{
  lines( c(dfr1$X[i],dfr1$X[i]), c(dfr1$Y[i], fitted[i]), col="red", lty=3)
}

par(opar)
  • Definiciones
  • \(SSY\) o \(SST\) es la suma total de cuadrados, ya que representa la variación máxima o total, respecto del modelo nulo (media)
  • \(SSE\) es la suma cuadrada de resíduos y representa el grado de error del modelo ajustado, o variación no explicada por el modelo
  • \(SSM\) o \(SSR\) es la suma de cuadrados del modelo de regresión y mide la diferencia entre SSY y SSE, indicando la ganancia en precisión de la estima al ajustar un modelo, con respecto a la media, o variación explicada por el modelo

Por tanto \(SSY = SSR + SSE\) representa al particionado de la suma de cuadrados de una regresión.

La siguient gráfica muestra \(SSR = SSY - SSE\) visualmente, es decir, las diferencias entre el valor promedio de \(Y\) y las estimas del modelo de regresión. Si \(SSR\) es grande, quiere decir que las predicciones hechas por el modelo son muy distintas a usar la media para predecir valores la variable de respuesta. Cuanto mayor sea \(SSR\), mejor será la capacidad del modelo de predecir la variable de respuesta.

## I. Gráfica de desviaciones de los datos de la media (para cálculo de SSR)
# 1. Graficamos los puntos
plot(dfr1$X, dfr1$Y, pch=21, bg="blue", ylab= "Y", xlab="X",
     main = "SSR (SSM) = H1 - H0")

# 2. Añadimos una recta, a la altura de la media de Y, como modelo nulo
abline(h=mean(dfr1$Y), col="red") 

abline(lm(dfr1$Y ~ dfr1$X), col="green") 
fitted <- predict(lm(dfr1$Y ~ dfr1$X))
# 3. Con un bucle for y la función lines() podemos representar las desviaciones 
#    de la media
for(i in 1:9)
{
  lines( c(dfr1$X[i],dfr1$X[i]), c(mean(dfr1$Y), fitted[i]), col="red", lty=3)
}

1.2.3.2 Cálculo manual del particionado de las sumas de cuadrados en el análisis de regresión: análisis de la varianza

La idea que subyace a la fórmula \(SSY = SSR + SSE\) es sencilla: tomamos la variación total de los datos originales (\(SSY\)) y la particionamos entre dos componentes que nos revelan el poder explicativo o predictor del modelo. La variación explicada por el modelo \(SSR\) se conoce como la suma de cuadrados del modelo (de regresión), y la variación no explicada (\(SSE\)) se conoce como suma de cuadrados de los errores

Ya vimos cómo calcular \(SSY\) y \(SSE\). Por tanto \(SSR = SSY - SSE = \left({\sum{y^2} - \frac{\left(\sum{y}\right)^2}{n}}\right) - \left(\sum{(y-a-bx)^2}\right)\)

Pero podemos hacerlo más sencillo, si en vez de tomar \(SSE =\sum{(y-a-bx)^2}\), usamos el siguiente atajo:

\[\displaystyle SSR = b*SSXY\]

Podemos así calcular fácilmente \(SSR\) usando los valores de \(b\) y \(SSXY\) previamente calculados:

\[\displaystyle SSR = -1.21667 X -73 = 88.81667\]

Teniendo \(SSY\) y \(SSR\) calculados, es trivial calcular \(SSE\):

\[\displaystyle SSE = SSY - SSR = 108.8889 -88.81667 = 20.07222\]

Ya tenemos los dos componentes (\(SSR\) y \(SSE\)) en los que se descompone la variación total de los datos (\(SSY\)). Hagamos un análisis de la varianza, la famosa ANOVA, para comparar estos dos componentes mediante un cociente \(F\) entre las varianzas asociadas a cada componente.

Estadísticos como \(F\) representan generalmente la cantidad de varianza sistemática (del modelo) / varianza no sistemática (de los datos), es decir \(F\) se basa en la razón de \(SSR/SSE\). Pero como las sumas de cuadrados dependen del número de puntos (diferencias sumadas), \(F\) usa la media de la suma de cuadrados, \(MS = promedio\ de\ suma\ de\ cuadrados\).

\[\displaystyle F = \frac{MS_M}{MS_E} = \frac{media\ de\ los\ cuadrados\ del\ modelo}{media\ cuadrados\ de\ los\ errores\ o\ resíduos}\]

donde \[\displaystyle MS_M = \frac{SSR}{g.l.} = \frac{SSR}{no.\ parámetros\ estimados\ en\ el\ modelo:\ b\ (1\ g.l.\ en\ ml\ simple)}\]

y \[\displaystyle MS_E=\frac{SSE}{g.l.} = \frac{SSE}{no.\ observaciones - no.\ parámetros\ libres:\ a,\ b\ (n-2)}\]

Construyamos la tabla de ANOVA correspondiente, paso a paso:

Ya tenemos calculados las sumas de cuadrados de cada fuente de variación, como se muestra abajo.

Fuente Suma de cuadrados Grados de libertad Media de cuadrados (Varianza) cociente \(F\)
Regresión (modelo) 88.817 (\(SSR\))
Error (resíduos) 20.072 (\(SSE\))
Total (datos) 108.889 (\(SSY\))

Llenemos ahora la columna de los grados de libartad (\(g.l.\)). Recordemos que, esencialmente, los grados de libartad representan el número de ‘entidades’ que están libres de variar cuando estimamos algún parámetro estadístico a partir de los datos. Los \(g.l\) generalmente se calculan como \(n-número\ de\ parámetros\ libres\ a\ estimar\ de\ los\ datos\). Por tanto:

  • \(g.l.(SSY)\) Dado que \(SSY = \sum{(y-\hat{y})^2}\), vemos que sólo necesitamos estimar 1 parámetro: \(\hat{y}\); \(g.l. = (n-1)\)
  • \(g.l.(SSE)\) Dado que \(SSE = \sum{(y -a -bx)^2}\), vemos que necesitamos estimar 2 parámetros: \(a\) y \(b\); \(g.l. = (n-2)\)
  • \(g.l.(SSR)\) en el modelo de regresión lineal simple hemos estimado 1 parámetro \(b\); \(g.l. = 1\)

Noten que la suma de \(g.l.(SSR) + g.l.(SSE) = g.l.(SSY\)), como se ve en la tabla.

Fuente Suma de cuadrados Grados de libertad Media de cuadrados (Varianza) cociente \(F\)
Regresión (modelo) 88.817 (\(SSR\)) 1
Error (resíduos) 20.072 (\(SSE\)) 7 (\(n-2\))
Total 108.889 (\(SSY\)) 8 (\(n-1\))

Necesitamos llenar la 4a. columna, reservada para las varianzas. Recoredemos que:

\[\displaystyle var = s^2 = \frac{suma\ de\ cuadrados}{grados\ de\ libertad}\]

Esta columna (4) es muy fácil de llenar, ya que tenemos los valors que necesitamos pre-calculados en las columnas 2 y 3. La \(var\ total\) no nos hace falta, pero sí el \(cociente-F = \frac{88.817}{2.86764}\), que incluímos tambien en la tabla de ANOVA, completándola.

Fuente Suma de cuadrados Grados de libertad Media de cuadrados (Varianza) cociente \(F\)
Regresión (modelo) 88.817 (\(SSR\)) 1 88.817 (var. regresión) 30.974
Error (resíduos) 20.072 (\(SSE\)) 7 (\(n-2\)) \(s^2\) = 2.86764 (var. error)
Total 108.889 (\(SSY\)) 8 (\(n-1\))

Calculemos ahora la significancia del \(cociente-F\). Recordemos que:

  • \(H_0: b = 0\), la pendiente de la recta de regresión es cero (el modelo no difiere de la media)
  • \(H_1: b\ne0\), usando una alternativa de doble cola asumimos que la pendiente es positiva o negativa

Podemos buscar el valor crítico de \(F\) en la tabla correspondiente, usando \(g.l._numerador = 1\) y \(g.l._denominador = 7\) y veríamos algo como ésto:

\(g.l.(denom)\) p 1 … (\(g.l.(num.)\))
1 0.05 161.45
1 0.01 4052.18
7 0.05 5.59
7 0.01 12.25

Dado que nuestro estadístico \(F\) es mucho mayor que el valor crítico de \(F_{(1,7)}\), rechazamos la hipótesis nula con \(p < .001\).

Pero estamos aprendiendo estadística usando R, ¿verdad? ¿O cargas contigo las tablas de valores críticos para múltiples estadísticos?

Recordemos que los intervalos de confianza se definen como quantiles del 95% ó 99% … Podemos calcular los quantiles para \(F\) usando la función \(qf(.95,glNum,glDenom)\). La \(p\) asociada al estadístico \(F\) la calculamos como \(1-pf(F,glNum,glDenom)\), como se muestra seguidamente:

glNum <- 1
glDenom <- 7
Fstat <- 30.974

CI95_F <- qf(.95, glNum, glDenom)
CI99_F <- qf(.99, glNum, glDenom)

p <- 1-pf(Fstat,glNum,glDenom)

cat("CI95% = ", CI95_F, "| CI99% = ", CI99_F, " | p = ", p, "\n")
## CI95% =  5.591448 | CI99% =  12.24638  | p =  0.0008460725

Por tanto, la probabilidad de observar un \(cociente-F\) de 30.974 bajo la \(H_0\) es muy poco probable, de hecho \(p < 0.001\), por lo que rechazamos \(H_0\) contundentemente.

1.2.3.3 Cálculo en R del particionado de las sumas de cuadrados en el análisis de regresión: \(summary.aov()\)

Podemos fácilmente comprobar en R lo que hemos hecho a mano en la sección anterior, usando la función \(summary.aov()\) a la que le pasamos el objeto que regresa \(lm()\):

# ajustemos un modelo lineal
fit1 <- lm(dfr1$Y ~ dfr1$X)

# summary.aov() nos da la tabla de ANOVA
summary.aov(fit1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dfr1$X       1  88.82   88.82   30.97 0.000846 ***
## Residuals    7  20.07    2.87                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.2.3.4 Estima de la significancia de los coeficientes de regresión y cálculo de los errores asociados

Aunque el resultado del análisis de la varianza sea altamente significativo, no es particularmente interesante. En el contexto del análisis de regresión lo que realmente nos interesa conocer es el tamaño de los efectos (estimas de \(a\) y \(b\)), su significancia y errores asociados a su estima.

  • Los errores asociados a la estima de los coeficientes se calculan así:

\[\displaystyle SE_b = \sqrt{\frac{s^2}{SSX}} = \sqrt{\frac{2.867}{60}} = 0.2186\]

\[\displaystyle SE_a = \sqrt{\frac{s^2 \sum{x^2}}{n * SSX}} = \sqrt{\frac{2.867 * 204}{9 * 60}} = 1.0408\]

  • Cálculo de la significancia de \(b\)

Ya vimos que en el peor de los escenarios \(b = 0\), es decir, que un cambio unitario en la variable predictora no se asocia con cambio alguno en la variable de respuesta. Por tanto si una variable predice significativamente una respuesta, su coeficiente de regresión tiene que ser significativamente distinto a 0.

Esta hipótesis la podemos evaluar con una \(prueba-t\), donde \(H_0: b = 0\).

Al igual que el \(cociente-F\), el estadístico-\(t\) se basaen el cociente entre la varianza explicada y la no explicada o error. El test será significativo si \(b\) es grande comparado con la magnitud del error en la estima:

\[\displaystyle t = \frac{b_{observada} - b{esperada}}{SE_b} = \frac{b_{observada} - 0}{SE_b}\]

En regresión lineal simple, los grados de libertad para los predictores (\(p\)) se calculan como: \(N-p-1\), donde \(N\) tamaño de muestra, y \(p\) es el número de predictores (1 en este caso: \(b\)).

  • Toda esta información sobre el modelo lineal nos la da la función \(summary()\):
fit1 <- lm(dfr1$Y ~ dfr1$X)

# la salida de summary() nos da los tamaños de los efectos, es decir, los 
# coeficientes de regresión (a, b), los errores asociados, su significancia (t-test) + 
# el resultado de la prueba F (ANOVA)
summary(fit1)
## 
## Call:
## lm(formula = dfr1$Y ~ dfr1$X)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4556 -0.8889 -0.2389  0.9778  2.8944 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.7556     1.0408  11.295 9.54e-06 ***
## dfr1$X       -1.2167     0.2186  -5.565 0.000846 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.693 on 7 degrees of freedom
## Multiple R-squared:  0.8157, Adjusted R-squared:  0.7893 
## F-statistic: 30.97 on 1 and 7 DF,  p-value: 0.0008461

1.2.3.5 Cálculo del grado de ajuste, \(R^2\)

Nos queda un asunto muy importante que considerar. Dos rectas de regresión pueden tener la misma pendiente y puntos de corte, pero derivarse de relaciones completamente diferentes entre las variables:

## Coeficientes para fit1:  0.5 2 Coeficientes para fit2:  0.5 2

Claramente el ajuste del modelo a los datos en la gráfica derecha es mucho mejor que en la izquierda. Llevado al extremo ideal, los puntos caerían todos sobre la recta, reflejando una \(dispersión = 0\). En el otro extremo, \(x\) podría no explicar ninguna de la variación observada en \(y\). En dicho caso las variables no estarían correlacionadas, el ajuste del modelo lineal sería 0 y la magnitud de la dispersión sería de 100%.

Combinando lo que hemos aprendido sobre \(SSY\), \(SSR\) y \(SSE\) debería quedar claro que podemos derivar una métrica que mida la fracción de la variación total en \(y\) explicada por la regresión. A esta medida es a la que llamamos coeficiente de determinación: \(R^2 = \frac{SSR}{SSY}\).

\(R^2\) varía entre 0, cuando la regresión no explica ninguna de la variación obserbada en \(y\) (\(SSE = SSY\); dispersión = 100%), y 1, cuando la regresión explica toda la variación en y (\(SSR = SSY\); dispersión = 0%). Gráficamente:

## * Coeficientes para fit1:  0 0 
## * Coeficientes para fit2:  0.5 2

Esta \(R^2\) corresponde al coeficiente de determinación que vimos en el tema de correlación. Generalmente lo expresaremos como porcentaje (\(R^2 * 100\)), y representa la cantidad de varianza en la variable de respuesta explicada por el modelo (\(SSR\)) relativa al total de la variación existente en los datos (\(SSY\) o \(SST\)). Por tanto, \(R^2\) representa el porcentaje de la variación en la variable de respuesta que puede ser explicada por el modelo.

Además, tratándose de regresión simple, podemos calcular el coeficiente de correlación de Pearson entre las dos variables como \(r=\sqrt{R^2}\).

1.2.4 Diagnóstico del modelo de regresión y validación de supuestos

Un último aspecto central en el análisis de regresión lineal es comprobar que el modelo no viole los supuestos de los que depende. Particularmente críticos son:

  1. linearidad la variable de respuesta está relacionada linealmente con las variables predictoras
  2. constancia de las varianzas (homocedasticidad) : la varianza de la variable de respuesta no varía con los niveles de la predictora
  3. normalidad de los errores
  4. ausencia de puntos marcadamente aberrantes (outliers) palanca o influyentes

Una primera opción para comprobar si el modelo cumple estos supuestos es mediante la generación de cuatro gráficas de checado del modelo a partir de la información que viene integrada en el objeto que regresa la función \(lm()\) usando la función \(plot()\)

opar <- par(no.readonly = TRUE)

# ajustamos modelo lineal
modelo <- lm(dfr1$Y ~ dfr1$X)

# imprimimos las gráficas en un arreglo de 2x2
par(mfrow = c(2,2))
plot(modelo)