Introducción

Este documento presenta una guía sobre la ejecución de la regresión lineal múltiple con el software R. El análisis de datos multifactorial sirve para expresar la relación entre una variable de interés (la respuesta) y un conjunto de variables predictoras relacionadas.

Considere ahora un ejemplo hecho con el software R que detalla paso a paso las instrucciones para realizar la regresión lineal múltiple:

Lectura de la base de datos

La base de datos está alojada en Github, para leerla se usa el código siguiente:

datos <- read.csv("https://raw.githubusercontent.com/fhernanb/Python-para-estadistica/master/03%20Regression/Regresi%C3%B3n%20lineal%20simple/softdrink.csv", header = T)

Observar la base de datos

En este paso, se ven las variables, las primeras observaciones y la estructura de la base de datos, entre otras cosas.

names(datos) # Ver nombre de las variables
## [1] "Obs" "y"   "x1"  "x2"
head(datos)  # Primeras observaciones
str(datos)   # Estructura de la base
## 'data.frame':    25 obs. of  4 variables:
##  $ Obs: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ y  : num  16.7 11.5 12 14.9 13.8 ...
##  $ x1 : int  7 3 3 4 6 7 2 7 30 5 ...
##  $ x2 : int  560 220 340 80 150 330 110 210 1460 605 ...

Note que el software R leyó la columna Obs de la base de datos como una variable, para eliminarla, se hace lo siguiente:

datos1 <- datos[,-1] # Elimina la primera columna
names(datos1)
## [1] "y"  "x1" "x2"

Coeficiente de correlación

Por consiguiente, se calcula el coeficiente de correlación de Pearson entre las variables regresoras.

cor(datos1$x1, datos1$x2)
## [1] 0.824215
  • Se puede observar que las variables regresoras (x1 y x2) están altamente correlacionadas.

Matriz de correlación

Para calcular dicha correlación, se emplea el siguiente código:

round(cor(datos1), 3) # Matriz de correlación redondeada a 3 decimales
##        y    x1    x2
## y  1.000 0.965 0.892
## x1 0.965 1.000 0.824
## x2 0.892 0.824 1.000

Matriz de dispersión

Esta matriz es una herramienta de exploración de datos que permite comparar varios subconjuntos de la base de datos para buscar patrones y relaciones. Se obtiene así:

# Opción básica
plot(datos1)

# Opción mejorada 
pairs(datos1, labels=c("y","x1","x2"), main='Matriz de dispersión', cex.main=0.8, cex = 1.5, pch = 20, bg="light blue", cex.labels = 1, font.labels=1)

Matriz de dispersión con correlaciones

Esta matriz también sirve para explorar los datos al comparar varios subconjuntos de la base de datos, pero adiciona la correlación entre cada uno de los subconjuntos. Para obtenerla con R, en este caso se creo una función y luego se usó la función pairs(), obsérvelo a continuación:

panel.cor <- function(x1, y, digits=2, prefix="", cex.cor){
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x1, y))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste(prefix, txt, sep="")
  if(missing(cex.cor))
  cex <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex)
}

pairs(datos1, lower.panel=panel.smooth, upper.panel=panel.cor, 
      main="Matriz de dispersión con correlaciones")

Ajustar el modelo de regresión

Ahora, se ajusta el modelo de regresión lineal para los datos en cuestión usando la función lm().

mod <- lm(y ~ x1 + x2, data = datos1)

Resultados del modelo ajustado

Para obtener un resumen del modelo ajustado se usa:

# Resumen del modelo ajustado
summary(mod)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = datos1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7880 -0.6629  0.4364  1.1566  7.4197 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.341231   1.096730   2.135 0.044170 *  
## x1          1.615907   0.170735   9.464 3.25e-09 ***
## x2          0.014385   0.003613   3.981 0.000631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.259 on 22 degrees of freedom
## Multiple R-squared:  0.9596, Adjusted R-squared:  0.9559 
## F-statistic: 261.2 on 2 and 22 DF,  p-value: 4.687e-16

  • Residuals: son las diferencias entre las respuestas observadas en cada combinación de valores de las variables regresoras y la correspondiente predicción de la respuesta calculada.

  • Estimate: corresponde a las estimaciones de los coeficientes del modelo. Los coeficientes de las variables regresoras siempre tienen interpretación; en este caso, el coeficiente de x1 expresa el efecto sobre la variable respuesta al controlar o ajustar la regresora x2; es decir, se asocia un aumento de una unidad en la variable regresora x1 con un aumento de 1.615907 en el ajuste de la variable respuesta al controlar la varaible regresora x2. En cambio, el intercepto no siempre puede ser interpretado, dependerá exclusivamente de la naturaleza del problema.

  • Std Error: representa el error estándar de cada coeficiente del modelo.

  • t value: es el valor que toma el estadístico de prueba.

  • Pr(>|t|): corresponde al valor-p para las pruebas de hipótesis de que los coeficientes son ceros. Si esta valor-p es menor que el nivel de significancia, que por defecto es del 5%, entonces se rechaza la hipótesis nula y por lo tanto, existe significancia en la regresión. A menudo, la prueba de hipótesis de que el intecerpto es cero es poco útil, ya que éste depende de la naturaleza del problema. En cambio, las pruebas de hipótesis de que los coeficientes de las variables regresoras son cero, son muy importantes; debido a que, si se prueba que el coeficiente es cero, entonces la variable regresora no tendría ningún aporte al modelo.

  • Signif codes: representa el nivel de significancia de las variables regresoras para la predicción del modelo de regresión. Entre más estrellas se le asigne a la regresora, más contribuye con el modelo.

  • Residual standard error: es la estimación de la desviación estandar. Este valor da una idea de cuán lejos están los valores del modelo ajustado a los valores observados de la variable respuesta.

  • Multiple R-squared: este valor expresa la proporción de la variación explicada por el modelo; es decir, por las variables explicativas. En este caso, aproximadamente el 96% de la variación en la variable respuesta puede ser explicada por el modelo; es decir, por las variables regresoras (x1 y x2).

  • F-statistic: equivale al estadístico F. Éste permite constrastar la hipótesis nula de que el valor poblacional R es cero, lo cual, equivale a contrastar la hipótesis de que los coeficientes de las variables regresoras son cero.

  • p-value: índica la importancia del modelo mediante una prueba de hipótesis donde la hipótesis nula es que todos los coeficientes del modelo son cero. En este ejemplo, prueba específicamente que los coeficientes de las regresoras (x1 y x2) son cero.

Similarmente, las estimaciones de los parámetros del modelo se obtienen al usar el siguiente cógido:

# Coeficientes
coef(mod)               # Opción 1
## (Intercept)          x1          x2 
##  2.34123115  1.61590721  0.01438483
mod$coef                # Opción 2
## (Intercept)          x1          x2 
##  2.34123115  1.61590721  0.01438483
# Varianza estimada
sigma(mod)^2            # Opción 1 
## [1] 10.62417
summary(mod)$sigma^2    # Opción 2
## [1] 10.62417
# Desviación estándar 
sigma(mod)              # Opción 1
## [1] 3.259473
summary(mod)$sigma      # Opción 2
## [1] 3.259473

Intervalos de confianza para los coeficientes del modelo

Se usará el comando confint() para hallar dichos intervalos, así:

confint(mod, conf.level=0.95) # I.C al 95% 
##                   2.5 %     97.5 %
## (Intercept) 0.066751987 4.61571030
## x1          1.261824662 1.96998976
## x2          0.006891745 0.02187791

  • Con una confianza del 95% se puede decir que el verdadero valor del coeficiente de la regresora x1 está entre 1.26 y 1.97. Observe que la estimación de dicho coeficiente fue 1.62.

  • Con una confianza del 95% se puede decir que el verdadero valor del coeficiente de la regresora x2 está entre 0.007 y 0.022. Observe que la estimación de dicho coeficiente fue 0.014.

Gráfico de dispersión 3D

Para realizar este gráfico es necesario instalar y cargar la libreria plot3D y plot3Drgl, las cuales sirven simultanéamente para dibujar un diagrama de dispersión en tres dimensiones y brindar interacción al gráfico; es decir, rotación o movimiento.

# install.packages("plot3D")
# install.packages("plot3Drgl")
library(plot3D)
library(plot3Drgl)

Luego, se obtiene el diagrama de dispersión en 3D.

grid.lines <- 40
x1.pred <- seq(min(datos1$x1), max(datos1$x1), length.out = grid.lines)
x2.pred <- seq(min(datos1$x2), max(datos1$x2), length.out = grid.lines)
x1x2 <- expand.grid(x1 = x1.pred, x2 = x2.pred)
y.pred <- matrix(predict(mod, newdata = x1x2), nrow = grid.lines, ncol = grid.lines)
fitpoints <- predict(mod)
scatter3D(x=datos1$x1, y=datos1$x2, z=datos1$y, theta=30, phi=8, pch=20)

Gráfico de dispersión con plano de regresión

grid.lines <- 40
x1.pred <- seq(min(datos1$x1), max(datos1$x1), length.out = grid.lines)
x2.pred <- seq(min(datos1$x2), max(datos1$x2), length.out = grid.lines)
x1x2 <- expand.grid(x1 = x1.pred, x2 = x2.pred)
y.pred <- matrix(predict(mod, newdata = x1x2), nrow = grid.lines, ncol = grid.lines)
fitpoints <- predict(mod)

scatter3D(x=datos1$x1, y=datos1$x2, z=datos1$y, theta=30, phi=8, pch=20, bty = "g", 
          colkey = TRUE, surf = list(x=x1.pred, y=x2.pred, z=y.pred, fit = fitpoints))

Nota: si desea ver la interacción de los gráfico anteriores, use el comando plotrgl() al final del código de cada gráfico de dispersión.

Fuente

Documento fue creado durante el curso Análisis de Regresión de la Universidad Nacional de Colombia sede Medellín. La base de datos utilizada son del ejemplo 3.1 del libro de Montgomery, Peck and Vining (2003). Además, dicha base se encuentra alojada en la cuenta de Github del profesor Freddy Hernández Barajas, https://github.com/fhernanb. Adicionalmente, se le agradece al compañero Felipe Bedoya por brindar aportes con algunos gráficos.