<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: { 
      equationNumbers: { 
            autoNumber: "AMS",
            formatNumber: function (n) {return '1.'+n}
      } 
  }
});
</script>

---
title: "Regresión lineal simple (1/4)"
author: Gonzalo Pérez
date: 4 de febrero de 2020
header-includes:
  - \usepackage{amsmath}
output: 
    html_document: 
    df_print: paged
    fig_caption: yes
    highlight: espresso
    incremental: yes
    keep_md: yes
    transition: slower
classoption: xcolor=table

---

```{r setup, include=FALSE}
rm(list = ls(all.names = TRUE))
gc()
library(reticulate)
knitr::opts_chunk$set(echo = F, warning = F, message = F, error = F, fig.height = 4, fig.width = 8)
library(xtable)
library(knitr)
library(tidyverse)
library(latex2exp)
options(digits=2)
set.seed(20202)
```

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

```{r, include=TRUE,  echo=TRUE}
library(ISLR)
Datos=read.csv("images/archivos/Advertising.csv", header=TRUE, sep="," )
head(Datos)
str(Datos)

```

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. 

- $Y$- Total de ventas del producto (sales).
- $X$- Total de gasto en publicidad en TV (TV).


```{r, echo=TRUE}
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$.

```{r, echo=FALSE}
fit1=lm(sales~TV, data=Datos)
plot(sales~TV, data=Datos)
abline(fit1, col="red")
attach(Datos)
segments(TV, sales, TV, predict(fit1))
detach(Datos)
```


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$**

```{r, echo=TRUE}
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))
```

La linea azul corresponde a la recta ajustada

```{r, echo=TRUE}
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}$$


```{r, echo=TRUE}
n=length(Datos$TV)
SSE=sum((Datos$sales-alpha-beta*Datos$TV)^2)/(n-2)
print(SSE)
```


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})}]  $$

```{r, echo=TRUE}

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

print(sqrt(var.alpha))
print(sqrt(var.beta))


#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)))

#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)))
```


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 en TV 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$**

.

.

.




```{r, echo=TRUE}

test=(beta-1)/sqrt(var.beta)
delta=.05
t= qt(1-delta, n-2,lower.tail = TRUE)
print(test)
print(t)
```

**Con los resultados anteriores concluir.**

.

.

.


## Usando la función lm() de R. 



```{r, echo=TRUE}
fit1=lm(sales~TV, data=Datos)
print(summary(fit1))
confint(fit1, level=.95)
```


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

