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.26

En 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-16

O 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-16

Dibujamos 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.002690607

3.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.972017

Por tanto:

source("02_funciones.R")
rls_ic(m1,0.05)
## $beta0_ic
## [1] 6.129719 7.935468
## 
## $beta1_ic
## [1] 0.04223072 0.05284256

En R utilizamos la funcion confint:

confint(m2, level = 0.95)
##                  2.5 %     97.5 %
## (Intercept) 6.12971927 7.93546783
## x           0.04223072 0.05284256

3.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-43

Por 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.258656

4.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.6118751

5 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.85686

En R:

xp = data.frame(x=90)
predict(m2, newdata = xp, level = 0.95, interval="confidence")
##        fit      lwr      upr
## 1 11.31089 10.76492 11.85686
xp = data.frame(TV = 90)
predict(m3, newdata = xp, level = 0.95, interval="confidence")
##        fit      lwr      upr
## 1 11.31089 10.76492 11.85686

Predecimos 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.760170

En R:

xp = data.frame(x=90)
predict(m2, newdata = xp, level = 0.95, interval="prediction")
##        fit      lwr      upr
## 1 11.31089 4.861613 17.76017

5.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.51265

6 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-16

Igual 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