Regresión Lineal Múltiple
Dr. Víctor Chung
4/13/2022
Ejemplo 1
Vamos a establecer un modelo económico para una cadena de hamburgueserías que llamamos Big Andy’s Burger Barn. Entre las decisiones importantes que toma la dirección de Big Andy’s se encuentran su política de precios para los diferentes productos y cuánto gastar en publicidad. Para evaluar el efecto de diferentes estructuras de precios y diferentes niveles de gasto en publicidad, Big Andy’s Burger Barn fija diferentes precios y gasta diferentes cantidades en publicidad en diferentes ciudades. A la dirección le interesa saber cómo cambian los ingresos por ventas cuando varía el nivel de gasto en publicidad. ¿Un aumento del gasto en publicidad conlleva un aumento de las ventas? En caso afirmativo, ¿es el aumento de las ventas suficiente para justificar el aumento del gasto en publicidad? La dirección también está interesada en la estrategia de precios. ¿Una reducción de los precios provocará un aumento o una disminución de los ingresos por ventas? Si una reducción del precio sólo conduce a un pequeño aumento de la cantidad vendida, los ingresos por ventas disminuirán (la demanda es inelástica al precio); una reducción del precio que conduzca a un gran aumento de la cantidad vendida producirá un aumento de los ingresos (la demanda es elástica al precio). Esta información económica es esencial para una gestión eficaz.
Modelo Económico
El primer paso es establecer un modelo económico en el que los ingresos por ventas dependen de una o varias variables explicativas. La hipótesis inicial es que los ingresos por ventas están relacionados linealmente con el precio y los gastos de publicidad. El modelo económico es
\[ VENTAS = \beta_0 + \beta_1 PRECIO + \beta_2 PUBLICIDAD \]
donde VENTAS representa los ingresos mensuales por ventas (totales) en una ciudad determinada, PRECIO representa el precio en esa ciudad y PUBLICIDAD es el gasto mensual en publicidad en esa ciudad. Tanto VENTAS como PUBLICIDAD se miden en términos de miles de dólares. Dado que las ventas en las ciudades más grandes tienden a ser mayores que las ventas en las ciudades más pequeñas, nos centramos en las ciudades más pequeñas con poblaciones comparables.
Dado que una hamburguesería vende varios productos -hamburguesas, patatas fritas y batidos- y que cada producto tiene su propio precio, no está claro qué precio debe utilizarse en la ecuacion (1). Lo que necesitamos es algún tipo de precio medio para todos los productos e información sobre cómo cambia este precio medio de una ciudad a otra. Para ello, la dirección ha construido un único índice de precios PRECIO, medido en dólares y centavos, que describe los precios globales en cada ciudad.
Los símbolos restantes en (1) son los parámetros desconocidos \(\beta_1\), \(\beta_2\) y \(\beta_3\) que describen la dependencia de las ventas (VENTAS), del precio (PRECIO) y la publicidad (PUBLICIDAD). Para ser más precisos sobre la interpretación de estos parámetros, pasamos del modelo económico de (1) a un modelo econométrico que hace suposiciones explícitas sobre la forma en que se generan los datos.
Modelo Econométrico
Cuando recojamos los datos de VENTAS, PRECIO y PUBLICIDAD de las franquicias en diferentes ciudades, las observaciones no satisfarán exactamente la relación lineal descrita en la ecuación (1). El comportamiento de los clientes de Andy en las distintas ciudades no será tal que los mismos precios y el mismo nivel de gasto en publicidad conduzcan siempre a los mismos ingresos por ventas. Otros factores que no aparecen en la ecuación y que probablemente afecten a las ventas son el número y el comportamiento de los establecimientos de comida rápida de la competencia, la naturaleza de la población de cada ciudad -su perfil de edad, ingresos y preferencias alimentarias- y la ubicación de las hamburgueserías de Andy’s -cerca de una autopista con mucho tráfico, en el centro de la ciudad, etc.-. Para tener en cuenta estos factores, incluimos un término de error \(e\) en la ecuación, de modo que el modelo se convierte en
\[ VENTAS = \beta_0 + \beta_1 PRECIO + \beta_2 PUBLICIDAD + e \]
Paquetes Utilizados
library("foreign")
library("olsrr")
library("mctest")
library("GGally")andy <- read.dta('andy.dta')
head(andy)## ventas precio publicidad
## 1 73.2 5.69 1.3
## 2 71.8 6.49 2.9
## 3 62.4 5.63 0.8
## 4 67.4 6.22 0.7
## 5 89.3 5.02 1.5
## 6 70.3 6.41 1.3Estimación del Modelo
mod_andy <- lm('ventas ~ precio + publicidad', data=andy)
summary(mod_andy)##
## Call:
## lm(formula = "ventas ~ precio + publicidad", data = andy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.4825 -3.1434 -0.3456 2.8754 11.3049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 118.9136 6.3516 18.722 < 2e-16 ***
## precio -7.9079 1.0960 -7.215 4.42e-10 ***
## publicidad 1.8626 0.6832 2.726 0.00804 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.886 on 72 degrees of freedom
## Multiple R-squared: 0.4483, Adjusted R-squared: 0.4329
## F-statistic: 29.25 on 2 and 72 DF, p-value: 5.041e-10El modelo estimado es:
\[ \hat{VENTAS} = 118.9136 - 7.9079PRECIO + 1.8626 PUBLICIDAD \]
A partir de la información de esta ecuación, se pueden construir fácilmente estimaciones de intervalo o pruebas de hipótesis para cada uno de los \(\beta_j\) de forma similar a un modelo univariado, pero con un cambio en el número de grados de libertad para la distribución \(t\). Al igual que antes, los valores \(t\) y los valores \(p\) de la tabla anterior se refieren a la prueba \(H_0:\beta_j = 0\) contra la alternativa \(H_1:\beta_j \neq 0\) para \(k = 1, 2, 3\).
Procedemos interpretando primero las estimaciones de la ecuación anterior. A continuación, para explicar el cambio de grados de libertad que surge al tener más de una variable explicativa, y para reforzar el material anterior, repasamos las propiedades de muestreo del estimador de mínimos cuadrados, seguido de la estimación por intervalos y las pruebas de hipótesis.
Además de proporcionar información sobre cómo cambian las ventas cuando cambia el precio o la publicidad, la ecuación estimada puede utilizarse para la predicción. Supongamos que el Gran Andy está interesado en predecir los ingresos por ventas para un precio de 5.50 dólares y un gasto en publicidad de 1200 dólares. Incluyendo lugares extradecimales para obtener un cálculo manual preciso, esta predicción es
ventas_p <- 118.9136 - 7.9079*5.5 + 1.8626*1.2
ventas_p## [1] 77.65527Estimador de la Varianza del Error
En el ejemplo de la cadena de hamburguesas, tenemos \(K = 3\). La estimación para nuestra muestra de datos en \[\hat \sigma^2 = \frac{\sum_{i=1}^{75} \hat e_i^2}{N-k} = \frac{1718.943}{75-3}=23.874\]
library('modelr')
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionandy %>%
add_predictions(mod_andy) %>%
summarise(MSE = sum((ventas - pred)^2)/mod_andy$df.residual)## MSE
## 1 23.87421Intervalos de Confianza para \(\beta_j\)
\[ \hat \beta_j - t_{1-\alpha/2,n-p}*ds(\hat \beta_j)\leq \beta_j \leq \hat \beta_j + t_{1-\alpha/2,n-p}*ds(\hat \beta_j) \]
confint(mod_andy)## 2.5 % 97.5 %
## (Intercept) 106.251852 131.575368
## precio -10.092676 -5.723032
## publicidad 0.500659 3.224510Esta estimación de intervalo sugiere que la disminución del precio en 1 dólar conducirá a un aumento de los ingresos medios entre 5723 y 10093 dólares. O, en términos de un cambio de precio cuya magnitud es más realista, una reducción del precio de 10 céntimos conducirá a un aumento de los ingresos medios entre 572 y 1009 dólares. Basándose en esta información, y en el coste de fabricar y vender más hamburguesas, Big Andy puede decidir si procede a una reducción de precios.
Estimamos que un aumento del gasto en publicidad de 1000 dólares conlleva un aumento de los ingresos medios por ventas de entre 501 y 3225 dólares. Este intervalo es relativamente amplio Este intervalo es relativamente amplio; implica que el gasto publicitario adicional podría ser poco rentable (el aumento de los ingresos es inferior a 1.000 dólares) o puede suponer un aumento de los ingresos de más de tres veces el coste de la publicidad. Otra forma de describir esta situación es decir que la estimación puntual b3 = 1,8626 no es no es muy fiable, ya que su error estándar (que mide la variabilidad de muestreo) es relativamente grande.
Intervalo de confianza la 90%
confint(mod_andy, level = 0.90)## 5 % 95 %
## (Intercept) 108.329917 129.497304
## precio -9.734101 -6.081608
## publicidad 0.724180 3.000989Intervalo de confianza para una combinación lineal de parámetros
Big Andy quiere que la próxima semana sea una semana de grandes ventas. Quiere aumentar los gastos de publicidad en 800 dólares y bajar el precio en 40 céntimos. El cambio en las ventas esperadas de la estrategia planificada por Andy es
\[\lambda = -0.4\beta_1+0.8\beta_2\]
library('multcomp')names(coef(mod_andy)) #Extraemos los nombres de los coeficientes (b0,b1,b2)## [1] "(Intercept)" "precio" "publicidad"mod.lh <- glht(mod_andy, linfct = c("-0.4*precio + 0.8*publicidad = 0"))
confint(mod.lh,level=0.9)##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = "ventas ~ precio + publicidad", data = andy)
##
## Quantile = 1.6663
## 90% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## -0.4 * precio + 0.8 * publicidad == 0 4.6532 3.4708 5.8356summary(glht(mod_andy, linfct = c("-0.4*precio + 0.8*publicidad = 0")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = "ventas ~ precio + publicidad", data = andy)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## -0.4 * precio + 0.8 * publicidad == 0 4.6532 0.7096 6.557 7.16e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Estimamos, con un 90% de confianza, que el aumento esperado de las ventas gracias a la estrategia de Big Andy se situará entre 3471 y 5835 dólares.
Prueba de Hipótesis
Prueba de Significancia para el Precio
Planteamiento de Hipótesis
\(H_0:\) \(\beta_1=0\)
\(H_1:\) \(\beta_1\neq 0\)
Nivel de Significancia
\(\alpha = 0.05\)
#####Función de Prueba
ttest <- function(reg, coefnum, val){
co <- coef(summary(reg))
tstat <- (co[coefnum,1]-val)/co[coefnum,2]
print('Valor t = ')
print(tstat)
print('Valor p = ')
2 * pt(abs(tstat), reg$df.residual, lower.tail = FALSE)
}ttest(mod_andy, 2,0)## [1] "Valor t = "
## [1] -7.215241
## [1] "Valor p = "## [1] 4.423997e-10#####Decisión
H0 se rechaza
#####Conclusión
Concluimos que los datos indican que los ingresos por ventas dependen del precio.
Prueba de Significancia para la Publicidad
Planteamiento de Hipótesis
\(H_0:\) \(\beta_2=0\)
\(H_1:\) \(\beta_2\neq 0\)
Nivel de Significancia
\(\alpha = 0.05\)
Función de Prueba
ttest(mod_andy, 3,0)## [1] "Valor t = "
## [1] 2.726283
## [1] "Valor p = "## [1] 0.008038182Decisión
H0 se rechaza
Conclusión
Los datos apoyan la conjetura de que los ingresos están relacionados con los gastos de publicidad.
Prueba de la Elasticidad de la Demanda
Planteamiento de Hipótesis
\(H_0:\) \(\beta_1 \geq 0\) (demanda inelástica)
\(H_1:\) \(\beta_2< 0\) (demanda elástica)
Nivel de Significancia
\(\alpha = 0.05\)
Función de Prueba
ttest2 <- function(reg, coefnum, val){
co <- coef(summary(reg))
tstat <- (co[coefnum,1]-val)/co[coefnum,2]
print('Valor t = ')
print(tstat)
print('Valor p = ')
pt(abs(tstat), reg$df.residual, lower.tail = FALSE)
}ttest2(mod_andy, 2,0)## [1] "Valor t = "
## [1] -7.215241
## [1] "Valor p = "## [1] 2.211999e-10Decisión
H0 se rechaza
Conclusión
Los datos apoyan la proposición de que una reducción del precio provocará un aumento de los ingresos por ventas.
Prueba de la eficacia de la publicidad
Planteamiento de Hipótesis
\(H_0:\) \(\beta_1 \leq 1\)
\(H_1:\) \(\beta_2> 1\)
Nivel de Significancia
\(\alpha = 0.05\)
Función de Prueba
ttest2(mod_andy, 3,1)## [1] "Valor t = "
## [1] 1.262573
## [1] "Valor p = "## [1] 0.1054083Decisión
H0 No se rechaza
Conclusión
La publicidad no es efectiva
Prueba de Hipótesis para una Combinación de Coeficientes
La asesora de marketing de Big Andy afirma que bajar el precio en 20 céntimos será más eficaz para aumentar los ingresos por ventas que aumentar el gasto en publicidad en 500 dólares. En otras palabras, afirma que \(-0.2\beta_1 > 0.5\beta_2\). Andy no desea aceptar esta proposición a menos que pueda ser verificada por datos pasados. Sabe que el cambio estimado en las ventas esperadas por la bajada de precios es
\[-0.2\hat \beta_1 = -0.2 × (-7.9079) = 1.5816\]
y que el cambio estimado en las ventas esperadas por la publicidad extra es
\[0.5 \hat \beta_2 = 0.5 × 1.8626 = 0.9319\]
por lo que la afirmación de la vendedora parece correcta. Sin embargo, quiere establecer si la diferencia \(1.5816 - 0.9319\) podría ser atribuible a un error de muestreo, o si constituye una prueba, a un nivel de significación del 5%, de que \(-0.2\beta_1 > 0.5 \beta_2\). Dado que \(-0.2\beta_1 > 0.5 \beta_2\) puede escribirse como \(-0.2 \beta_1 - 0.5\beta_2> 0\), estamos probando una hipótesis sobre la combinación lineal \(-0.2\beta_1 - 0.5\beta_2\).
Planteamiento de Hipótesis
\(H_0:\) \(-0.2\beta_1 - 0.5\beta_2 \leq 0\) (la afirmación del comercializador no es correcta)
\(H_1:\) \(-0.2\beta_1 - 0.5\beta_2 > 0\) (la afirmación del comercializador es correcta)
Nivel de Significancia
\(\alpha = 0.05\)
Función de Prueba
\[\begin{equation} t = \frac{-0.2\hat \beta_1 - 0.5\hat \beta_2}{ds(-0.2\hat \beta_1 - 0.5\hat \beta_2)} \sim t_{72} \end{equation}\]
\(ds(-0.2\hat \beta_1 - 0.5\hat \beta_2) = \sqrt{\hat{var}(-0.2\hat \beta_1 - 0.5\hat \beta_2)}\)
\(ds(-0.2\hat \beta_1 - 0.5\hat \beta_2) = [(-0.2)^2 \hat \sigma^2_{\hat \beta_1} + (-0.5)^2 \hat \sigma^2_{\hat \beta_2} + 2 \times (-0.2)(-0.5) \hat \sigma_{\hat \beta_1, \hat \beta_2}]^{1/2}\)
\(ds(-0.2\hat \beta_1 - 0.5\hat \beta_2) = [0.04 \times 1.2012 + 0.25 \times 0.4668 + 2\times (-0.0197)]^{1/2}\) \(ds(-0.2\hat \beta_1 - 0.5\hat \beta_2)= 0.4010\)
\[\begin{equation} t = \frac{1.58158 - 0.9319}{0.4010} = 1.622 \end{equation}\]
prueba <- glht(mod_andy, linfct = c("-0.2*precio - 0.5*publicidad = 0"))
summary(prueba)##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = "ventas ~ precio + publicidad", data = andy)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## -0.2 * precio - 0.5 * publicidad == 0 0.6503 0.4010 1.622 0.109
## (Adjusted p values reported -- single-step method)El valor p unilateral es la mitad del bilateral, siempre que la combinación lineal estimada tenga el signo previsto. Para evitar confusiones, recomendamos calcular el valor crítico adecuado o el valor p para las pruebas unilaterales, como se hace a continuación.
pt(abs(1.622), mod_andy$df.residual, lower.tail = FALSE)## [1] 0.05458727Decisión
H0 No se rechaza
Conclusión
A un nivel de significación del 5%, no hay pruebas suficientes para respaldar la afirmación del vendedor.
Evaluación de los residuales
Gráfica de Residuales
Normalidad
Gráfico Q-Q
ols_plot_resid_qq(mod_andy)
Prueba de Normalidad
\(H_0:\) Los residuales siguen una distribución normal
ols_test_normality(mod_andy)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.988 0.7028
## Kolmogorov-Smirnov 0.0673 0.8630
## Cramer-von Mises 5.6577 0.0000
## Anderson-Darling 0.3035 0.5643
## -----------------------------------------------Gráfico residuales vs ajustados
ols_plot_resid_fit(mod_andy)
library("ggfortify")autoplot(mod_andy,which=1:6,ncol = 2,label.size = 2)
Ejemplo 2
La producción de la minería española entre los años 1969 y 1984 expresada en unidades monetarias constantes de 1984 toma los valores de \(X_t\) de la base de datos 1-5.wf1. El empleo del factor trabajo en la producción se expresa mediante la variable \(W_t\) que cuantifica los millones de horas/hombre trabajdas. Para medir el stock de capital o riqueza se utiliza la variable \(C_t\) que representa la potencia instaladas en miles de caballos a vapor. Se trata de estimar las elasticidades del trabajo y el capital respecto de la producción de la industria minera en el periodo 1964-84 considerando como modelo la función de producción de Cobb-Douglas.
\[ X_t = KW_t^aC_t^b \]
library('hexView')ej2 <- readEViews('1-5.wf1')## Warning in readEViews("1-5.wf1"): Skipping boilerplate variable
## Warning in readEViews("1-5.wf1"): Skipping boilerplate variablehead(ej2)## Date ANOS CT T WT XT
## 1 1969-01-01 1969 1141 1 193.5 179.2
## 2 1970-01-01 1970 1241 2 182.8 181.0
## 3 1971-01-01 1971 1357 3 171.7 183.1
## 4 1972-01-01 1972 1465 4 163.4 184.9
## 5 1973-01-01 1973 1562 5 143.4 185.8
## 6 1974-01-01 1974 1742 6 140.4 220.8Aplicamos logarimos para tranformar el modelo a una estructura lineal:
\[ logX_t = logK + alogW_t + blogC_t + e_t \]
Estimación del Modelo
modelo2 <- lm(log(XT) ~ log(WT) + log(CT), data=ej2)
summary(modelo2)##
## Call:
## lm(formula = log(XT) ~ log(WT) + log(CT), data = ej2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16167 -0.04696 -0.01689 0.03870 0.16617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.2946 2.2985 -3.174 0.00733 **
## log(WT) 0.9412 0.3146 2.992 0.01040 *
## log(CT) 1.0641 0.1120 9.503 3.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0938 on 13 degrees of freedom
## Multiple R-squared: 0.9139, Adjusted R-squared: 0.9006
## F-statistic: 68.97 on 2 and 13 DF, p-value: 1.198e-07El modelo estimado es:
\[ logX_t = -7.2946 + 0.9412 logW_t + 1.0641C_t \]
Puede escribirse como:
\[ X_t = 0.00068W_t^{0.94}C_t^{1.06} \]
Evaluación de los residuales
ej2$residual <- residuals(modelo2)
ej2$predicc <- predict(modelo2)Normalidad
Gráfico Q-Q
ols_plot_resid_qq(modelo2)
Prueba de Normalidad
\(H_0:\) Los residuales siguen una distribución normal
ols_test_normality(modelo2)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.9628 0.7124
## Kolmogorov-Smirnov 0.1257 0.9348
## Cramer-von Mises 4.468 0.0000
## Anderson-Darling 0.2964 0.5486
## -----------------------------------------------Gráfico residuales vs ajustados
ols_plot_resid_fit(modelo2)
Gráfico residuales vs variable dependiente
ggplot(ej2,aes(residual,log(XT))) + geom_point()