Regresión lineal múltiple

Comparamos:

Regresión lineal simple (una variable independiente o regresora)

$$ y = \beta_{0} + \beta_{1}x+\varepsilon $$

Regresión lineal múltiple ($n$ variables independientes o regresoras)

$$ y = \beta_{0} + \beta_{1}x_{1}+ \beta_{2}x_{2}+\cdots+ \beta_{n}x_{n}+\varepsilon $$

---

La implementación de una regresión lineal múltiple es muy sencilla:

df <- iris
head(df)

A data.frame: 6 × 5

	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
	<dbl>	<dbl>	<dbl>	<dbl>	<fct>
1	5.1	3.5	1.4	0.2	setosa
2	4.9	3.0	1.4	0.2	setosa
3	4.7	3.2	1.3	0.2	setosa
4	4.6	3.1	1.5	0.2	setosa
5	5.0	3.6	1.4	0.2	setosa
6	5.4	3.9	1.7	0.4	setosa

# Definimos las variables

# Variables regresoras
x1 <- df$Sepal.Length
x2 <- df$Sepal.Width
x3 <- df$Petal.Width

# Variable objetivo
y <- df$Petal.Length

# implementacion
fit <- lm(y ~ x1 + x2 + x3)
summary(fit)

Call:
lm(formula = y ~ x1 + x2 + x3)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.99333 -0.17656 -0.01004  0.18558  1.06909 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.26271    0.29741  -0.883    0.379    
x1           0.72914    0.05832  12.502   <2e-16 ***
x2          -0.64601    0.06850  -9.431   <2e-16 ***
x3           1.44679    0.06761  21.399   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.319 on 146 degrees of freedom
Multiple R-squared:  0.968,	Adjusted R-squared:  0.9674 
F-statistic:  1473 on 3 and 146 DF,  p-value: < 2.2e-16

---

Análisis de supuestos sobre una regresión lineal simple

# Variable regresora
x <- df$Petal.Width

# Variable objetivo
y <- df$Petal.Length

# implementacion
fit <- lm(y ~ x)
summary(fit)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33542 -0.30347 -0.02955  0.25776  1.39453 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08356    0.07297   14.85   <2e-16 ***
x            2.22994    0.05140   43.39   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4782 on 148 degrees of freedom
Multiple R-squared:  0.9271,	Adjusted R-squared:  0.9266 
F-statistic:  1882 on 1 and 148 DF,  p-value: < 2.2e-16

plot(x=x, y=y)

1. Homocedasticidad: Emplearemos la prueba de hipótesis Breuch-Pagan

• $H_{0}$: varianza constante (homocedasticidad) vs
• $H_{1}$: no hay varianza constante.

Buscamos no rechazar ($p-value>\alpha$).

• p-value > 0.05: nos quedamos quedamos con $H_{0}$
• p-value < 0.05: nos quedamos quedamos con $H_{1}$

# Librería necesaria
library(lmtest)

# bptest(fit): breush-pagan
bptest(fit)

Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

	studentized Breusch-Pagan test

data:  fit
BP = 23.897, df = 1, p-value = 1.016e-06

dado que

$$ p-value = 1.016e-06 = 0.000001016<\alpha = 0.05 $$$$ 1.016e-06=1.016x10^{-6} $$

entonces se rechaza $H_{0}$, es decir, se rechaza la homocedasticidad por lo cual este supuesto no se está cumpliendo.

1. Linealidad: las variables interactúan mediante una relación lineal.

# buscamos que la recta roja sea lo mas parecida a 
# la recta horizontal
plot(fit, 1)

El supuesto de la linealidad tampoco se está cumpliendo

1. Los errores (valor predicho - valor real) se distribuyen normal. Para comprobar empleamos la prueba de hipótesis Shapiro-Wilk. Buscamos no rechazar $H_{0}$ (los errores se distribuyen normal).

p-value > 0.05: nos quedamos con $H_{0}$

# Utilizamos la librería broom para calcular los errores de forma automática
library(broom)

# Guardamos los errores en columnas de nuestro dataframe de los datos 
Datosfit <- augment(fit)

shapiro.test(Datosfit$.std.resid)

	Shapiro-Wilk normality test

data:  Datosfit$.std.resid
W = 0.99019, p-value = 0.3823

Este supuesto sí se está cumpliendo

1. Las observaciones son independientes entre sí. Para ello emplearemos la prueba de Durbin-Watson donde buscamos no rechazar $H_{0}$ (las observaciones son independientes entre sí)

p-value > 0.05: nos quedamos con $H_{0}$

#Prueba para autocorrelación de orden 1

dwtest(fit)

	Durbin-Watson test

data:  fit
DW = 1.4296, p-value = 0.0001623
alternative hypothesis: true autocorrelation is greater than 0

p -value < 0.05: nos quedamos con $H_{1}$

library(randtests)
runs.test(Datosfit$.std.resid)

	Runs Test

data:  Datosfit$.std.resid
statistic = -3.7664, runs = 50, n1 = 74, n2 = 69, n = 143, p-value =
0.0001656
alternative hypothesis: nonrandomness

#Prueba de rachas
library(lawstat)
runs.test(Datosfit$.std.resid, plot.it = TRUE)

Attaching package: ‘lawstat’


The following object is masked from ‘package:randtests’:

    runs.test

	Runs Test - Two sided

data:  Datosfit$.std.resid
Standardized Runs Statistic = -4.2601, p-value = 2.043e-05

Este supuesto tampoco se cumple.