Modelo de regresión lineal múltiple
El análisis de regresión múltiple es uno de los modelos estadísticos más ampliamente utilizados. Cuando se usa una única variable predictora suelen proporcionarse descripciones inadecuadas del fenómeno, ya que, suelen tenerse series de variables clave que afectan a la variable de respuesta de forma importante y distintiva, un modelo más complejo, que contenga variables predictoras adicionales, suele ser más útil.
Modelo de regresión lineal múltiple.
Se define el modelo de regresión lineal múltiple como se muestra en la Ecuación 1:
\[\begin{align} y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ...+ \beta_{p-1}x_{i,~p-1} + \epsilon_i \end{align} \tag{1}\]
Donde:
\(\beta_0,~\beta_1,...,~\beta_{p-1}\) son parámetros.
\(x_{i1},~x_{i2},...,x_{i,~p-1},\) son constantes conocidas.
\(\epsilon_i\) los términos del error, son independientes y \(N(0,~\sigma^2)\)
\(i=1,2,..,n\)
De manera compacta se puede expresar como se muestra en la Ecuación 2:
\[\begin{align} y_i = \sum_{k=0}^{p-1} \beta_kx_{ik} + \epsilon_i~~;~~ x_{i0} = 1 \end{align} \tag{2}\]
La respuesta esperada \(E(y)\) se muestra en la Ecuación 3:
\[\begin{align} E(y) = \beta_0 + \beta_1x_1 + \beta_2x_2+...+\beta_{p-1}x_{p-1} \end{align} \tag{3}\]
Por lo tanto, el modelo de regresión lineal general con términos del error \(\epsilon_i\) normales, implica que las observaciones \(y_i\) son variables normales e independientes con media E(y_i) dada en Ecuación 3 y con varianza constante \(\sigma^2\).
Modelos de primer orden
Cuando \(x_1, x_2,..., x_{p-1}\) representan \(p-1\) diferentes variables predictoras, el modelo de regresión lineal mostrado en Ecuación 1 se conoce como modelo de primer orden, en el cual no existe interacción entre las variables predictoras.
Variables predictoras cualitativas.
El modelo de regresión lineal múltiple de Ecuación 1 permite, no solo variables predictoras cuantitativas, también puede incluir variables predictoras cualitativas.
Considere un análisis de regresión para predecir la duración de estancia de hospitalización \(y\), basada en la edad \(x_1\) y en el género \(x_2\) del paciente. Se puede definir \(x_2\) como sigue, Ecuación 4:
\[\begin{align} x_2= \left\{ \begin{array}{lcc} 1 & ~si~paciente~es~mujer\\ \\ 0 & si ~paciente~es~hombre \end{array} \right. \end{align} \tag{4}\]
El modelo de regresión lineal de primer orden quedaría como sigue, Ecuación 5:
\[\begin{align} y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \epsilon_i \end{align} \tag{5}\]
Donde:
\[\begin{align} x_{i1} &= edad~del~paciente\\ \\ x_2&= \left\{ \begin{array}{lcc} 1 & ~si~paciente~es~mujer\\ \\ 0 & si ~paciente~es~hombre \end{array} \right. \end{align}\]
La función respuesta para el modelo Ecuación 4 se muestra en la Ecuación 6:
\[\begin{align} E(y) = \beta_0 + \beta_1x_{1} + \beta_2x_{2} \end{align} \tag{6}\]
Para los pacientes hombres, \(x_2 = 0\) el modelo sería, Ecuación 7.
\[\begin{align} E(y) = \beta_0 + \beta_1x_{1} \end{align} \tag{7}\]
Para los pacientes mujeres: \(x_2=1\), el modelo sería, Ecuación 8.
\[\begin{align} E(y) = (\beta_0 + \beta_2) + \beta_1x_{1} \end{align} \tag{8}\]
Ecuación 7 y Ecuación 8 representan rectas paralelas con interceptos diferentes.
En general, representamos una variable cualitativa con \(c\) clases mediante \(c - 1\) variables indicadoras. Por ejemplo, si en el ejemplo de la estancia hospitalaria hay que añadir la variable cualitativa estado de discapacidad como otra variable predictora, se puede representar de la siguiente manera mediante las dos variables indicadoras \(X_3\) y \(X_4\)
\[\begin{align} x_3 &= \left\{ \begin{array}{lcc} 1 & ~si~paciente~no~está~discapacitado\\ \\ 0 & en ~caso~contrario \end{array} \right.\\ \\ x_4 &= \left\{ \begin{array}{lcc} 1 & ~si~paciente~~está~parcialmente~discapacitado\\ \\ 0 & en ~caso~contrario \end{array} \right. \end{align}\]
El modelo de primer orden con edad, género y condición de discapacidad como variables predictoras se muestra en la Ecuación 9:
\[\begin{align} y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \beta_3x_{i3} + \beta_4x_{i4} + \epsilon_i \end{align} \tag{9}\]
Donde:
\[\begin{align} x_{i1} &= edad\\ \\ x_{i2} &= \left\{ \begin{array}{lcc} 1 & ~si~paciente~es~mujer\\ \\ 0 & si~paciente ~es~hombre \end{array} \right.\\ \\ x_{i3} &= \left\{ \begin{array}{lcc} 1 & ~si~paciente~no~está~discapacitado\\ \\ 0 & en ~caso~contrario \end{array} \right.\\ \\ x_{i4} &= \left\{ \begin{array}{lcc} 1 & ~si~paciente~~está~parcialmente~discapacitado\\ \\ 0 & en ~caso~contrario \end{array} \right. \end{align}\]
Modelo de regresión lineal en su forma matricial
Se puede expresar el modelo de Ecuación 1 en su forma matricial definiendo las matrices \(Y\), \(X\), \(\beta\) y \(\epsilon\), Ecuación 10:
\[\underset{n\text{x}1} Y = \begin{bmatrix} y_1 \\ y_2 \\ .\\ .\\ .\\ y_n \end{bmatrix} \underset{n\text{x}p} X = \begin{bmatrix} 1 & x_{11} & x_{12} &...&x_{1,p-1} \\ 1 & x_{21} & x_{22} &...&x_{2,p-1} \\ .&.&.&...&.\\ .&.&.&...&.\\ .&.&.&...&.\\ 1 & x_{n1} & x_{n2} &...&x_{n,p-1} \end{bmatrix} \underset{p\text{x}1} \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ .\\ .\\ .\\ \beta_{p-1} \end{bmatrix} \underset{n\text{x}1} \epsilon = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ .\\ .\\ .\\ \epsilon_{n} \end{bmatrix} \tag{10}\]
Donde:
\(Y\): matriz de variable respuesta.
\(\beta\): vector de parámetros.
\(X\) matriz relacionada con variables represoras-
\(\epsilon\): vector de variables aleatorias independientes con valor esperado \(E(\epsilon) = 0\) y matriz de varianza-covarianza:
\[\underset{n\text{x}n} X = \begin{bmatrix} \sigma^2 & 0 & 0 &...&0 \\ 0 & \sigma^2 & 0 &...& 0 \\ . & . &. & ...& .\\ 0 & 0 & 0 & ...& \sigma^2 \end{bmatrix} = \sigma^2I\]
Estimación de coeficientes de regresión mediante mínimos cuadrados
Se define \(Q\) como:
\[\begin{align} Q = \sum_{i=1}^n (y_i - \beta_0 - \beta_1x_{i1} - ... - \beta_{p-1}x_{i,~p-1})^2 \end{align}\]
Los estimadores por mínimos cuadrados son los valores de \(\beta_0, \beta_1,..,\beta_{p-1}\) que minimizan \(Q\). Se denota al vector de estimadores de los coeficientes de regresión \(\hat{\beta}\), como aparecen en Ecuación 11$
\[\underset{p\text{x}1} {\hat{\beta}} = \begin{bmatrix} \hat{\beta_0}\\ \hat{\beta_1}\\ .\\ .\\ .\\ \hat{\beta_{p-1}}\\ \end{bmatrix} \tag{11}\]
Los estimadores se calculan como se muestra en Ecuación 12
\[\begin{align} \hat{\beta} = (X^{'}X)^{-1} X^{'}Y \end{align} \tag{12}\]
Ejemplo estimación de parámetros regresión lineal múltiple.
Dwaine Studios, Inc. gestiona estudios de fotografía en \(21\) ciudades intermedias. Estos estudios están especializados en retratos de niños. La empresa está considerando expandirse a otras ciudades intermedias y desea investigar si las ventas en una ciudad pueden predecirse a partir del número de personas de \(16\) años o menos en la comunidad y el ingreso per cápita en la comunidad. En la Tabla 1 se muestran los datos de estas variables correspondientes al año más reciente de las \(21\) ciudades en las que opera actualmente Dwaine Studios.
Observación \(i\) | \(x_1\): Población menor 16 años (miles de personas) | \(x_2\): Ingreso per cápita (miles de dólares) | \(y\): ventas (miles de dólares) |
---|---|---|---|
1 | 68.5 | 16.7 | 174.4 |
2 | 45.2 | 16.8 | 164.4 |
3 | 91.3 | 18.2 | 244.2 |
4 | 47.8 | 16.3 | 154.6 |
5 | 46.9 | 17.3 | 181.6 |
6 | 66.1 | 18.2 | 207.5 |
7 | 49.5 | 15.9 | 152.8 |
8 | 52.0 | 17.2 | 163.2 |
9 | 48.9 | 16.6 | 145.4 |
10 | 38.4 | 16.0 | 137.2 |
11 | 87.9 | 18.3 | 241.9 |
12 | 72.8 | 17.1 | 191.1 |
13 | 88.4 | 17.4 | 232.0 |
14 | 42.9 | 15.8 | 145.3 |
15 | 52.5 | 17.8 | 161.1 |
16 | 85.7 | 18.4 | 209.7 |
17 | 41.3 | 16.5 | 146.4 |
18 | 51.7 | 16.3 | 144.0 |
19 | 89.6 | 18.1 | 232.6 |
20 | 82.7 | 19.1 | 224.1 |
21 | 52.3 | 16.0 | 166.5 |
En el Figura 1 se observa la representación gráfica de los datos.
Si ajustamos un plano de regresión, se observa en la Figura 2 A continuación se muestra paso a paso la estimación de los parámetros en la matriz \(\hat{\beta}\)
Definición de matrices \(X\) y \(Y\):
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)) Y
Transpuesta de la matriz \(X\)
<- t(X) XT
\(X^{'}X\)
<-XT %*% X XTX
\((X^{'}X)^{-1}\)
<-solve(XTX) INVXTX
\(X^{'}Y\)
<- XT%*%Y XTY
\(\hat{\beta}\)
<- INVXTX%*%XTY
betas betas
[,1]
[1,] -68.85707
[2,] 1.45456
[3,] 9.36550
De los resultados obtenidos en la matriz \(\hat{\beta}\) se tiene la función de regresión estimada que se muestra en Ecuación 13:
\[\begin{align} \hat{y} = -68,85707 + 1,45456x_1 + 9,3655x_2 \end{align} \tag{13}\]
Ejemplo estimación de parámetros regresión lineal múltiple - RStudio
Para estimar los parámetros usando RStudio
seguimos la siguiente secuencia:
#Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Ajustar el modelo de regresión
<- lm(ventas ~ poblacion + ingreso)
modelo
# Extraer información del modelo
summary(modelo)
Call:
lm(formula = ventas ~ poblacion + ingreso)
Residuals:
Min 1Q Median 3Q Max
-18.4239 -6.2161 0.7449 9.4356 20.2151
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -68.8571 60.0170 -1.147 0.2663
poblacion 1.4546 0.2118 6.868 2e-06 ***
ingreso 9.3655 4.0640 2.305 0.0333 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.01 on 18 degrees of freedom
Multiple R-squared: 0.9167, Adjusted R-squared: 0.9075
F-statistic: 99.1 on 2 and 18 DF, p-value: 1.921e-10
Del procedimiento anterior en la sección de coeficientes se obtienen los datos de los coeficientes de regresión estimados así:
\(\hat{\beta_0}= -68,8571\)
\(\hat{\beta_1} = 1,4546\)
\(\hat{\beta_2} = 9,3655\)
Llegando a la misma función de regresión estimada de Ecuación 13.
Valores ajustados \(\hat{Y}\) y residuales \(e\)
Se denota a la matriz de valores ajustados como \(\hat{Y}\) como se muestra en la Ecuación 14.
\[ \underset {n\text{x}1} {\hat{Y}} = \begin{bmatrix} \hat{y_1} \\ \hat{y_2} \\ .\\ .\\ .\\ \hat{y_n} \\ \end{bmatrix} \tag{14}\]
los valores estimados o ajustados \(\hat{Y}\) se calculan siguiendo la Ecuación 15:
\[\begin{align} \hat{Y} = X\hat{\beta} \end{align} \tag{15}\]
Se denota a la matriz de residuales como \(e\), como se muestra en la Ecuación 16.
\[ \underset {n\text{x}1} e = \begin{bmatrix} e_1 \\ e_2 \\ .\\ .\\ .\\ e_n \\ \end{bmatrix} \tag{16}\]
Los residuales \(e\) se calculan siguiendo la Ecuación 17
\[\begin{align} e &= Y - \hat{Y}\\ e &= Y - X\hat{\beta} \end{align} \tag{17}\]
Ejemplo cálculo de valores ajustados \(\hat{Y}\) y residuales \(e\).
Dwaine Studios, Inc. gestiona estudios de fotografía en \(21\) ciudades intermedias. Estos estudios están especializados en retratos de niños. La empresa está considerando expandirse a otras ciudades intermedias y desea investigar si las ventas en una ciudad pueden predecirse a partir del número de personas de \(16\) años o menos en la comunidad y el ingreso per cápita en la comunidad. En la Tabla 1 se muestran los datos de estas variables correspondientes al año más reciente de las \(21\) ciudades en las que opera actualmente Dwaine Studios.
# Definición de matrices
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Transpuesta de X
<- t(X)
XT
# Transpuesta de X por X
<-XT %*% X
XTX
# Inversa de X transpuesta por X
<-solve(XTX)
INVXTX
# X transpuesta por Y
<- XT%*%Y
XTY
# Betas
<- INVXTX%*%XTY
betas
# Valores ajustados
<- X%*%betas
ajustados
# Residuales
<- Y-ajustados residuales
Ejemplo cálculo de ajustados \(\hat{Y}\) y residuales \(e\) en RStudio
Para obtener los valores ajustados \(\hat{Y}\) y los residuales \(e\) usando RStudio
seguimos la siguiente secuencia:
#Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Ajustar el modelo de regresión
<- lm(ventas ~ poblacion + ingreso)
modelo
# Extraer información de ajustados desde el modelo
<- modelo$fitted.values
ajustados ajustados
1 2 3 4 5 6 7 8
187.1841 154.2294 234.3963 153.3285 161.3849 197.7414 152.0551 167.8666
9 10 11 12 13 14 15 16
157.7382 136.8460 230.3874 197.1849 222.6857 141.5184 174.2132 228.1239
17 18 19 20 21
145.7470 159.0013 230.9870 230.3161 157.0644
# Extraer información de residuales desde el modelo
<- modelo$residuals
residuales residuales
1 2 3 4 5 6
-12.7841146 10.1705737 9.8036764 1.2714690 20.2150722 9.7585779
7 8 9 10 11 12
0.7449178 -4.6666316 -12.3381967 0.3539791 11.5126289 -6.0849209
13 14 15 16 17 18
9.3142995 3.7815611 -13.1132116 -18.4238900 0.6530062 -15.0013134
19 20 21
1.6129777 -6.2160615 9.4356009
Gráficamente podemos observar el comportamiento de \(Y\), \(\hat{Y}\) en la Figura 3
Inferencia sobre el modelo - Análisis de varianza para el modelo de regresión lineal múltiple
Sumas de cuadrados \(SS\)
Suma de cuadrados total \(SST\)
La suma de cuadrados del total \(SST\) se calcula como se muestra en la Ecuación 18.
\[\begin{align} SST = Y^{'} \left[ I - \left( \frac{1}{n} \right) J \right]Y \end{align} \tag{18}\]
La matriz \(J\) es una matriz \(n\text{x}n\) definida en la con todos los elementos \(\left[a_{ij}\right] =1\)
La matriz \(H\) está definida en la Ecuación 19
\[\begin{align} H \underset {n\text{x}n} = X(X^{'}X)^{-1}X^{´} \end{align} \tag{19}\]
A continuación se calculará, paso a paso, la suma de cuadrados total \(SST\) para el ejemplo de Dwaine Studios, Inc con datos de Tabla 1
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Matriz J
<- matrix(c(rep(1,21*21)), ncol=21, nrow=21)
J
# 1/n J
<- (1/21)*J
J_sobre_n
# Matriz Identidad n*n
<- diag(21)
I
# Identidad menos i/n H
<- I-J_sobre_n
I_menos_J_sobre_n
# SST
<- t(Y)%*%I_menos_J_sobre_n%*%Y
SST SST
[,1]
[1,] 26196.21
Suma de cuadrados del error \(SSE\)
La suma de cuadrados del error \(SSE\) se calcula como se muestra en la Ecuación 20.
\[\begin{align} SSE = Y^{'} (I-H)Y \end{align} \tag{20}\]
A continuación se calculará, paso a paso, la suma de cuadrados del error \(SSE\) para el ejemplo de Dwaine Studios, Inc con datos de Tabla 1
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso # Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Matriz H
<- X%*%solve(t(X)%*%X)%*%t(X)
H
#Matriz identidad
<- diag(21)
I
#Matriz identidad menos H
<- I-H
I_menos_H
#SSE
<- t(Y)%*%I_menos_H%*%Y
SSE SSE
[,1]
[1,] 2180.927
Suma de cuadrados del modelo o suma de cuadrados de la regresión \(SSR\)
La suma de cuadrados del modelo \(SSR\) se calcula como se muestra en la Ecuación 21.
\[\begin{align} SSR = Y^{'} \left[H- \left( \frac{1}{n} \right)J \right]Y \end{align} \tag{21}\]
A continuación se calculará, paso a paso, la suma de cuadrados del modelo \(SSR\) para el ejemplo de Dwaine Studios, Inc con datos de Tabla 1
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso # Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Matriz H
<- X%*%solve(t(X)%*%X)%*%t(X)
H
# Matriz J
<- matrix(c(rep(1,21*21)), ncol=21, nrow=21)
J
#Matriz J sobre n
<- J*(1/21)
J_sobre_n
#Matriz H - J sobre n
<- H-J_sobre_n
H_menos_J_sobre_n
#SSR
<- t(Y)%*%H_menos_J_sobre_n%*%Y
SSR SSR
[,1]
[1,] 24015.28
Grados de libertad \(DF\)
Grados de libertad del total \(DF_T\)
Los grados de libertad del total \(DF_T\) se calculan como se muestra en la Ecuación 22
\[\begin{align} DF_T = n-1 \end{align} \tag{22}\]
Grados de libertad del error \(DF_{E}\)
Los grados de libertad del error \(DF_E\) se calculan como se muestra en la Ecuación 23
\[\begin{align} DF_E = n-p \end{align} \tag{23}\]
Grados de libertad del modelo \(DF_{R}\)
Los grados de libertad del modelo \(DF_R\) se calculan como se muestra en la Ecuación 24
\[\begin{align} DF_R = p-1 \end{align} \tag{24}\]
Cuadrados medios \(MS\)
Cuadrados medios del total \(MST\)
Los cuadrados medios del total \(MST\) se calculan como se muestra en la Ecuación 25
\[\begin{align} MST = \frac{SST}{DF_T} \end{align} \tag{25}\]
Cuadrados medios del error \(MSE\)
Los cuadrados medios del error \(MSE\) se calculan como se muestra en la Ecuación 26
\[\begin{align} MSE = \frac{SSE}{DF_E} \end{align} \tag{26}\]
Cuadrados medios del modelo \(MSR\)
Los cuadrados medios del modelo \(MSR\) se calculan como se muestra en la Ecuación 27
\[\begin{align} MSR = \frac{SSR}{DF_R} \end{align} \tag{27}\]
A continuación se calculará, paso a paso, los cuadrados medios para el ejemplo de Dwaine Studios, Inc con datos de Tabla 1, teniendo los datos de \(SST\), \(SSE\) Y \(SSR\) calculados anteriormente
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Grados de libertad
<- length(ventas)-1
DFT <- length(ventas) - 3
DFE <- 3-1
DFR
# Cuadrados medios
<- SST/DFT
MST <- SSE/DFE
MSE <- SSR/DFR
MSR
MST
[,1]
[1,] 1309.81
MSE
[,1]
[1,] 121.1626
MSR
[,1]
[1,] 12007.64
Planteamiento de las hipótesis
Para la realización del análisis de varianza, ANOVA, que representa la significancia global del modelo se plantean las hipótesis como se muestra en la Ecuación 28.
\[\begin{align} H_0&: \beta_1 = \beta_2 =...=\beta_{p-1} = 0 \\ \\ H_1&: \beta_k \neq 0~ para~al~menos~un~k~; k=1,2,...,p-1 \end{align} \tag{28}\]
Estadístico de prueba \(F_0\)
El estadístico de prueba \(F_0\) se calcula como se muestra en la Ecuación 29.
\[\begin{align} F_0 = \frac{MSR}{MSE} \end{align} \tag{29}\]
Rechazo \(H_0\) si:
\[ F_0 > F_{1-\alpha,~p-1,~n-p} \tag{30}\]
A continuación se calculará el estadístico de prueba y el estadístico de referencia para el ejemplo de Dwaine Studios, Inc con datos de Tabla 1, teniendo los datos de \(MSR\) Y \(MSE\) calculados anteriormente.
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Estadístico de prueba F_0
<- MSR/MSE
F_0 F_0
[,1]
[1,] 99.1035
#Cuantil teórico para alpha=0.05
qf(0.05, 2, 18, lower.tail = FALSE)
[1] 3.554557
Para un nivel de significancia \(\alpha = 0.05\). El estadístico de prueba \(F_0= 99.1035\) y el estadístico teórico \(F_{\alpha, ~p-1,~n-p} = 3.55\) por lo que existe evidencia estadística suficiente para rechazar \(H_0\) por lo que al menos un \(\beta_k\neq0\). Podemos concluir que las ventas están relacionadas con la población y el ingreso per capita.
Análisis de varianza en RStudio
A continuación se realiza el ANOVA en RStudio
# Vectores de variables
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
#Modelo
<- lm(ventas~poblacion + ingreso)
modelo
# Análisis de varianza
<- aov(modelo)
anova summary(anova)
Df Sum Sq Mean Sq F value Pr(>F)
poblacion 1 23372 23372 192.896 4.64e-11 ***
ingreso 1 643 643 5.311 0.0333 *
Residuals 18 2181 121
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelo)
Call:
lm(formula = ventas ~ poblacion + ingreso)
Residuals:
Min 1Q Median 3Q Max
-18.4239 -6.2161 0.7449 9.4356 20.2151
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -68.8571 60.0170 -1.147 0.2663
poblacion 1.4546 0.2118 6.868 2e-06 ***
ingreso 9.3655 4.0640 2.305 0.0333 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.01 on 18 degrees of freedom
Multiple R-squared: 0.9167, Adjusted R-squared: 0.9075
F-statistic: 99.1 on 2 and 18 DF, p-value: 1.921e-10
#Cuantil teórico para alpha=0.05
qf(0.05, 2, 18, lower.tail = FALSE)
[1] 3.554557
Las conclusiones son las misma: “Para un nivel de significancia \(\alpha = 0.05\). El estadístico de prueba \(F_0= 99.1035\) y el estadístico teórico \(F_{\alpha, ~p-1,~n-p} = 3.55\) por lo que existe evidencia estadística suficiente para rechazar \(H_0\) por lo que al menos un \(\beta_k\neq0\). Podemos concluir que las ventas están relacionadas con la población y el ingreso per capita.”
Inferencias sobres los \(\beta\)
Los estimadores \(\hat{\beta}\) obtenidos mediante mínimos cuadrados (o máxima verosimilitud) son insesgados, esto es, Ecuación 31:
\[ E \left[ \hat{\beta} \right] = \beta \tag{31}\]
La matriz de varianza-covarianza \(\sigma^2(\hat{\beta})\) se muestra en la Ecuación 32.
\[\begin{align} \sigma^2 (\hat{\beta}) = \begin{bmatrix} \sigma^2 (\hat{\beta_0} ) ~& ~ \sigma(\hat{\beta_0}, \hat{\beta_1})~ &~ ...~ &~ \sigma(\hat{\beta_0}, \hat{\beta_{p-1}}) \\ \sigma(\hat{\beta_1}, \hat{\beta_0}) & \sigma^2 (\hat{\beta_1}) ~& ~ ...~ &~ \sigma(\hat{\beta_1}, \hat{\beta_{p-1}}) \\ . & . & ...&.\\ . & . & ...&.\\ . & . & ...&.\\ \sigma(\hat{\beta_{p-1}}, \hat{\beta_0}) & \sigma (\hat{\beta_{p-1}}, \hat{\beta_1}) ~& ~ ...~ &~ \sigma^2 ( \hat{\beta_{p-1}}) \end{bmatrix} \end{align} \tag{32}\]
La Ecuación 32 está dada por la Ecuación 33:
\[ \sigma^2(\hat{\beta}) = \sigma^2(X^{'}X)^{-1} \tag{33}\]
La matriz de varianza covarianza se puede estimar usando la matriz \(S^2(\hat{\beta})\), Ecuación 34
\[\begin{align} S^2 (\hat{\beta}) = \begin{bmatrix} S^2 (\hat{\beta_0} ) ~& ~ S(\hat{\beta_0}, \hat{\beta_1})~ &~ ...~ &~ S(\hat{\beta_0}, \hat{\beta_{p-1}}) \\ S(\hat{\beta_1}, \hat{\beta_0}) & S^2 (\hat{\beta_1}) ~& ~ ...~ &~ S(\hat{\beta_1}, \hat{\beta_{p-1}}) \\ . & . & ...&.\\ . & . & ...&.\\ . & . & ...&.\\ S\hat{\beta_{p-1}}, \hat{\beta_0}) & S (\hat{\beta_{p-1}}, \hat{\beta_1}) ~& ~ ...~ &~ S^2 ( \hat{\beta_{p-1}}) \end{bmatrix} \end{align} \tag{34}\]
La Ecuación 34 está dada por la Ecuación 35:
\[ S^2(\hat{\beta}) = MSE(X^{'}X)^{-1} \tag{35}\]
Pruebas concercientes a los \(\beta_k\)
Para el modelo de regresión lineal múltiple mostrado en Ecuación 1 se tiene que:
\[\begin{align} \frac{\hat{\beta_k} - \beta_k}{S(\hat{\beta_k})} \sim t_{n-p}~~;~~ k=0,1,2,...,p-1 \end{align} \tag{36}\]
Para todos los \(\beta_k\) parámetros se plantea el siguiente procedimiento:
Planteamiento de hipótesis
\[\begin{align} H_0: \beta_k = 0 \\ \\ H_1: \beta_k \neq 0 \end{align} \tag{37}\]
Estadístico de prueba
Dado lo mostrado en Ecuación 36 el estadístico de prueba se muestra en Ecuación 38:
\[\begin{align} t_{0_k} = \frac{\hat{\beta_k}}{S(\beta_k)} \end{align} \tag{38}\]
Estadístico teórico o de referencia
El estadístico teórico para la decisión corresponde a un cuantil de la distribución \(t\) con probabilidad \(1-\alpha\) y \(n-p\) grados de libertad. Por lo tanto, rechazo \(H_0\) si:
\[\begin{align} |t_{0_k}| > t_{1-\frac{\alpha}{2}~~;~~n-p} \end{align} \tag{39}\]
Inferencias sobre \(\beta_k\) ejemplo.
Se tomará el ejemplo de Dwaine Studios, Inc. En la Tabla 1 se muestran los datos de del problema.
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# MSE para el cálculo de la matriz de varianza-covarianza estimada
<- X%*%solve(t(X)%*%X)%*%t(X)
H <- diag(21)
I <- I-H
I_menos_H <- t(Y)%*%I_menos_H%*%Y
SSE <- 18
DFE <-SSE/DFE
MSE <- MSE[1,1]
MSE_E #MSE como escalar
#X transpuesta
<- t(X)
XT
# X transpuesta por X
<- XT %*% X
XTX
# Inversa de XT*X
<- solve(XTX)
INV_XTX
# Matriz de carianza covarianza estimada
<- MSE_E*INV_XTX
VAR_COV VAR_COV
[,1] [,2] [,3]
[1,] 3602.03467 8.74593958 -241.4229923
[2,] 8.74594 0.04485151 -0.6724426
[3,] -241.42299 -0.67244260 16.5157558
Teniendo la matriz de varianza-covarianza se pueden realizar las pruebas de hipótesis para cada uno de los \(\beta_k\).
Para \(\beta_0\)
\[\begin{align} H_0: \beta_0 =0 \\ H_1: \beta_0 \neq 0 \end{align}\]
El estadístico de prueba se calcularía de la siguiente manera:
\[\begin{align} t_{0_{\beta_0}}= \frac{\hat{\beta_0}}{S(\hat{\beta_0})} \end{align}\]
\(\hat{\beta_0}\) se ha estimado mediante mínimos cuadrados y resultó ser igual a \(\hat{\beta_0}=-68,8571\).
\(S^2(\hat{\beta_0})\) se puede tomar de la matriz varianza-covarianza y con él, se puede calcular \(S(\hat{\beta_0})\):
VAR_COV
[,1] [,2] [,3]
[1,] 3602.03467 8.74593958 -241.4229923
[2,] 8.74594 0.04485151 -0.6724426
[3,] -241.42299 -0.67244260 16.5157558
<- sqrt(VAR_COV[1,1])
s_beta_0 s_beta_0
[1] 60.01695
Por lo que \(S(\hat{\beta_0})=60.01695\) y el estadístico de prueba:
\[\begin{align} t_{0_{\beta_0}}= \frac{-68,8571}{60.01695}=-1,147294 \end{align}\]
El cuantil teórico para la prueba sobre \(\beta_0\), para un nivel de significancia \(\alpha=0.05\) corresponde a:
\[\begin{align} t_{1-\frac{0.05}{2}~, 18} = 2,100922 \end{align}\]
<- qt(1-(0.05/2),18)
cuantil_beta_0 cuantil_beta_0
[1] 2.100922
Como:
\[\begin{align} |t_{0_{\beta_0}}| = 1,147294 \ngtr t_{1-\frac{0.05}{2}~, 18} = 2,100922 \end{align}\]
No existe evidencia estadística suficiente para rechazar \(H_0\) por lo que \(\beta_0 =0\). Esto implica que el modelo no requiere un término independiente significativo (intercepto) y que el valor esperado de la variable respuesta puede depender solo de los regresores.
De manera similar se analiza la significancia de cada \(\beta_k\)
Intervalo de confianza para \(\beta_k\)
Un intervalo de confianza al \((1-\alpha)\%\) para \(\beta_k\) se estima como se muestra en la Ecuación 40
\[\begin{align} \hat{\beta_k} \pm t_{1-\frac{\alpha}{2}~,~n-p}~S(\hat{\beta_k}) \end{align} \tag{40}\]
Ejemplo de cálculo de intervalo de confianza para \({\beta_k}\)
Se tomará el ejemplo de Dwaine Studios, Inc. En la Tabla 1 se muestran los datos de del problema. Se calculará el intervalo de confianza para \(\hat{\beta_1}\)
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Matriz X, Y
<- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3,16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16), nrow = length(ventas), ncol = 3)
X<- matrix(c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5))
Y
# Betas estimados
<- solve(t(X)%*%X)%*%(t(X)%*%Y)
BETAS_EST <- BETAS_EST[2]
BETA_1_EST
# MSE para el cálculo de la matriz de varianza-covarianza estimada
<- X%*%solve(t(X)%*%X)%*%t(X)
H <- diag(21)
I <- I-H
I_menos_H <- t(Y)%*%I_menos_H%*%Y
SSE <- 18
DFE <-SSE/DFE
MSE <- MSE[1,1]
MSE_E
# Matriz de varianza covarianza
<- MSE_E*(solve(t(X)%*%X))
S_2_BETA
# Desviación típica para BETA 1
<- sqrt(S_2_BETA[2,2])
DS_BETA1
# Cuantil teórico alpha =0.05
<- qt(1-(0.05/2),18)
CUANTIL
# Límite inferior
<- BETA_1_EST - (CUANTIL*DS_BETA1)
LIM_INF
#Límite superior
<- BETA_1_EST + (CUANTIL*DS_BETA1)
LIM_SUP
#Intervalo
LIM_INF
[1] 1.009623
LIM_SUP
[1] 1.899497
Ejemplo de cálculo de intervalo de confianza para \(\beta_k\) en RStudio
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Modelo
<- lm(ventas~poblacion+ingreso)
modelo
# Intervalo de confianza
confint(modelo)
2.5 % 97.5 %
(Intercept) -194.9480130 57.233867
poblacion 1.0096226 1.899497
ingreso 0.8274411 17.903560
Predicción de una nueva observación \(Y_{hn}\)
Un intervalo de confianza al \((1-\alpha)\) para una nueva observación \(Y_{hn}\) correspondiente a \(X_h\) se calcula como se muestra en Ecuación 41.
\[\begin{align} \hat{Y_h} \pm t_{1-\frac{\alpha}{2}~,~n-p}~S(pred) \end{align} \tag{41}\]
Donde:
\[\begin{align} S^2(pred) = MSE(1 + X_{hr}(X^{'}X)^{-1}X_{hr}^{'}) \end{align} \tag{42}\]
# Datos
<- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
ventas <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
poblacion <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)
ingreso
# Modelo
<- lm(ventas~poblacion+ingreso)
modelo
# Intervalo de confianza
predict(modelo, newdata = data.frame(poblacion=c(70, 75, 82,90), ingreso=c(17,18,19,17.5)), interval="prediction")
fit lwr upr
1 192.1756 168.0690 216.2822
2 208.8139 184.7072 232.9206
3 228.3613 202.4683 254.2544
4 225.9495 200.1574 251.7417