MT1021: EL MODELO DE REGRESIÓN LINEAL en R

1 Introducción y ejemplo

En el modelo de regresión lineal hay:

• Una variable numérica (que denotaremos por \(Y\)), objeto de estudio, cuyo valor depende en parte de otras. Por eso se le llama variable dependiente, respuesta, explicada, regresada, etc. (depende de la disciplina donde se aplica).
• Una o varias variables numericas (que denotaremos por \(X\) o bien \(X_1\), \(X_2\), …, \(X_p\)), que parecen influir sobre \(Y\). Su valor puede ser aleatorio o controlado por el experimentador. Por eso se les llama variables independientes, de control, explicativas, regresoras, etc.
• Unos parámetros que expresan la dependencia de \(Y\) respecto de \(X\). Se denotan por \(\beta_0\), \(\beta_1\), …, \(\beta_p \in \mathbb{R}\), y \(\sigma^2 > 0\).
• La relación de dependencia de \(Y\) respecto de las variables independientes se puede escribir con la ecuación: \(Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots \beta_p X_p + \text{N}(0, \sigma^2)\).
• Esto significa que en cada muestro, el valor de \(Y\) es la combinación lineal de los valores de las variables independientes MÁS un término de ERROR ALEATORIO, que sale siempre del mismo modelo normal \(\text{N}(0, \sigma^2)\) de forma independiente de otros muestreos y de los valores que tomen las \(X\)’s.

Ejercicio evaluable 1: Muestrear 20 datos del modelo lineal representado por la formula \(Y = 2.5 - 0.6 X + E\) con \(E \sim \text{N}(0, 0.5)\). Para ello elige \(X\) que tome 2 veces cada valor entero desde 1 hasta 10.

1. Completa el código que falta para realizar el muestreo y representarlo gráficamente. Copia y pega en un editor de texto, completa lo que falta y pásalo a la R Console para ver el resultado:

# CAMBIA EL PARÁMETRO eval DE ARRIBA A TRUE PARA COMPILAR !
xi = rep(x=1:10, each=2) # 20 muestreos de la X
ei =                     # 20 muestreos del error (normal)
yi =                     # valores de Y a partir de la fórmula del modelo
xx = data.frame(X=xi, Y=yi) # muestra conjunta (en hoja de datos)
xx # visualizamos datos
plot(xx, ylim=c(-6,3)) # muestra en gráfica
abline(a=2.5, b=-0.6) # recta "real" del modelo lineal

2. Observa la muestra en el gráfico, junto a la verdadera recta que relaciona \(Y\) con \(X\) (sin el error aleatorio).

3. Vuelve a pegar indefinidamente, sin perder de vista la ventana gráfica, para ver el efecto del muestreo. Observa cómo los datos fluctúan de manera aleatoria en torno a la recta \(Y = 2.5 - 0.6X\).

2 Representación gráfica de los datos para intuir si el modelo lineal es creíble

• Si se trata de regresión simple (una sola variable independiente X), lo lógico es representar la nube de puntos \((x_i, y_i)\) para ver si el modelo lineal es razonable.
• Si hay más variables independientes, la nube de puntos es tridimensional (o mayor), y no se aprecia bien. En ese caso puede “suponer correcto” el modelo de regresión lineal y chequear más tarde, cuando se observen las estimaciones de los errores, si la suposición ha sido aceptable o no.

3 Regresiones no lineales

El modelo de regresión lineal se puede adaptar para ajustar datos a muchas familias de funciones. Por ejemplo:

• Regresión polinómica de \(Y\) sobre \(X\): por ejemplo, para grado 3, el modelo \(Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3 + E\) aproxima datos \((x_i, y_i)\) a un polinomio de grado 3 usando regresión múltiple de \(Y\) sobre 3 variables, que en realidad es la variable original \(X\) y sus potencias, es decir \(X_1 = X\), \(X_2 = X^2\) y \(X_3 = X^3\).
• Regresión exponencial de \(Y\) sobre \(X\): transformando una o varias variables (\(\log(Y)\), \(\log(X)\)) y usando el modelo lineal, se llega a otras posibles relaciones funcionales entre ellas.

4 Estimación de parámetros del modelo de regresión lineal

Se procede a estimar los parámetros del modelo, es decir, calcular \(\widehat{\beta}_0^{LS}\), \(\widehat{\beta}_1^{LS}\), …, \(\widehat{\beta}_p^{LS}\) y \(\widehat{\sigma}^2_{ML}\) (por mínimos cuadrados y máxima verosimilitud).

• Función: lm( formula, data, subset,... )
• Argumentos:
  • data: hoja de datos con columnas,
  • formula: expresión de la dependencia usando las etiquetas de las variables: y~x1+x2+... (los puntos suspensivos indican posibles variables independientes si hay más),
  • subset: expresión logica para seleccionar la parte de los datos que entran en juego
• Devuelve: una lista con componentes:
  • $coefficients: vector de coeficientes (\(\widehat{\beta}_0^{\text{LS}}\), \(\widehat{\beta}_1^{\text{LS}}\), etc. ).
  • $residuals: los residuos \(y_i - \widehat{y}_i\) (valor real de \(Y\) menos valor que predice la recta o hiperplano de regresión para las variables independientes).
  • $fitted.values: \(\widehat{y}_i\) (la predicción de Y para cada dato de variable independiente de la muestra).

No devuelve la estimación de \(\sigma^2\), pero recuerda que es la varianza, con denominador \(n\), de los residuos. Así que se puede calcular a partir del objeto que devuelve la función lm().

Ejemplo 1: Carga el conjunto de datos con data(iris). En ellos hay datos de 150 flores sobre la longitud y anchura de pétalo y sépalo, y la especie a la que pertenece. Supongamos que la “longitud de pétalo” puede depender del resto de variables numéricas según el modelo lineal. Estima todos los parámetros del modelo lineal (recuerda que la estimación de \(\sigma^2\) es la varianza, no la cuasivarianza, de los residuos) en dos casos:

1. Usando todos los datos (guarda el objeto de la función con nombre regflores)
2. Usando sólo los datos de las flores de la especie versicolor (guarda el objeto de la función con nombre regversi)

Ejercicio evaluable 2: Usando como datos la última de las simulaciones del ejercicio anterior (variable xx), supondremos que siguen el modelo lineal (lo cual es cierto porque así se han generado los datos). Escribe las estimaciones de los parámetros del modelo y comenta lo mucho o poco que se parecen a los valores verdaderos de los parametros (ya que en este ejemplo simulado los conocemos realmente).

4 Resumen de resultados del modelo lineal: modelo ajustado, contrastes sobre los coeficientes, bondad de ajuste

• Función: summary(object)
• Argumentos:
  • object: objeto resultado de aplicar la función lm().
• Devuelve:
  • Residuos (\(y_i - \widehat{y}_i\)) y un resumen estadístico de los mismos,
  • Tabla con una fila para cada parámetro (\(\beta_0\) , \(\beta_1\), etc.) y con columnas:
    • Estimación puntual (por mínimos cuadrados),
    • Error estándar de la estimación,
    • Valor del estadístico t de Student,
    • p-valor de la estimación (para contrastar si puede ser “0” el verdadero valor del parámetro),
  • Error estándar de los residuos (con sus grados de libertad),
  • Coeficientes R-cuadrado de bondad de ajuste (el normal y el ajustado).
  • Estadístico F para contrastar si Y no depende de la(s) variable(s) independiente(s) (i.e. todos los parámetros nulos excepto \(\beta_0\)), sus grados de libertad y el p-valor del mismo.

Ejemplo 2: Calcula el coeficiente de bondad de ajuste R-cuadrado para los dos modelos lineales del Ejemplo 1, y contrasta si es posible, en ambos casos, que la “longitud de pétalo” NO dependa en absoluto de las variables que entraban en el modelo, usando un nivel de significación del 5%.

Ejercicio evaluable 3: Calcula el coeficiente R-cuadrado para el modelo ajustado en el Ejercicio evaluable 2, y contrasta si es posible que \(Y\) NO dependa de \(X\) (a pesar de que sabemos seguro que sí), usando un nivel de significación del 10%.

5 Contrastes de hipótesis sobre los valores de los coeficientes del modelo lineal (solo para contrastar con el valor “0”)

El p-valor de cada uno de esos contrastes se halla en la fila correspondiente al parámetro, que aparece en la tabla de coeficientes de la salida de la función summary(). Con dicho p-valor se debe decidir el contraste según el nivel de significación \(\alpha\) utilizado (rechazar \(H_0\) si \(\text{p-valor} < \alpha\)).

Ejemplo 3: Comenta para los modelos del Ejemplo 1, si es posible que alguno de los coeficientes sea nulo, y que la variable a la que multiplica desaparezca del modelo. Usa una significación del 10%.

Ejercicio evaluable 4: ¿Es creible que el término independiente, del modelo lineal que ajusta los datos del Ejercicio evaluable 1, sea \(0\)? Usa una significación del 5% para contestar.

6 Intervalos de confianza (IC) para los parámetros del modelo lineal

• Función: confint(object,param,level)
• Argumentos:
  • object: objeto resultado de aplicar la función lm().
  • param: vector con el nombre de las variables independientes (o etiquetas de columnas) para cuyos coeficientes se desea ver el intervalo de confianza. Si se omite se toman todas las variables independientes.
  • level: nivel de confianza para el intervalo.
• Devuelve: una matriz con los intervalos de confianza de cada coeficiente.

Ejercicio evaluable 5: Para el modelo lineal que ajusta los datos del Ejercicio evaluable 1, escribe un intervalo de confianza para el verdadero valor de cada parámetro, con un nivel de confianza del 80%.

7 Predicciones con el modelo lineal (con eventual intervalo de confianza para la predicción o para la media de las predicciones)

• Función: predict(object, newdata, interval, level)
• Argumentos:
  • object: objeto resultado de aplicar la función lm()
  • newdata: hoja de datos con los nuevos datos de la(s) variable(s) independiente(s) para los que se desea hacer predicciones.
  • interval: si se quiere obtener el IC para la media (poner confidence) o para la predicción (poner prediction
  • level: nivel de confianza para el intervalo.
• Devuelve: como mínimo las predicciones solicitadas
  • $coefficients: vector de coeficientes (\(\widehat{\beta}_0^{\text{LS}}\), \(\widehat{\beta}_1^{\text{LS}}\), etc.)
  • $residuals: los residuos \(y_i - \widehat{y}_i\).

Ejemplo 4: ¿Qué predicción de “longitud de pétalo” corresponde a una flor de especie desconocida, cuya anchura de pétalo es 1.55, longitud de sépalo 5.3 y anchura de sépalo 3.1? ¿Y si se sabe que la flor es de la especie virginica? Haz además un intervalo de confianza para cada predicción usando un nivel de confianza del 95%.

Ejercicio evaluable 6: Para el modelo lineal que ajusta los datos del Ejercicio evaluable 1, realiza una predicción y su intervalo de confianza al 80%, para un dato cuya X sea 5.5.

8 Comprobaciones gráficas de la adecuación del modelo lineal a los datos

• Función: plot( object )
• Argumentos:
  • object: objeto resultado de aplicar la función lm().
• Devuelve: una sucesión de 4 gráficas, de las que las 2 primeras se pueden interpretar como sigue:
  1. Residuos vs valores ajustados (\(y_i - \widehat{y}_i\) vs \(\widehat{y}_i\)): permiten contrastar visualmente la aleatoriedad de los errores (si se observa un patrón no aleatorio de izquierda a derecha se puede dudar de la independencia y/o de la varianza constante del error).
  2. Q-q plot normal: se puede contrastar la hipótesis de normalidad de los residuos (si se alejan demasiado de la diagonal se puede dudar de la normalidad).

Ejemplo 7: Muestra los gráficos diagnósticos con comentarios para decidir si los datos parecen cumplir las condiciones del modelo lineal ajustado para cada uno de los casos del Ejemplo 1.

Ejercicio evaluable 7: Muestra los gráficos junto a los comentarios para decidir si los datos parecen cumplir con las hipótesis del modelo lineal (y por tanto nos podemos fiar de los cálculos de los ejercicios anteriores)

9 Ejemplo sorprendente

Hay un clásico conjunto de datos en R, que se puede cargar con:

data(anscombe)

Es una hoja de datos con columnas para reflexionar sobre el peligro de uso del modelo lineal “a ciegas”. Las gráficas de las cuatro parejas de variables son:

Se observa que se trata de:

• Cuatro dependencias de \(Y\) sobre \(X\) muy distintas (patrones de las nubes de puntos),
• La misma recta de regresión y coeficiente R-cuadrado (bondad de ajuste).

Moralejas de este ejemplo:

• No confiar a ciegas en los estadísticos: simplifican demasiado los datos.
• Graficar los datos cuando se pueda, para conocer mejor lo que pasa.
• R-cuadrado alto (cerca de “1”) significa datos globalmente “cerca” de la recta, pero no indica necesariamente que el modelo lineal sea adecuado. Podria haber patrones raros en los residuos, que desacrediten su uso.

Por ejemplo, la falta de “patron” en los errores (independencia y varianza común):

descartaría por evidencias de patrón no aleatorio el uso del modelo lineal en los ejemplos 2 y 3, mientras que el 4 quedaría bajo sospecha porque no hay datos suficientes para apreciar patrones. Respecto a la normalidad de los residuos:

## Warning: not plotting observations with leverage one:
##   8

ningún gráfico levanta sospechas, porque la muestra es tan pequeña que lo normal puede oscilar mucho. Únicamente el ejemplo 3 un residuo se aleja mucho de lo normal, pero es sólo uno (puede ser error de medida al tomar el dato o no).

10 Ejercicios

Ejercicio evaluable 8: Abre una sesión de R e introduce la línea

data(women)

• Inspecciona la variable women, que almacena alturas y pesos de mujeres americanas. En realidad no son datos de mujeres, sino promedios de pesos de mujeres con alturas determinadas, pero no se indica de cuántas mujeres por cada altura.
• Para poder usarlo como ejemplo, vamos a suponer que cada dato es la altura y peso de una persona. De esta manera:
• Realiza el gráfico de peso vs altura e indica si la relación lineal parece razonable o sospechosamente descabellada.
• Escribe la ecuación de la recta de regresión que ajusta los datos.
• Suponiendo correcto el modelo lineal, ¿se puede rechazar, estadísticamente y usando una significación del 5%, que el peso depende de la altura?
• Escribe los intervalos de confianza al 99% para los verdaderos valores de la ordenada en el origen (intercept) y la pendiente.
• Suponiendo correcto el modelo lineal, ¿se puede rechazar, estadísticamente y usando una significación del 95%, que la verdadera ordenada en el origen (intercept) vale 0?
• Suponiendo correcto el modelo lineal, ¿se puede rechazar, estadísticamente y usando una significación del 95%, que la ordenada en el origen vale -100? (Ayuda: haz el IC.)
• Aporta un indicador de la bondad del ajuste de los datos a la recta e interprétalo.
• Pronostica el peso (en libras) de una mujer de 74 pulgadas, y acompáñalo con un intervalo de confianza al 80%.
• ¿Los datos son incompatibles con las hipótesis del modelo de regresión lineal? Usa los dos gráficos habituales, y realiza también el gráfico de los residuos vs la variable independiente. Interprétalos todos.

Ejercicio evaluable 9: Supongamos que se pretende analizar el efecto de un medicamento sobre el nivel de colesterol en sangre. Se experimenta con pacientes de características muy similares, administrando distintas dosis del medicamento a cada uno. Usa el modelo de regresión lineal múltiple para ajustar los datos del fichero colesterol.txt a un modelo polinómico de hasta grado 4.

• Representa la nube de puntos de los datos.
• Escribe la ecuación del polinomio de regresión de grado 4 por mínimos cuadrados.
• Según los datos y usando una significación del 5%, ¿se puede admitir que el polinomio de la regresión no alcanza el grado máximo, y que por tanto se puede buscar un grado menor? ¿Por qué?
• Repite el apartado 1, bajando el grado del polinomio, mientras el apartado 2 siga teniendo respuesta afirmativa.
• Una vez encontrado el grado para el que no hay argumentos estadísticos que permitan “bajar”, realiza la predicción del posible valor de Y para X = “últimas dos cifras de tu DNI”, y escribe el intervalo de confianza al 90% para dicha predicción.
• Copia las gráficas que permiten evaluar la adecuación del modelo de regresión lineal, e interprétalas.

Ejercicio evaluable 10: Define una función en lenguaje R, con las siguientes características:

• Que se llame modeloLineal.
• Que admita como argumentos:
  • x: una matriz u hoja de datos cuyas columnas se puede interpretar como muestras de variables independientes.
  • beta: un vector de dimensión “número de columnas de x más 1”, que será el vector de parámetros del modelo lineal.
  • sigma2: un valor no negativo, que será el parámetro \(\sigma^2\) del modelo.
  • seed: un valor entero que sirva de semilla para poder replicar las simulaciones que calcule.
• Que devuelva una hoja de datos con todas las columnas de x a la que se añade una columna nueva, la de los valores simulados de la \(Y\) (calculados a partir de las x y los parámetros del modelo).

En definitiva, un ejemplo de invocación de esa función puede ser

• modeloLineal(x=1:10, beta=c(2,0.5), sigma2=1.5, seed=1000)

# Escribe aquí tu codigo y compila con RStudio CTRL+SHIFT+K
# verás la salida de R