Regresión Lineal Simple
Javier Cara
Curso 2018-19
1 Modelo
Datos:
\((x_i,y_i), i = 1,\ldots,n\)
Modelo:
\(y_i = \beta_0 + \beta_1 x_i + u_i, \ u_i \sim N(0,\sigma^2)\)
Parámetros:
\(\beta_0, \beta_1, \sigma^2\)
Hipótesis:
- Normalidad: \(u_i \sim Normal \Rightarrow y_i \sim Normal\)
- Varianza de los datos constante (homocedasticidad): \(Var(y_i) = \sigma^2\)
- Linealidad: \(E(y_i) = \beta_0 + \beta_1 x_i\)
2 Estimación del modelo
Vamos a utilizar un método que se denomina mínimos cuadrados.
datos = read.csv("Advertising.csv")
str(datos)## 'data.frame': 200 obs. of 5 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ TV : num 230.1 44.5 17.2 151.5 180.8 ...
## $ radio : num 37.8 39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 ...
## $ newspaper: num 69.2 45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 ...
## $ sales : num 22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...plot(datos$TV, datos$sales, xlab = "TV", ylab = "Sales", col = "red", pch = 20)
Se define la suma de residuos al cuadrado o residual sum of squares (RSS):
\(RSS = (y_1 - \hat \beta_0 - \hat \beta_1 x_1)^2 + (y_2 - \hat \beta_0 - \hat \beta_1 x_2)^2 + \cdots + (y_p - \hat \beta_0 - \hat \beta_1 x_p)^2 = \sum \limits _{i=1}^n e_i^2\)
El metodo de minimos cuadrados minimiza esta cantidad.
\(\frac{\partial RSS}{\partial \hat \beta_1} = 0 \Rightarrow \hat \beta_1 = \frac{\sum _{i=1}^n (x_i - \bar x)(y_i - \bar y)}{\sum _{i=1}^n (x_i - \bar x)^2} = \frac{cov(x,y)}{s_x^2}\)
\(\frac{\partial RSS}{\partial \hat \beta_0} = 0 \Rightarrow \hat \beta_0 = \bar y - \beta_1 \bar x\)
\(\hat \sigma ^2 = \frac{\sum _{i=1}^n (y_i - \hat y)^2}{n-2} = \frac{RSS}{n-2}\)
Esos parámetros los podemos calcular creando una función:
source("02_funciones.R")
x = datos$TV
y = datos$sales
m1 = rls(y,x)
str(m1)## List of 5
## $ x : num [1:200] 230.1 44.5 17.2 151.5 180.8 ...
## $ y : num [1:200] 22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...
## $ beta0_e: num 7.03
## $ beta1_e: num 0.0475
## $ sigma_e: num 3.26En R se hace:
m2 = lm(y ~ x)
summary(m2)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3860 -1.9545 -0.1913 2.0671 7.2124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.032594 0.457843 15.36 <2e-16 ***
## x 0.047537 0.002691 17.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 198 degrees of freedom
## Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
## F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16O también:
m3 = lm(sales ~ TV, data = datos)
summary(m3)##
## Call:
## lm(formula = sales ~ TV, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3860 -1.9545 -0.1913 2.0671 7.2124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.032594 0.457843 15.36 <2e-16 ***
## TV 0.047537 0.002691 17.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 198 degrees of freedom
## Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
## F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16Dibujamos la solución:
plot(datos$TV, datos$sales, xlab = "TV", ylab = "Sales", col = "red", pch = 20)
abline(m2, col="blue", lwd = 1)
3 Precisión de los coeficientes estimados
Los estimadores de los coeficientes son variables aleatorias. Se puede demostrar que:
\(\hat \beta_0 \sim N \left( \beta_0, SE(\hat \beta_0) \right)\)
\(\hat \beta_1 \sim N \left( \beta_1, SE(\hat \beta_1) \right)\)
Donde SE significa Standard Error y mide la precisión de los estimadores.
\(SE(\hat \beta_0) = \frac{\sigma^2}{n}\left(1 + \frac{\bar x^2}{s_x^2} \right)\)
\(SE(\hat \beta_1) = \frac{\sigma^2}{ns_x^2}\)
source("02_funciones.R")
rls_se(m1)## $beta0_se
## [1] 0.4578429
##
## $beta1_se
## [1] 0.0026906073.1 Intervalos de confianza
Se calcula con la formula:
\(\hat \beta_1 \pm t_{n-2;\alpha/2} \cdot SE(\hat \beta_1)\),
donde \(t_{n-2;\alpha/2}\) es el valor de una t-student con n-2 grados de libertad que deja una probabilidad de \(\alpha/2\) a su derecha.
n = length(y)
curve(dt(x,df=n-2), from = -3, to = 3, xlab = "x", ylab = "t-student con 298 grados de libertad")
\(t_{198;0.05/2}\) =
qt(0.975,df=n-2)## [1] 1.972017Por tanto:
source("02_funciones.R")
rls_ic(m1,0.05)## $beta0_ic
## [1] 6.129719 7.935468
##
## $beta1_ic
## [1] 0.04223072 0.05284256En R utilizamos la funcion confint:
confint(m2, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 6.12971927 7.93546783
## x 0.04223072 0.052842563.2 Contraste de hipotesis
También estamos interesados en resolver preguntas del tipo: \(\beta_1 = 0\)?
\(H_0 : \beta_1 = 0\)
\(H_1 : \beta_1 \neq 0\)
Para ello, utilizamos que \(\hat \beta_1\) es una variable aleatoria.
\(\hat \beta_1 \sim N \left( \beta_1, SE(\hat \beta_1) \right)\)
Suponiendo que \(H_0\) es cierta, se cumple:
\(t_0 = \frac{\hat \beta_1}{SE(\hat \beta_1)} \sim t_{n-2}\)
Si el valor obtenido es muy probable, es que la hipótesis que hemos adoptado es plausible. Esa probabilidad se denomina p-valor.
m1se = rls_se(m1)
t0 = m1$beta1_e/m1se$beta1_se
( pvalor = pt(t0, df = n-2, lower.tail = F) )## [1] 7.336949e-43Por tanto, es muy poco probable que ese valor corresponda a una t con 198 grados de libertad, lo que indica que la hipótesis que hemos adoptado no es cierta, luego \(\beta_1 \neq 0\).
4 Precisión del modelo
Habitualmente se utilizan dos varialbes:
4.1 Residual Standard Error (RSE)
Se define como
\(RSE = \sqrt{ \frac{\sum _{i=1}^n (y_i - \hat y)^2}{n-2} } = \sqrt{ \frac{RSS}{n-2} }\)
Si recordamos, RSE es una estimación de \(\sigma\). Como
\(y_i = \beta_0 + \beta_1 x_i + u_i, \ u_i \sim N(0,\sigma^2)\)
cuanto mayor sea \(\sigma\), mayor es la variabilidad de los datos con respecto a la recta.
# RSE
m3_sm = summary(m3) # consultar str(m3)
m3_sm$sigma ## [1] 3.2586564.2 \(R^2\)
El problema de RSE es que depende de las unidades de los datos.
Se puede demostrar que
\(\sum _{i=1}^n (y_i - \bar y)^2 = \sum _{i=1}^n (\hat y_i - \bar y)^2 + \sum _{i=1}^n (y_i - \hat y_i)^2\)
VT = VE + VNE (se observa que VNE = RSE)
Se define:
\(R^2 = \frac{VE}{VT} = \frac{VT-VNE}{VT} = 1-\frac{VNE}{VT}\)
- \(R^2 \approx 0\). Los datos se alejan de la recta.
- \(R^2 \approx 1\). Los datos se acercan a la recta.
m3_sm$r.squared## [1] 0.61187515 Predicción
5.1 Cálculo de la predicción
Para un valor \(x_p\), el valor medio predicho por el modelo es (Predicción puntual):
\(\hat m_p = \hat \beta_0 + \hat \beta_1 x_p\)
Se puede demostrar que:
\(\hat m_p \sim N \left( m_p, \sigma^2v_{pp} \right)\)
\(m_p = \beta_0 + \beta_1 x_p\)
\(v_{pp} = \frac{1}{n}\left(1 + \frac{(x_p - \bar x)^2}{s_x^2} \right)\)
Por tanto, el Intervalo de Confianza es:
\(m_p \in \hat m_p \pm t_{n-2;\alpha/2}*RSE*\sqrt{v_{pp}}\)
source("02_funciones.R")
yp = rls_pred(m1,90,0.05)
c(yp$mpp,yp$mpp1,yp$mpp2)## [1] 11.31089 10.76492 11.85686En R:
xp = data.frame(x=90)
predict(m2, newdata = xp, level = 0.95, interval="confidence")## fit lwr upr
## 1 11.31089 10.76492 11.85686xp = data.frame(TV = 90)
predict(m3, newdata = xp, level = 0.95, interval="confidence")## fit lwr upr
## 1 11.31089 10.76492 11.85686Predecimos todos los datos
pred1 = rls_pred(m1,x,0.05)
plot(x,y)
abline(a=m1$beta0_e, b=m1$beta1_e)
points(x,pred1$mpp, col = "red")
También se puede calcular un intervalo para la predicción de nuevos valores de la variable y, lo que se conoce como intervalo de predicción:
\(y_p = \beta_0 + \beta_1 x_p + u_p = m_p + u_p \Rightarrow\)
\(Var(y_p) = Var(m_p) + Var(u_p) = \sigma^2 + \sigma^2v_{pp} = \sigma^2(1 + v_{pp})\)
Por tanto, el Intervalo de Predicción es:
\(\hat y_p \in \hat m_p \pm t_{n-2;\alpha/2}*RSE*\sqrt{1+v_{pp}}\)
c(yp$mpp,yp$ypp1,yp$ypp2)## [1] 11.310891 4.861613 17.760170En R:
xp = data.frame(x=90)
predict(m2, newdata = xp, level = 0.95, interval="prediction")## fit lwr upr
## 1 11.31089 4.861613 17.760175.2 Precisión de la estimación
Se calcula el mean squared error:
\(MSE = \frac{1}{n} \sum \limits _{i=1}^n (y_i - \hat y_i)^2\)
Esta es una fórmula general, aplicable a todos los modelos. En el caso de regresión simple MSE es similar a RSE\(^2\).
yp = predict(m3, newdata = datos)
n = nrow(datos)
(MSE = 1/n*sum((y-yp)^2))## [1] 10.512656 Diagnosis
El modelo de regresión simple se basa en las siguientes hipótesis:
- Normalidad: \(u \sim Normal \Rightarrow y_i \sim Normal\)
- Varianza de los datos constante (homocedasticidad): \(Var(y_i) = \sigma^2\)
- Linealidad: \(E(y_i) = \beta_0 + \beta_1 x_i\)
Estas hipótesis se comprueban analizando los residuos:
\(e_i = y_i - (\beta_0 + \beta_1 x_i)\)
Los residuos en R se calculan con la función residuals. Por ejemplo, el gráfico valores predichos - residuos:
plot(yp,residuals(m3), xlab = "Fitted values", ylab = "Residuals")
Este gráfico, junto con otros adicionales, se puede obtener con:
par(mfrow=c(2,2))
plot(m2)
- Normalidad: se comprueba con el gráfico superior derecho. Los residuos se deben ajustar a la recta punteada.
- Linealidad: se comprueba con el gráfico superior izquierdo. Los residuos deben distribuirse simétricamente alrededor del eje x. Es equivalente a que la línea roja del gráfico no tenga curvatura.
- Homocedasticidad: se comprueba con el gráfico superior izquierdo. Los residuos deben tener ancho constante, no pueden aumentar con los valores predichos.
En este caso no se cumple homocedasticidad. Lo que se hace es utilizar una transformación cóncava de la variable y (generalmente log(y) \(\sqrt{y}\)) :
m4 = lm(sqrt(sales) ~ TV, data = datos)
par(mfrow=c(2,2))
plot(m4)
Parece que la raiz cuadrada mejora la hipótesis de homocedasticidad. En cuanto a la bondad del nuevo modelo, apenas cambia
summary(m4)##
## Call:
## lm(formula = sqrt(sales) ~ TV, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44955 -0.25244 0.01172 0.31823 0.85003
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7098601 0.0605679 44.74 <2e-16 ***
## TV 0.0065781 0.0003559 18.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4311 on 198 degrees of freedom
## Multiple R-squared: 0.633, Adjusted R-squared: 0.6312
## F-statistic: 341.5 on 1 and 198 DF, p-value: < 2.2e-16Igual que el MSE de las predicciones.
yp_raiz = predict(m4, newdata = datos)
n = nrow(datos)
(MSE = 1/n*sum((y-(yp_raiz)^2)^2))## [1] 10.87166