Ejemplo de motivación.

Supongamos que nos realizan una consulta sobre cómo promover las ventas de determinado producto. Consideremos los siguientes datos que incluyen el total de ventas de ese producto (sales, en miles de unidades) en 200 diferentes mercados, junto con los gastos (en miles de dólares) realizados en publicidad en los tres siguientes medios: TV, radio y periódico (newspaper).

library(ISLR)
Datos=read.csv("images/archivos/Advertising.csv", header=TRUE, sep="," )
head(Datos)
##   X    TV radio newspaper sales
## 1 1 230.1    38        69  22.1
## 2 2  44.5    39        45  10.4
## 3 3  17.2    46        69   9.3
## 4 4 151.5    41        58  18.5
## 5 5 180.8    11        58  12.9
## 6 6   8.7    49        75   7.2
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 ...

El gasto en publicidad en televisión (TV) es una variable que se puede controlar, de manera que a partir de ésta se podría analizar si a mayor gasto en publicidad en TV hay mayores ventas. En caso afirmativo podríamos sugerir aumentar el gasto en publicidad en TV, en caso contrario podríamos sugerir realizar el gasto en otro tipo de publicidad. Veamos como se ven los datos.

require(ggplot2)
require(ggiraph)
require(ggiraphExtra)

ggPoints(aes(x=TV,y=sales),smooth=FALSE, data=Datos,interactive=TRUE)

Se observa que a mayor gasto en publicidad en TV mayor número de ventas, aunque también se observa que hay mayor variabilidad al incrementarse el gasto en publicidad.

¿Se podría considerar que la relación es lineal?

Modelo de regresión lineal simple.

En el modelo de regresión lineal simple se asume lo siguiente:

\[y_i=\alpha+\beta x_i + \epsilon_i, i=1,..., n,\]

donde \(E(\epsilon_i)=0, \; V(\epsilon_i)=\sigma^2 \;\; \text{y} \;\; Cov(\epsilon_i, \epsilon_j)=0, \; i\neq j \; \; \forall i,j = 1,...,n.\)

En general, se realiza el supuesto adicional siguiente para realizar inferencias sobre los parámetros:

\[y_i \sim N(\mu_i, \sigma^2),\] con \(\mu_i=\alpha+\beta x_i\) y \(Cov(y_i, y_j)=0, \; i\neq j \; \; \forall i,j = 1,...,n.\)

El supuesto de normalidad es el que hace a este modelo un caso particular de los modelos lineales generalizados. Pues el objetivo se reduce a realizar inferencias a partir de \(E(Y|x_1,..., x_n)\).

Estimación

La estimación de los parámetros en el modelo de regresión lineal simple se puede realizar a partir de los estimadores obtenidos por el método de mínimos cuadrados o bien por máxima verosimilitud.

En el caso de mínimos cuadrados, se buscan los valores de \(\alpha\) y \(\beta\) que minimizan la suma de cuadrados de los errores definidos como \(e_i=y_i-\alpha-\beta x_i, \; i=1,...,n\). Es decir, minimizar \(\sum_{i=1}^n (y_i-\alpha-\beta x_i)^2\).

En el caso de máxima verosímilitud, se debe maximizar la log-verosimilitud:

\[l = log(L) =-\frac{n}{2}log(2\pi) - \frac{n}{2} log(\sigma^2)-\frac{1}{2 \sigma^2}\sum_{i=1}^n (y_i-\alpha-\beta x_i)^2. \] Los estimadores de \(\alpha\) y \(\beta\) obtenidos con ambos métodos coinciden y son:

\[\begin{align} \widehat{\alpha} = \overline{y} - \widehat{\beta}\overline{x} \end{align}\]

\[\begin{align} \widehat{\beta} = \frac{\sum_{i=1}^{n}(x_{i}-\overline{x})(y_{i}-\overline{y})}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} \end{align}\]

Nota: para encontrar los estimadores es necesario derivar con respecto a cada uno de los parámetros y resolver el sistema de ecuaciones que se obtiene al igualar las derivadas igual a cero

Ambos estimadores son insesgados, es decir, \(E(\widehat{\alpha})=\alpha\) y \(E(\widehat{\beta})=\beta\).

Para obtener las propiedades de los estimadores basta escribirlos como combinación lineal de los valores de \(y_i\)

xbar=mean(Datos$TV)
SSx=sum((Datos$TV-xbar)^2)
ybar=mean(Datos$sales)
SSy=sum((Datos$sales-ybar)^2)
SSxy=sum((Datos$sales-ybar)*(Datos$TV-xbar))
beta=SSxy/SSx
alpha=ybar-beta*xbar
options(digits=5)
print (c(alpha, beta))
## [1] 7.032594 0.047537

La linea azul corresponde a la recta ajustada

ggPoints(aes(x=TV,y=sales),method="lm", data=Datos,interactive=TRUE)

Las varianzas de los estimadores de \(\alpha\) y \(\beta\) son:

\[\begin{align} Var(\hat{\beta}) = \frac{\sigma^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} \end{align}\]

\[\begin{align} Var(\hat{\alpha}) = \sigma^2 [\frac{\sum_{i=1}^{n} x_i^2}{n\sum_{i=1}^{n}(x_{i}-\overline{x})^2}] \end{align}\]

Un estimador insesgado de \(\sigma^2\) es

\[\begin{align} \widehat{\sigma}^2= \dfrac{\sum_{i=1}^n (y_i-\widehat{\alpha}-\widehat{\beta} x_i)^2}{n-2} \end{align}\]

n=length(Datos$TV)
SSE=sum((Datos$sales-alpha-beta*Datos$TV)^2)/(n-2)
print(SSE)
## [1] 10.619

Usando esta estimación también se pueden obtener estimaciones de las varianzas de los parámetros \(\alpha\) y \(\beta\).

A partir de esto es posible observar que:

\[\dfrac{\hat{\beta} -\beta}{\sqrt{ \widehat{Var}(\hat{\beta})}} \sim t_{n-2}\]

Intervalos de confianza sobre \(\alpha\) y \(\beta\).

De donde es posible obtener un intervalo al \((1-\delta)\)% de confianza para \(\beta\).

\[[\hat{\beta} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\beta})}, \; \hat{\beta} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\beta})}] \]

De manera similar, un intervalo de confianza para \(\alpha\) es:

\[[\hat{\alpha} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\alpha})}, \; \hat{\alpha} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\alpha})}] \]

var.alpha=SSE*sum(Datos$TV^2)/(n*SSx)
var.beta=SSE/SSx

print(sqrt(var.alpha))
## [1] 0.45784
print(sqrt(var.beta))
## [1] 0.0026906
#Intervalos al 95%
delta=.05
#alpha
print(c(alpha-qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.alpha),
        alpha+qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.alpha)))
## [1] 6.1297 7.9355
#beta
print(c(beta-qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.beta),
        beta+qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.beta)))
## [1] 0.042231 0.052843

También es posible realizar pruebas de hipótesis sobre combinaciones lineales de los parámetros.

Sea \(b_0\) una constante.

Pruebas de hipótesis sobre \(\alpha\) y \(\beta\).

a)

\[H_0: \beta = b_0 \] vs \[H_a: \beta \neq b_0 \]

b)

\[H_0: \beta \leq b_0 \] vs \[H_a: \beta > b_0 \]

En ambos casos, la estadística de prueba es

\[t^*= \dfrac{\hat{\beta} -b_0}{\sqrt{ \widehat{Var}(\hat{\beta})}}\]

La regla de decisión para una significancia \(\delta\) en el caso a) es:

Se rechaza \(H_0\) si \(|t^*|>t_{n-2, 1-\frac{\delta}{2}}\)

Para el caso b) es:

Se rechaza \(H_0\) si \(t^*>t_{n-2, 1-\delta}\)

**Nota. Cuando el objetivo es estudiar si hay una relación lineal entre la variable dependiente y la independiente se considera \(b_0=0\) con el caso a), mientras que si, por ejemplo, se quiere estudiar si la pendiente es positiva se considera \(b_0=0\) con el caso b).

Las pruebas de hipótesis asociadas al parámetro \(\alpha\) tienen las mimas reglas de decisión, sólo que en ese caso

\[t^*= \dfrac{\hat{\alpha} -a_0}{\sqrt{ \widehat{Var}(\hat{\alpha})}}\]

¿En el ejemplo se puede concluir que por cada mil dólares gastados en publicidad el aumento de ventas supera las 1000 unidades?

Escribir las hipótesis asociadas, estadística de prueba y regla de decisión. Usar \(\delta=.05\)

.

.

.

test=(beta-1)/sqrt(var.beta)
delta=.05
t= qt(1-delta, n-2,lower.tail = TRUE)
print(test)
## [1] -354
print(t)
## [1] 1.6526

Con los resultados anteriores concluir.

.

.

.

Usando la función lm() de R.

fit1=lm(sales~TV, data=Datos)
print(summary(fit1))
## 
## Call:
## lm(formula = sales ~ TV, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.386 -1.955 -0.191  2.067  7.212 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.03259    0.45784    15.4   <2e-16 ***
## TV           0.04754    0.00269    17.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.26 on 198 degrees of freedom
## Multiple R-squared:  0.612,  Adjusted R-squared:  0.61 
## F-statistic:  312 on 1 and 198 DF,  p-value: <2e-16
confint(fit1, level=.95)
##                2.5 %   97.5 %
## (Intercept) 6.129719 7.935468
## TV          0.042231 0.052843

Identificar lo que se calcula en cada salida de R. Considerar la prueba de hipótesis tipo a) con \(a_0=0\) y \(b_0=0\).