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.

Tabla 1: Datos ejemplo ventas de 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.

Figura 1: Grafíco de dispersión

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}\)

Figura 2: Plano ajustado

Definición de matrices \(X\) y \(Y\):

X<- 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)
Y <- 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))

Transpuesta de la matriz \(X\)

XT<- t(X) 

\(X^{'}X\)

XTX<-XT %*% X 

\((X^{'}X)^{-1}\)

INVXTX<-solve(XTX)

\(X^{'}Y\)

XTY<- XT%*%Y

\(\hat{\beta}\)

betas<- INVXTX%*%XTY
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Ajustar el modelo de regresión
modelo <- lm(ventas ~ poblacion + ingreso)

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

# Transpuesta de X
XT<- t(X) 

# Transpuesta de X por X
XTX<-XT %*% X 

# Inversa de X transpuesta por X
INVXTX<-solve(XTX)

# X transpuesta por Y
XTY<- XT%*%Y

# Betas
betas<- INVXTX%*%XTY

# Valores ajustados
ajustados <- X%*%betas

# Residuales
residuales <- Y-ajustados

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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Ajustar el modelo de regresión
modelo <- lm(ventas ~ poblacion + ingreso)

# Extraer información de ajustados desde el modelo
ajustados <- modelo$fitted.values
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
residuales <- modelo$residuals
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

Figura 3: Plano ajustado, ajustados y residuales

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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Matriz X, Y
X<- 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)
Y <- 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))

# Matriz J
J <- matrix(c(rep(1,21*21)), ncol=21, nrow=21)

# 1/n J
J_sobre_n <- (1/21)*J

# Matriz Identidad n*n 
I <- diag(21)

# Identidad menos i/n H
I_menos_J_sobre_n <- I-J_sobre_n

# SST
SST <- t(Y)%*%I_menos_J_sobre_n%*%Y
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)
# Matriz X, Y
X<- 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)
Y <- 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))

# Matriz H
H <- X%*%solve(t(X)%*%X)%*%t(X)

#Matriz identidad
I <- diag(21)

#Matriz identidad menos H
I_menos_H <- I-H

#SSE 
SSE <- t(Y)%*%I_menos_H%*%Y
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)
# Matriz X, Y
X<- 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)
Y <- 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))

# Matriz H
H <- X%*%solve(t(X)%*%X)%*%t(X)

# Matriz J
J <- matrix(c(rep(1,21*21)), ncol=21, nrow=21)

#Matriz J sobre n
J_sobre_n <- J*(1/21)

#Matriz H - J sobre n
H_menos_J_sobre_n <- H-J_sobre_n

#SSR
SSR <- t(Y)%*%H_menos_J_sobre_n%*%Y
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Matriz X, Y
X<- 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)
Y <- 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))

# Grados de libertad
DFT <- length(ventas)-1
DFE <- length(ventas) - 3
DFR <- 3-1

# Cuadrados medios
MST <- SST/DFT
MSE <- SSE/DFE
MSR <- SSR/DFR

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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Matriz X, Y
X<- 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)
Y <- 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))

# Estadístico de prueba F_0

F_0 <- MSR/MSE
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

#Modelo
modelo <- lm(ventas~poblacion + ingreso)

# Análisis de varianza
anova <- aov(modelo)
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Matriz X, Y
X<- 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)
Y <- 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))

# MSE para el cálculo de la matriz de varianza-covarianza estimada

H <- X%*%solve(t(X)%*%X)%*%t(X)
I <- diag(21)
I_menos_H <- I-H
SSE <- t(Y)%*%I_menos_H%*%Y
DFE <- 18
MSE <-SSE/DFE
MSE_E <- MSE[1,1]
#MSE como escalar

#X transpuesta
XT <- t(X)

# X transpuesta por X
XTX <- XT %*% X

# Inversa de XT*X
INV_XTX <- solve(XTX)

# Matriz de carianza covarianza estimada
VAR_COV <- MSE_E*INV_XTX
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
s_beta_0 <- sqrt(VAR_COV[1,1])
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}\]

cuantil_beta_0 <- qt(1-(0.05/2),18)
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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Matriz X, Y
X<- 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)
Y <- 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))

# Betas estimados
BETAS_EST <- solve(t(X)%*%X)%*%(t(X)%*%Y)
BETA_1_EST <- BETAS_EST[2]

# MSE para el cálculo de la matriz de varianza-covarianza estimada
H <- X%*%solve(t(X)%*%X)%*%t(X)
I <- diag(21)
I_menos_H <- I-H
SSE <- t(Y)%*%I_menos_H%*%Y
DFE <- 18
MSE <-SSE/DFE
MSE_E <- MSE[1,1]

# Matriz de varianza covarianza
S_2_BETA <- MSE_E*(solve(t(X)%*%X))

# Desviación típica para BETA 1
DS_BETA1 <- sqrt(S_2_BETA[2,2])

# Cuantil teórico alpha =0.05
CUANTIL <- qt(1-(0.05/2),18)

# Límite inferior
LIM_INF <- BETA_1_EST - (CUANTIL*DS_BETA1)

#Límite superior
LIM_SUP <- BETA_1_EST + (CUANTIL*DS_BETA1)

#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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Modelo
modelo <- lm(ventas~poblacion+ingreso)

# 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
ventas <- 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)
poblacion <- 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)
ingreso <- 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)

# Modelo
modelo <- lm(ventas~poblacion+ingreso)

# 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