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.
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.\)
En ocasiones nos interesa proporcionar intervalos de confianza o de predicción de forma simultánea. Por ejemplo, se tiene interés en saber cuáles serían las ventas en dos mercados A y B con gastos en publicidad \(x_{h1}\) y \(x_{h2}\). Es decir, para cada mercado nos interesa saber sus ventas de forma simultánea y no sólo el total de ventas en ambos mercados.
Todos los intervalos anteriores se desarrollaron para un parámetro en particular \(\theta\) o bien una función unidimensional de nuevas observaciones \(y_{h1},..., y_{hk}\).
Existen tres métodos comúnmente usados para calcular intervalos de confianza/predicción simultáneos.
Método de Bonferroni Supongamos que se tiene interés en k predicciones o bien k intervalos de confianza, entonces se calculan los intervalos tal como se han visto pero en lugar de usar el valor de \(t_{n-2, 1-\frac{\delta}{2}}\) se usan los de \(t_{n-2, 1-\frac{\delta}{2k}}\).
Método de Scheffé
Se calculan los intervalos tal como se han visto pero en lugar de usar el valor de \(t_{n-2, 1-\frac{\delta}{2}}\) se usan los de \((2 {F}_{2 , n-2 , 1-\delta})^{1/2}\).
Se calculan los intervalos tal como se han visto pero en lugar de usar el valor de \(t_{n-2, 1-\frac{\delta}{2}}\) se usan los de \((k {F}_{k , n-2 , 1-\delta})^{1/2}\).
Método del máximo módulo
Se calculan los intervalos tal como se han visto pero en lugar de usar el valor de \(t_{n-2, 1-\frac{\delta}{2}}\) se usan los de \(u_{k,n-2, 1-\delta}\).
fit1=lm(sales~TV, data=Datos)
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.3 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
Por ejemplo. Supongamos que deseamos saber las ventas para tres mercados en donde se realizarán los siguientes gastos en publicidad:
(xh1=min(Datos$TV))
## [1] 0.7
(xh2=mean(Datos$TV))
## [1] 147
(xh3=max(Datos$TV))
## [1] 296
Estamos hablando de predicción. Entonces, la función predict() nos podría dar los intervalos de Bonferroni, pero deberíamos ajustar el nivel.
newdata <- data.frame(TV = c(xh1, xh2, xh3))
delta=.05
level1= 1-delta/length(newdata$TV)
pred.w.clim <- predict(fit1, newdata, interval = "prediction", level = level1)
options(digits=10)
head(pred.w.clim)
## fit lwr upr
## 1 7.065869197 -0.8788634712 15.01060187
## 2 14.022500000 6.1348555990 21.91014440
## 3 21.122453773 13.1753534593 29.06955409
Se pueden hacer los cálculos a mano. Por ejemplo, uno de esos intervalos de predicción.
delta=.05
ttab=qt(1-delta/(2*length(newdata$TV)), n-2,lower.tail = TRUE)
yfith=alpha+beta*xh1
Var.yfithpred=SSE*(1+1/n+(xh1-xbar)^2/SSx)
print(c(yfith-ttab*sqrt(Var.yfithpred),
yfith+ttab*sqrt(Var.yfithpred)))
## [1] -0.8788634712 15.0106018660
El paquete investr también se puede usar para intervalos de predicción/confianza simultáneos.
library(investr)
pred.w.clim2=predFit(fit1, newdata, interval = "prediction", level = 0.95, adjust = c("Bonferroni"), k=length(newdata$TV))
head(pred.w.clim2)
## fit lwr upr
## 1 7.065869197 -0.8788634712 15.01060187
## 2 14.022500000 6.1348555990 21.91014440
## 3 21.122453773 13.1753534593 29.06955409
También se pueden obtener los intervalos por el método de Scheffé, por ejemplo.
library(investr)
pred.w.climScheffe=predFit(fit1, newdata, interval = "prediction", level = 0.95, adjust = c("Scheffe"), k=length(newdata$TV))
head(pred.w.climScheffe)
## fit lwr upr
## 1 7.065869197 -2.212128932 16.34386733
## 2 14.022500000 4.811170550 23.23382945
## 3 21.122453773 11.841690666 30.40321688
A mano, uno de los intervalos de Scheffé
delta=.05
Ftab=qf(1-delta,length(newdata$TV) ,n-2,lower.tail = TRUE)
yfith=alpha+beta*xh1
Var.yfithpred=SSE*(1+1/n+(xh1-xbar)^2/SSx)
print(c(yfith-sqrt(length(newdata$TV)*Ftab)*sqrt(Var.yfithpred),
yfith+sqrt(length(newdata$TV)*Ftab)*sqrt(Var.yfithpred)))
## [1] -2.212128932 16.343867327
En general, la banda de intervalos simultáneos llamada Working-Hotelling band se puede graficar como sigue (sólo es de confianza!).
plotFit(fit1, interval = "confidence", adjust = "Scheffe")
Comparación con los intervalos de confianza individuales:
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)
pred.w.climScheffe=predFit(fit1, newdata, interval = "confidence", level = 0.95, adjust = c("Scheffe"))
matplot(newdata$TV, cbind(pred.w.clim, pred.w.climScheffe[,2:3]),
lty = c(1,2,2, 1, 1), type = "l", ylab = "Nuevas ventas dado un gasto de TV", xlab = "TV")
points(Datos$TV, Datos$sales, cex=.5)
Nota: no se pueden obtener bandas de predicción sobre \(x\), pues el método de Scheffé y de Bonferroni dependen de \(k\), a diferencia de bandas de confianza donde el método de Scheffé no depende de \(k\).
Nota 2. Se pueden obtener intervalos de confianza de \(k\) combinaciones lineales usando el paquete multcomp por el método del máximo módulo; los de Bonferroni y Scheffé se pueden calcular dando el valor del cuantil asociado.
library(multcomp)
newdata <- data.frame(TV = c(xh1, xh2, xh3))
#Para alpha
K=matrix(c(1,1,1,xh1, xh2, xh3), ncol=2, nrow=3)
m2=0
prueba1=glht(fit1, linfct=K, rhs=m2)
#Máximo modulo
confint(prueba1)
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = sales ~ TV, data = Datos)
##
## Quantile = 2.3618755
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 7.0658692 5.9883428 8.1433956
## 2 == 0 14.0225000 13.4782724 14.5667276
## 3 == 0 21.1224538 20.0283481 22.2165595
#Bonferroni
confint(prueba1, calpha=qt(1-.05/(2*length(newdata$TV)), n-2,lower.tail = TRUE) )
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = sales ~ TV, data = Datos)
##
## Quantile = 2.4144918
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 7.0658692 5.9643384 8.1674000
## 2 == 0 14.0225000 13.4661485 14.5788515
## 3 == 0 21.1224538 20.0039743 22.2409332
#Scheffé
confint(prueba1, calpha=sqrt(2*qf(1-delta,2 ,n-2,lower.tail = TRUE) ) )
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = sales ~ TV, data = Datos)
##
## Quantile = 2.4663813
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 7.0658692 5.9406656 8.1910728
## 2 == 0 14.0225000 13.4541920 14.5908080
## 3 == 0 21.1224538 19.9799373 22.2649703
library(investr)
pred.w.climBon=predFit(fit1, newdata, interval = "confidence", level = 0.95, adjust = c("Bonferroni"),k=length(newdata$TV) )
head(pred.w.climBon)
## fit lwr upr
## 1 7.065869197 5.964338438 8.167399957
## 2 14.022500000 13.466148460 14.578851540
## 3 21.122453773 20.003974340 22.240933207
pred.w.climScheffe=predFit(fit1, newdata, interval = "confidence", level = 0.95, adjust = c("Scheffe") )
head(pred.w.climScheffe)
## fit lwr upr
## 1 7.065869197 5.940665596 8.191072799
## 2 14.022500000 13.454191988 14.590808012
## 3 21.122453773 19.979937257 22.264970290