--- title: "Regresión lineal simple (2/4)" author: Gonzalo Pérez date: 10 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=FALSE} require(ggplot2) require(ggiraph) require(ggiraphExtra) ggPoints(aes(x=TV,y=sales),method="lm", data=Datos,interactive=TRUE) ``` ## 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_i$ y $y_j$ variables independientes para $\; i\neq j \; \; \forall i,j = 1,...,n.$ ## Estimaciones para combinaciones lineales de los parámetros. En la clase pasada observamos que los estimadores de $\alpha$ y $\beta$ 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} ```{r, echo=FALSE} 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 ``` ```{r, echo=FALSE} n=length(Datos$TV) SSE=sum((Datos$sales-alpha-beta*Datos$TV)^2)/(n-2) ``` ```{r, echo=FALSE} var.alpha=SSE*sum(Datos$TV^2)/(n*SSx) var.beta=SSE/SSx ``` ```{r, echo=FALSE} fit1=lm(sales~TV, data=Datos) ``` En general hay dos usos de la regresión: - Describir un fenómeno (en general, en promedio). - Predecir (nuevas observaciones). Hay más incertidumbre cuando se predice que cuando se describe un fenómeno en general. ### Descripción del fenómeno En este caso nos interesa realizar inferencias sobre aspectos relacionados con combinaciones lineales de los parámetros. Por ejemplo: - Inferencias sobre los parámetros. Aquí la idea fue describir, por ejemplo, el cambio que en general se espera en la variable $Y$ cuando se aumenta en una unidad la variable $X$. - Otro aspecto a considerar es $E(Y|X=x)=\alpha+\beta x$. Es decir, cual sería el **valor esperado** o el **promedio** que se espera cuando el valor de la variable independiente es $X=x$. En general, supongamos que se tiene interés en estimar $$\theta=z_0 \alpha + z_1 \beta,$$ donde $z_0$ y $z_1$ son constantes, no ambas iguales a cero. Consideremos $\widehat{\theta}=z_0 \widehat{\alpha} + z_1 \widehat{\beta}$. Las propiedades siguientes se cumplen: - De acuerdo con el Teorema de Gauss-Markov, el estimador $z_0 \widehat{\alpha} + z_1 \widehat{\beta}$ es el de mínima varianza entre los estimadores lineales e insesgados. - $E(\widehat{\theta})=E(z_0 \widehat{\alpha} + z_1 \widehat{\beta})=\theta=z_0 \alpha + z_1 \beta$ - $V(\widehat{\theta})=V(z_0 \widehat{\alpha} + z_1 \widehat{\beta})=\sigma^2[\dfrac{z_0^2}{n}+\dfrac{(z_1-z_0\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]$ - $z_0 \widehat{\alpha} + z_1 \widehat{\beta} \sim N(z_0 \alpha + z_1 \beta, \sigma^2[\dfrac{z_0^2}{n}+\dfrac{(z_1-z_0\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} )$ - Un estimador insesgado de $V(\widehat{\theta})$ es $\hat{V}(\widehat{\theta})= \widehat{Var}(z_0 \widehat{\alpha} + z_1 \widehat{\beta})=\hat{\sigma}^2[\dfrac{z_0^2}{n}+\dfrac{(z_1-z_0\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]$ - Un intervalo de confianza para $\theta$ es $$[\hat{\theta} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\theta})}, \; \hat{\theta} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\theta})}] $$ - Considerando la siguiente prueba de hipótesis $$H_0: z_0 \alpha + z_1 \beta = c_0 $$ vs $$H_a: z_0 \alpha + z_1 \beta \neq c_0 $$ La estadística de prueba es $$t^*= \dfrac{z_0 \widehat{\alpha} + z_1 \widehat{\beta} -c_0}{\sqrt{ \widehat{Var}(z_0 \widehat{\alpha} + z_1 \widehat{\beta})}}$$ La regla de decisión para una significancia $\delta$ es: Rechazar $H_0$ si $|t^*|> t_{n-2, 1-\frac{\delta}{2}}$. Casos particulares son: - $E(Y|X=x)=\alpha+\beta x$, donde $z_0=1$ y $z_1=x$. - $E(Y|X=x+z_1)-E(Y|X=x)= \alpha+(x+z_1)\beta-(\alpha+x\beta)=z_1\beta$, donde $z_0=0$. - $E(Y|X=0)=\alpha$, donde $z_0=1$ y $z_1=0$. En R, se pueden obtener los intervalos de confianza para $E(Y|X=x)$ de forma muy sencilla usando la función predict(). ```{r, echo=TRUE} fit1=lm(sales~TV, data=Datos) newdata <- data.frame(TV = seq(min(Datos$TV), max(Datos$TV), 5)) pred.w.clim <- predict(fit1, newdata, interval = "confidence", level = 0.95) options(digits=10) head(pred.w.clim) matplot(newdata$TV, cbind(pred.w.clim), lty = c(1,2,2), type = "l", ylab = "Estimation of E(sales|TV)", xlab = "TV") points(Datos$TV, Datos$sales, cex=.5) ``` También se pueden hacer los cálculos a mano ```{r, echo=TRUE} delta=.05 ttab=qt(1-delta/2, n-2,lower.tail = TRUE) x=min(Datos$TV) yfit=alpha+beta*x Var.yfit=SSE*(1/n+(x-xbar)^2/SSx) print(c(yfit-ttab*sqrt(Var.yfit), yfit+ttab*sqrt(Var.yfit))) ``` O usar la función asociada a pruebas de hipótesis del modelo lineal general usando matrices. Esta prueba es más general y sirve para cualquier valor de $z_0$ y $z_1$. ```{r, echo=TRUE} library(multcomp) #Para alpha K=matrix(c(1,0), ncol=2, nrow=1) m=0 prueba1=glht(fit1, linfct=K, rhs=m) confint(prueba1) #Para beta K=matrix(c(0,1), ncol=2, nrow=1) m=0 prueba2=glht(fit1, linfct=K, rhs=m) confint(prueba2) #Para E(Y|X), con X el mínimo valor observado en los datos K=matrix(c(1,min(Datos$TV)), ncol=2, nrow=1) m=0 prueba3=glht(fit1, linfct=K, rhs=m) confint(prueba3) ``` Además la función glht permite hacer pruebas de hipótesis. Por ejemplo repetiremos la prueba de la sesión pasada $$H_0: \beta \leq 1 $$ vs $$H_a: \beta > 1 $$ ```{r, echo=TRUE} #Para beta K=matrix(c(0,1), ncol=2, nrow=1) m=1 summary(glht(fit1, linfct=K, rhs=m, alternative="greater")) ``` Por default se realiza la prueba $$H_0: \beta = 1 $$ vs $$H_a: \beta \neq 1 $$ Lo que se modifica es la forma de calcular el p-value ```{r, echo=TRUE} #Para beta K=matrix(c(0,1), ncol=2, nrow=1) m=1 summary(glht(fit1, linfct=K, rhs=m)) ``` Hagamos una prueba donde la variable independiente es **newspaper**. Así podremos comprobar la equivalencia en el cálculo de los p-values. ```{r, echo=TRUE} #Para beta fit1=lm(sales~newspaper, data=Datos) K=matrix(c(0,1), ncol=2, nrow=1) m=0 #Pruebas simultaneas (sólo hay una, así que debe coincidir con la global) summary(glht(fit1, linfct=K, rhs=m)) #Prueba global summary(glht(fit1, linfct=K, rhs=m),test=Ftest()) #La prueba de la tabla anova coincide en este caso con la prueba global y con la individual summary(fit1) #Resultafo de la prueba de la tabla anova. anova(fit1) ### La siguiente no es equivalente, recae en supuestos asintóticos summary(glht(fit1, linfct=K, rhs=m),test=Chisqtest()) ``` Otra forma de hacer algunas pruebas de hipótesis es usando la función anova(). Para esto es necesario ajustar dos modelos, el completo o sin restrincciones y uno restringido. Por ejemplo, para la prueba $$H_0: \beta = 1 $$ vs $$H_a: \beta \neq 1 $$ ```{r, echo=TRUE} fit1=lm(sales~TV, data=Datos) fit2=lm(sales~offset(TV), data=Datos) anova(fit1, fit2) ``` ### Predicción