1 Modelo

Datos:

\((y_i,x_{1i},x_{2i},\cdots,x_{ni}), i = 1,\ldots,n\)

Modelo:

\(y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_k x_{ki} + u_i, \ u_i \sim N(0,\sigma^2)\)

Parámetros:

\(\beta_0, \beta_1, \beta_2, \cdots, \beta_k, \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_{1i} + \beta_2 x_{2i} + \cdots + \beta_k x_{ki}\)

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

Se define la suma de residuos al cuadrado o residual sum of squares (RSS):

\(RSS = \sum \limits _{i=1}^n e_i^2\)

donde \(e_i = y_i - (\hat \beta_0 + \hat \beta_1 x_{1i} + \hat \beta_2 x_{2i} + \cdots + \hat \beta_k x_{ki})\)

El metodo de minimos cuadrados minimiza esta cantidad.

\(\frac{\partial RSS}{\partial \hat \beta_i} = 0 \Rightarrow \hat \beta_i\)

\(\hat \sigma ^2 = \frac{\sum _{i=1}^n (y_i - \hat y)^2}{n-k-1} = \frac{RSS}{n-k-1}\)

Esos parámetros los podemos calcular con R:

m1 = lm(sales ~ TV + radio + newspaper, data = datos)
summary(m1)
## 
## Call:
## lm(formula = sales ~ TV + radio + newspaper, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8277 -0.8908  0.2418  1.1893  2.8292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.938889   0.311908   9.422   <2e-16 ***
## TV           0.045765   0.001395  32.809   <2e-16 ***
## radio        0.188530   0.008611  21.893   <2e-16 ***
## newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16
  • \(\hat \beta_0\) = 2.94. Como p-valor < \(\alpha\), es un parámetro significativo.
  • \(\hat \beta_1\) = 0.05. Si mantenemos la inversión en radio y newspaper constante, un incremento de 1000 $ en TV, por término medio supone un aumento en las ventas de 0.05*1000 = 50 unidades. Según el pvalor, es un parámetro significativo.
  • \(\hat \beta_2\) = 0.19. Si mantenemos la inversión en TV y newspaper constante, un incremento de 1000 $ en TV, por término medio supone un aumento en las ventas de 190 unidades. Según el pvalor, es un parámetro significativo.
  • \(\hat \beta_3\) = -0.001. Según el pvalor, es un parámetro NO significativo, luego no influye en las ventas. Sin embargo, si analizamos la regresión simple de newspaper:
m2 = lm(sales ~ newspaper, data = datos)
summary(m2)
## 
## Call:
## lm(formula = sales ~ newspaper, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.2272  -3.3873  -0.8392   3.5059  12.7751 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.35141    0.62142   19.88  < 2e-16 ***
## newspaper    0.05469    0.01658    3.30  0.00115 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.092 on 198 degrees of freedom
## Multiple R-squared:  0.05212,    Adjusted R-squared:  0.04733 
## F-statistic: 10.89 on 1 and 198 DF,  p-value: 0.001148

Parece un resultado contradictorio. Esto es debido a la correlación entre radio y newspaper:

cor(datos[,-1])
##                   TV      radio  newspaper     sales
## TV        1.00000000 0.05480866 0.05664787 0.7822244
## radio     0.05480866 1.00000000 0.35410375 0.5762226
## newspaper 0.05664787 0.35410375 1.00000000 0.2282990
## sales     0.78222442 0.57622257 0.22829903 1.0000000

Como vemos, la correlación entre radio y newspaper es 0.35, lo que indica que en los mercados donde se invierte en radio también se invierte en newspaper.

3 Intervalos de confianza

confint(m1, level = 0.99)
##                   0.5 %     99.5 %
## (Intercept)  2.12757072 3.75020802
## TV           0.04213632 0.04939297
## radio        0.16613095 0.21092909
## newspaper   -0.01630884 0.01423386

4 Contraste general de regresión

Existe un contraste que nos informa si existe relación entre la variable respuesta y los regresores.

\(H_0 : \beta_1 = \beta_2 = \cdots = \beta_k = 0\)

\(H_1 : \text{al menos un } \beta_j \text{ es distinto de cero}\)

Este contraste se resuelve utilizando el estadístico F. El valor de este estadístico y el pvalor se pueden consultar en la tabla obtenida conla función summary()

summary(m1)
## 
## Call:
## lm(formula = sales ~ TV + radio + newspaper, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8277 -0.8908  0.2418  1.1893  2.8292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.938889   0.311908   9.422   <2e-16 ***
## TV           0.045765   0.001395  32.809   <2e-16 ***
## radio        0.188530   0.008611  21.893   <2e-16 ***
## newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

Como pvalor < \(\alpha\), se rechaza la hipótesis nula, luego hay algún regresor distinto de cero, luego hay relación entre sales y los regresores.

5 Precisión del modelo

Residual standard error RSE:

m1_sm = summary(m1)
m1_sm$sigma
## [1] 1.68551

Como vemos, disminuye con respecto a la regresión simple.

R cuadrado:

m1_sm$r.squared
## [1] 0.8972106

R cuadrado ajustado (mejor para regresión simple):

m1_sm$adj.r.squared
## [1] 0.8956373

Este valor también ha mejorado considerablemente.

6 Predicción

Valor medio predicho - intervalo de confianza del valor medio predicho:

xp = data.frame(TV = 50, radio = 40, newspaper = 60)
predict(m1, newdata = xp, level = 0.95, interval="confidence")
##        fit      lwr      upr
## 1 12.70607 12.18819 13.22396

Valor medio predicho - intervalo de predicción:

predict(m1, newdata = xp, level = 0.95, interval="prediction")
##        fit      lwr      upr
## 1 12.70607 9.341907 16.07024

7 Hipótesis del modelo

par(mfrow=c(2,2))
plot(m1)

Como vemos, no se cumple linealidad. Una solución simple consiste en usar transformaciones no-lineales de las X. Las más comunes son: \(\log(X), \sqrt{X}, X^2\).

Podemos comprobar que ninguna de ellas corrige la linealidad:

m3 = lm(sales ~ TV + radio + newspaper + I(TV^2), data = datos)
par(mfrow=c(2,2))
plot(m3)

8 Extensiones del modelo lineal

8.1 Términos de interacción

m4 = lm(sales ~ TV * radio + newspaper, data = datos)
summary(m4)
## 
## Call:
## lm(formula = sales ~ TV * radio + newspaper, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2929 -0.3983  0.1811  0.5957  1.5009 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.728e+00  2.533e-01  26.561  < 2e-16 ***
## TV          1.907e-02  1.509e-03  12.633  < 2e-16 ***
## radio       2.799e-02  9.141e-03   3.062  0.00251 ** 
## newspaper   1.444e-03  3.295e-03   0.438  0.66169    
## TV:radio    1.087e-03  5.256e-05  20.686  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9455 on 195 degrees of freedom
## Multiple R-squared:  0.9678, Adjusted R-squared:  0.9672 
## F-statistic:  1466 on 4 and 195 DF,  p-value: < 2.2e-16

El modelo mejora considerablemente el R\(^2\). Luego invertir dinero en radio también mejora la inversión en TV.

\(sales = \beta_0 + \beta_1*TV + \beta_2*radio + \beta_3*newspaper + \beta_4*TV*radio + u\)

\(sales = \beta_0 + (\beta_1 + \beta_4*radio)*TV + \beta_2*radio + \beta_3*newspaper + u\)

La linealidad también ha mejorado:

par(mfrow=c(2,2))
plot(m4)

8.2 Términos no lineales

m5 = lm(sales ~ TV * radio + newspaper + I(TV^2), data = datos)
summary(m5)
## 
## Call:
## lm(formula = sales ~ TV * radio + newspaper + I(TV^2), data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8886 -0.3043  0.0003  0.3736  1.2637 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.079e+00  1.958e-01  25.943  < 2e-16 ***
## TV           5.100e-02  2.226e-03  22.916  < 2e-16 ***
## radio        3.321e-02  6.020e-03   5.518 1.09e-07 ***
## newspaper    3.294e-03  2.170e-03   1.518    0.131    
## I(TV^2)     -1.103e-04  6.880e-06 -16.031  < 2e-16 ***
## TV:radio     1.078e-03  3.457e-05  31.199  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6217 on 194 degrees of freedom
## Multiple R-squared:  0.9862, Adjusted R-squared:  0.9858 
## F-statistic:  2764 on 5 and 194 DF,  p-value: < 2.2e-16

Este es el mejor modelo. Comprobamos los residuos:

par(mfrow=c(2,2))
plot(m5)