Regresión Lineal Simple y Múltiple
Regresión Lineal Simple
La regresión es una técnica estadística utilizada para encontrar una relación numérica entre dos variables. Entonces, se trata de hallar el parámetro que explique tal relación. Aquí, se desarrollan modelos de regresión lineal simple y múltiple. Posteriormente se explicarán las consideraciones de significancia y bondad de ajuste. Este tipo de regresión solo maneja una variable independiente, a causa de esto el modelo por estimar únicamente tendrá dos parámetros: \[ Y = B_o+B_1X+e... (1) \] Donde:
- \(X\) es la variable independiente.
- \(Y\) es la variable dependiente.
- \(e\) es el error de estimación.
- \(B_o\) es el intercepto u ordenada en el origen.
- \(B_1\) es la pendiente de la regresión.
Entonces, la esperanza de (1) es: \[ E(Y)= \widehat{Y} = E(\beta_o) + E(\beta_1X) + E(e) \\ \widehat{Y} = \beta_0 + \beta_1X \] Por el método de mínimos cuadrados ordinarios: \[ \widehat{\beta_1}=\frac{n\sum XY - \sum X \sum Y}{n\sum X^2-(\sum X)^2} \\ \widehat{\beta_0}=\frac{\sum Y \sum X^2 - \sum X \sum XY}{n \sum X^2 - (\sum X)^2} \] Luego, \(\widehat{y}=\widehat{\beta_0}+\widehat{\beta_1X}\) es la recta de regresión lineal de “y” sobre “x”.
Función LM de R
Para estimar un modelo de regresión lineal simple en R, debe usarse el comando “lm()”. \[ reg.simple<-lm(ingresos~empleo, data=bankloan) \] Donde:
- ingresos es la variable dependiente.
- empleo es la variable independiente.
- bankloan es la data usada.
- Reg.simple es el modelo de regresión lineal simple.
Ejemplo
Una empresa requiere saber si existe una relación confiable entre las inversiones en publicidad y las ventas que se obtienen, a fin de calcular estimar las ventas que recibirían si hicieran una inversión en publicidad de S/ 8000, para el onceavo periodo de producción. Aproximar con el modelo de regresión lineal.
Datos
Solución:
Para estimar el modelo de regresión lineal, hay que reconocer la variable independiente y la variable dependiente:
- Variable Independiente X: Inversión en publicidad.
- Variable dependiente Y: Ventas
Resolviendo para \(\widehat{\beta_1}\): \[ \widehat{\beta_1}=\frac{10\sum_{i=1}^{10}X_iY_i-\sum_{i=1}^{10}X_i*\sum_{i=1}^{10}Y_i}{10 \sum_{i=1}^{10}X_1^2-(\sum_{i=1}^{10}X_i)^2} \\ \widehat{\beta_1}=\frac{10(X_1Y_1+...+X_{10}Y_{10})-(X_1+...+X_{10})(Y_1+...+Y_{10})}{10(X_1^2+...+X_{10}^2)-(X_1+...+X_{10})^2} \\ \widehat{\beta_1}=\frac{10*10893520600-55800*1862480}{10*326625800-55800^2} = \frac{5008822000}{152618000} = 32.8193 \] Resolviendo para \(\widehat{\beta_0}\): \[ \widehat{\beta_0}=\frac{\sum_{i=1}^{10}Y_i*\sum_{i=1}^{10}X_i^2-\sum_{i=1}^{10}X_i*\sum_{i=1}^{10}X_iY_i}{10\sum_{i=1}^{10}X_i^2-(\sum_{i=1}^{10}X_i)^2} \\ \widehat{\beta_0}=\frac{(Y_1+...+Y_{10})*(X_1^2+...+X_{10}^2)-(X_1+...+X_{10})(X_1Y_1+...+X_{10}Y_{10})}{10(X_1^2+...+X_{10}^2)-(X_1+...+X_{10})^2} \\ \widehat{\beta_0}=\frac{1862480*326625800-55800*10893520600}{10*326625800-55800^2}=\frac{475570504000}{152618000}=3116.0839 \\ \] Hallado el B1 y el B0, la ecuación de la regresión quedaría de la siguiente forma: \[ \widehat{y}=3116.0840+32.8193X \] Se quiere evaluar cuando la inversión sea S/ 8000.00 \[ 3116.0840+32.8193(8000)=265670.4840 \] Entonces, para una inversión de S/ 8000.00 en el onceavo período, se tendrá una venta aproximada de S/ 265,670.4840.
En R, el ejemplo sería:
Inversion<-c(5000,5570,4350,7900,6800,5400,6900,3900,4200,5780)
Venta<-c(160000,189380,139200,260700,217600,183600,234600,136500,138600,202300)
Data<-data.frame(Inversion,Venta)
Data
## Inversion Venta
## 1 5000 160000
## 2 5570 189380
## 3 4350 139200
## 4 7900 260700
## 5 6800 217600
## 6 5400 183600
## 7 6900 234600
## 8 3900 136500
## 9 4200 138600
## 10 5780 202300
Hacemos uso de la función “lm”:
##
## Call:
## lm(formula = Venta ~ Inversion, data = Data)
##
## Coefficients:
## (Intercept) Inversion
## 3116.08 32.82
De una forma más detallada:
##
## Call:
## lm(formula = Venta ~ Inversion, data = Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8687.6 -5599.5 785.3 4637.9 9488.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3116.084 9641.674 0.323 0.755
## Inversion 32.819 1.687 19.454 5.06e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6591 on 8 degrees of freedom
## Multiple R-squared: 0.9793, Adjusted R-squared: 0.9767
## F-statistic: 378.4 on 1 and 8 DF, p-value: 5.064e-08
Con la inversión de 8000:
beta0<-reg.simple$coefficients[1]
beta1<-reg.simple$coefficients[2]
paste('El valor de B0 es:', beta0)
## [1] "El valor de B0 es: 3116.08397436736"
## [1] "El valor de B1 es: 32.8193397895399"
X<-8000
paste('Se espera que la venta para una inversión de S/',X,"sea de S/", format(beta0+X*beta1,digits = 10))
## [1] "Se espera que la venta para una inversión de S/ 8000 sea de S/ 265670.8023"
Para ver la gráfica de dicha regresión:
plot(Inversion, Venta, pch=22, bg="red",
main="Venta VS Inversión",xlab="Inversión en publicación (S/.)",ylab="Venta (S/.)")
abline(reg=reg.simple, col="blue", lwd="1.5")
recta<-reg.simple$coefficients
text(6800, 250000, bquote(paste('y=', .(format(recta[1],digits=6)),
'+',.(format(recta[2],digits=4)),'x')))
grid(10,10,lwd=2)
Para ver los residuos de esta regresión de manera gráfica:
plot(predict(reg.simple),residuals(reg.simple), pch= 21, bg="green",
main = "Gráfico de Residuos", xlab="Predicción (S/.)", ylab="Residuo del Modelo (S/.)")
grid(10, 10, lwd=15)
lines(predict(reg.simple), residuals(reg.simple),type = "h", col="green", lwd=2)
abline(h=0, col="red", lty=3, lwd=1.5)
Regresión Lineal Múltiple
Sea 𝑌 una variable de respuesta cuantitativa, y al menos una variable de predicción \(x_i\) es cuantitativa. Para estos casos, el modelo de regresión lineal múltiple suele ser muy útil: \[ y=\beta_0+\beta_1x_1+\beta_2x_2+...+\beta_kx_k+e \] Donde:
- \(y\) es el regresando.
- \(x_1,x_2,x_k\) son los regresores.
- \(e\) reprenta el error o perturbación aleatoria.
- Los parámetros \(\beta_0,\beta_1...,\beta_k\) son fijos y desconocidos.
- Supóngase que se cuenta con una muestra aleatoria de \(n\) datos:
\[ y_1=\beta_0+\beta_1x_{1,1}+\beta_2x_{2,1}+...+\beta_kx_{k,1}+e_1 \\ y_2=\beta_0+\beta_1x_{1,2}+\beta_{2,2}+...+\beta_kx_{k,2}+e_2 \\ \vdots \\ y_n=\beta_0+\beta_1x_{1,n-1}+\beta_2x_{2,n-1}+...+\beta_kx_{k,n-1}+e_n \] Llevando ello a su forma matricial, se obtiene:
Todo conjuntamente en notación matricial sería:
Expresado de forma compacta: \[ y=X\beta+e \] Donde:
- \(y\) es un vector nx1
- \(X\) es una matriz n*k
- \(\beta\) es un vector (k+1)*1
- \(e\) es un vector n*1
Por el método de mínimos cuadrados: \[
X'X\widehat{\beta}=X'y \\
\] Lo que da:
Ejemplo:
Un estudio analiza la relación entre la composición de un cemento tipo Portland y el calor desprendido durante la fase de fraguado. La muestra está formada por 13 cementos. Aproximar una función 𝒚 que dependa de las variables 𝒙𝟏, 𝒙𝟐, 𝒙𝟑, 𝒙𝟒.
Prueba de hipótesis:
- 𝐻0: 𝛽𝑖 = 0, 𝑖 = 0,1,2,3,4
- 𝐻1: 𝛽𝑖 ≠ 0, para al menos un 𝑖
La matriz X serí: Dicha matriz inversa es:
Mientras tanto, la matriz “y” es:
Para hallar \([X'X]^{-1}\), debe demostrarse que X’X es inversible, hallando la determinante. \[ det(X'X)=27717211679.97 \not = 0 \]
Por lo tanto X’X sería:
Y la inversa de esa matriz, es decir, de \([X'X]^{-1}\):
Mientras que el cálculo del Beta se hace de la siguiente forma:
Dando como resultado:
En R sería así:
x1<-c(7,1,11,11,7,11,3,1,2,21,1,11,10)
x2<-c(26,29,56,31,52,55,71,31,54,47,40,66,68)
x3<-c(6,15,8,8,6,9,17,22,18,4,23,9,8)
x4<-c(60,52,20,47,33,22,6,44,22,26,34,12,12)
y<-c(78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3,109.4)
Data<-data.frame(x1,x2,x3,x4,y)
Data
## x1 x2 x3 x4 y
## 1 7 26 6 60 78.5
## 2 1 29 15 52 74.3
## 3 11 56 8 20 104.3
## 4 11 31 8 47 87.6
## 5 7 52 6 33 95.9
## 6 11 55 9 22 109.2
## 7 3 71 17 6 102.7
## 8 1 31 22 44 72.5
## 9 2 54 18 22 93.1
## 10 21 47 4 26 115.9
## 11 1 40 23 34 83.8
## 12 11 66 9 12 113.3
## 13 10 68 8 12 109.4
Hacemos la regresión:
##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1750 -1.6709 0.2508 1.3783 3.9254
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.4054 70.0710 0.891 0.3991
## x1 1.5511 0.7448 2.083 0.0708 .
## x2 0.5102 0.7238 0.705 0.5009
## x3 0.1019 0.7547 0.135 0.8959
## x4 -0.1441 0.7091 -0.203 0.8441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.446 on 8 degrees of freedom
## Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736
## F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07
Vemos los residuos:
plot(predict(reg.multiple),residuals(reg.multiple),pch=21, bg="red",
main="Gráfica de Residuos", xlab="Predicción", ylab="Residuo del Modelo")
grid(10, 10, lwd=1.5)
lines(predict(reg.multiple),residuals(reg.multiple), type="h", col="red", lwd=2)
abline( h=0, col="blue", lty=3, lwd=1.5)
Significancia y Bondad de Ajuste
La hipótesis nula 𝐻0 indica que la variable de respuesta no está influenciada por las variables. La hipótesis alternativa 𝐻1 indica que existe algún tipo de influencia de estas variables independientes. Si la significación es pequeña, el intervalo de confianza no contiene el valor cero. Esto se considera como una indicación de que esa variable es interesante en el modelo.
Si contiene al cero (no significativa), posiblemente sea preferible eliminar del modelo para simplificar. Sin embargo, se debe tener cuidado con las variables confusoras, aquellas que si se eliminan y los otros coeficientes cambian notoriamente siguen siendo significativas. Este tipo de variables deben conservase, aunque sus coeficientes sean significativos.
El término R cuadrado es el porcentaje de reducción de la incertidumbre cuando son conocidas las variables independientes. Tiende a estimar de forma optimista el ajuste de la regresión lineal. Siempre aumenta como el número de efectos que se incluyen en el modelo. El término R cuadrado ajustado intenta corregir esta estimación excesiva (𝑅2 > 𝑅2 ajustado).
Supuestos del Modelo de Regresión Lineal
Los supuestos que se asumen a la hora de realizar un modelo de regresión lineal son los siguientes:
Linealidad
Si no hay linealidad, se dice que hay un error de especificación.
Independencia
Los datos de las variables explicativas deben ser independientes de los residuos (principalmente para datos de series de tiempo), si no, se generará un problema de autocorrelación.
Homocedasticidad
Este supuesto afirma que debe existir igual varianza entre los residuos y los pronósticos. Implica que la variación de los residuos sea uniforme en todo el rango de valores de los pronósticos.
Normalidad
Este supuesto afirma que los residuos deben seguir una distribución normal.
No Colinealidad
No debe existir colinealidad. Esta puede ser:
Colinealidad perfecta: Si una de las variables independientes tiene una relación lineal con otras independientes.
Colinealidad parcial: Si entre variables independientes existen altas correlaciones.
Ejemplos
Ejemplo 1: Regresión Lineal Simple por Mínimos Cuadrado I
Se determinará el costo de una jornada laboral de 44 horas utilizando la siguiente información:
Introducimos los datos en R:
x<-c(10,15,18,20,24,30,35,40,45,48)
y<-c(500,800,850,950,1100,1300,1550,1700,1900,2100)
Data<-data.frame(x,y)
Data
## x y
## 1 10 500
## 2 15 800
## 3 18 850
## 4 20 950
## 5 24 1100
## 6 30 1300
## 7 35 1550
## 8 40 1700
## 9 45 1900
## 10 48 2100
Realizamos la regresión
##
## Call:
## lm(formula = y ~ x, data = Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.576 -31.846 -1.855 15.680 61.661
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 142.0495 29.9693 4.74 0.00146 **
## x 39.7527 0.9633 41.27 1.31e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38 on 8 degrees of freedom
## Multiple R-squared: 0.9953, Adjusted R-squared: 0.9947
## F-statistic: 1703 on 1 and 8 DF, p-value: 1.309e-10
La estimación indica que, con el R^2 de 0.9953 quiere decir que la variable dependiente está explicada al 99.63% por la variable independiente, horas laboradas. Además, por cad ahora laboral, el costo aumentaría en 39.75 (el beta de la variable X). Por otro lado, vemos que las variables son estadísticamente significativas.
Se busca evaluar para cuando las horas laborables son 44: \[ 142.0495+39.75(44)=1891.05 \] Entonces, el costo de las 44 horas laborales es S/ 1891.05. De manera gráfica sería como lo elaborado anteriormente, y le agregaremos intervalos de confianza.
nuevo.x<-data.frame(x=seq(0,50))
plot(Data$x, Data$y, pch=22, bg="red", main="Costo de Jornada vs Horas Laboradas",
xlab="Horas laborales (h)", ylab="Costo de hora laboral (S/.)")
abline(reg1, col="blue", lwd=1.5)
#Intervalos de confianza
int.confianza<- predict(reg1, nuevo.x, interval='confidence')
lines(nuevo.x$x, int.confianza[,2], lty=2, col='darkgreen')
lines(nuevo.x$x, int.confianza[,3], lty=2, col='darkgreen')
int.confianza2<- predict(reg1, nuevo.x, interval='prediction')
lines(nuevo.x$x, int.confianza2[,2], lty=2, col='red')
lines(nuevo.x$x, int.confianza2[,3], lty=2, col='red')
Gráfico de residuos:
plot(predict(reg1),residuals(reg1), pch=21, bg='green',
main="Gráfico de residuos", xlab="Predicción", ylab="Residuo del Modelo")
grid(10 ,10, lwd=1.5)
lines(predict(reg1), residuals(reg1), type="h", col="green", lwd=2)
abline(h =0, col="red", lty=3, lwd=1.5)
Prueba de Normalidad
Ejemplo 3: Regresión Lineal y Prueba de Hipótesis
Se tiene una base de datos que muestra información de 935 individuos. Con ello se estimará una ecuación de salarios y se aplicará una prueba de hipótesis para el valor de cierta variable explicativa en el modelo.
Los datos:
## # A tibble: 935 x 12
## Salario.Mes Horas.Sem Coef.Int A.Educ A.Exper A.Empl Edad Casado Num.Herm
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 769 40 93 12 11 2 31 1 1
## 2 808 50 119 18 11 16 37 1 1
## 3 825 40 108 14 11 9 33 1 1
## 4 650 40 96 12 13 7 32 1 4
## 5 562 40 74 11 14 5 34 1 10
## 6 1400 40 116 16 14 2 35 1 1
## 7 600 40 91 10 13 0 30 0 1
## 8 1081 40 114 18 8 14 38 1 2
## 9 1154 45 111 15 13 1 36 1 2
## 10 1000 40 95 12 16 16 36 1 1
## # ... with 925 more rows, and 3 more variables: Ord.Nacim <chr>,
## # Edu.Madre <chr>, Edu.Padre <chr>
Con ello, se hace el modelo de regresión lineal:
modelo<-lm(log(Salario.Mes)~ Horas.Sem + Coef.Int + A.Educ + A.Exper+ Edad, data= MiData)
summary(modelo)
##
## Call:
## lm(formula = log(Salario.Mes) ~ Horas.Sem + Coef.Int + A.Educ +
## A.Exper + Edad, data = MiData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.87570 -0.22024 0.01763 0.26211 1.25343
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0349059 0.1778614 28.308 < 2e-16 ***
## Horas.Sem -0.0047976 0.0017460 -2.748 0.00612 **
## Coef.Int 0.0060561 0.0009745 6.214 7.77e-10 ***
## A.Educ 0.0516386 0.0075798 6.813 1.72e-11 ***
## A.Exper 0.0126788 0.0038632 3.282 0.00107 **
## Edad 0.0150958 0.0048458 3.115 0.00189 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3831 on 929 degrees of freedom
## Multiple R-squared: 0.1768, Adjusted R-squared: 0.1724
## F-statistic: 39.91 on 5 and 929 DF, p-value: < 2.2e-16
La interpretación se los resultados de la regresión sería lo siguiente:
- Estos resultados muestran todos los coeficientes significativos de forma individual, excepto el que acompaña a la variable Horas.Sem, promedio de horas trabajadas durante el mes.
- Este modelo tiene una particularidad: es lineal. A excepción de la variable dependiente log (Salario.Mes), se ha transformado la variable a logaritmos para que se suavice la dispersión de los datos, entonces la interpretación de los datos será diferente.
- Según lo expuesto en el párrafo anterior, se entenderá en referencia a los resultados el aumento en una unidad del coeficiente intelectual (variable Coef.Int), que generará un aumento en 0.60561 % de los ingresos mensuales, en función del modelo planteado.
Prueba de Hipótesis:
Evaluar el modelo cuando el aumento en un año de educación genera un aumento en 10 % del salario mensual, es decir, se tendrá la hipótesis de que el valor del coeficiente para la variable A.educ es 0.1. Para esto, se debe usar el comando posregresión llamado linearHypothesis( ), del paquete car. La sintaxis será la siguiente:
## Linear hypothesis test
##
## Hypothesis:
## A.Educ = 0.1
##
## Model 1: restricted model
## Model 2: log(Salario.Mes) ~ Horas.Sem + Coef.Int + A.Educ + A.Exper +
## Edad
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 930 142.34
## 2 929 136.36 1 5.9753 40.708 2.785e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Como se puede observar, se rechaza la hipótesis nula, cuanto más se aproxime al valor del coeficiente, más cercano será a la unidad. Además, con un valor se significancia menor al 0.05, se rechaza la hipótesis nula y se acepta la hipótesis alterna.
Bibliografía
- The Comprehensive R Archive Network. Recuperado el 24 de Julio de 2018 de: https://cran.r-project.org/
- Guisande, C.; Vaamonde, A. y Barreiro, A. (2013) Tratamiento de datos con R, STATISTICA y SPSS. Vigo, España: Díaz de Santos.
- Olive, D. (2010) Multiple Linear and 1D Regression. Illinois, EEUU: Southern Illinois University Department of Mathematics.