Examen 2, versión B. Regresión lineal múltiple

El examen se deberá subir al classroom antes de las 11:59 PM del 17 de enero de 2022.

Favor de argumentar con detalle las respuestas.

NOTA. En caso de que se identifiquen respuestas iguales en otros examenes, se procederá a la anulación de los examenes involucrados.

NOTA. Incluir el(los) nombre(s) completo(s) de la(s) persona(s) que está(n) resolviendo los ejercicios.

Usar una confianza de 95% o una significancia de .05 en los casos en donde no se requiera otro nivel de forma explícita. En el caso de realizar alguna transformación de las variables, se tiene que hacer explícita la variable que se usa y la interpretación en las pruebas de hipótesis o intervalos de confianza.

1. (1.5 puntos)

Considere el modelo de regresión

\[ y=\beta_{0} + \beta_{1}x_1 +...+ \beta_px_p + \epsilon, \] y los estimadores obtenidos por mínimos cuadrados escritos en forma matricial \[ \widehat{\bf{\beta}} = \mathbf{ (X^{t}X)^{-1} X^{t}} \bf{y}.\]

Encuentre $V(\widehat{\bf{y}})$, donde $\widehat{\bf{y}}=\bf{X}\widehat{\bf{\beta}}$. Deberá quedar en función de la matriz proyección $\bf H$.
Encuentre $\sum_{i=1}^n V(\widehat{y}_i)$. Notar que esto es equivalente a la traza de la matriz $V(\widehat{\bf{y}})$.

Demostración: $a)$ Para este inciso basta con considerar la siguiente secuencia de igualdades:

\[\begin{align*} V[\hat{\mathbf{Y}}]&=V(\mathbf{X}\hat{\beta})=V[\mathbf{H}y]\\ &\overset{(1)}{=}\mathbf{H}V[\mathbf{y}]\mathbf{H}^{t}=\mathbf{H}(\sigma^{2}\mathbb{I})\mathbf{H}^{t}\\ &\overset{(2)}{=}\sigma^{2}\left(\mathbf{H}\mathbf{H}\right)\\ &\overset{(3)}{=}\sigma^{2}\mathbf{H} \end{align*}\]

donde

$(1)$ Se utilizó una propiedad de la varianza

$(2)$ Se usó que $\bf{H}$ es simétrica. Además podemos conmutar $\sigma^{2}$ con las matrices puesto que ésta es una constante.

$(3)$ Se usó que $\bf{H}$ es idempotente.

$b)$ Para ello notamos que

\[\begin{align*} \sum_{i=1}^{n}V(\hat{y_{i}})&=\sum_{i=1}^{n}V(\mathbf{H} y_{i})\overset{(1)}{=}\sum_{i=1}^{n}\mathbf{H}V(y_{i})\mathbf{H}^{t}\overset{(2)}{=}\sum_{i=1}^{n}\mathbf{H}V(y_{i})\mathbf{H}\\ &=\sum_{i=1}^{n}\mathbf{H}V(\sigma^{2}\mathbb{I})\mathbf{H}=\sum_{i=1}^{n}\sigma^{2}\mathbf{H}\mathbf{H}\\ &\overset{(3)}{=}\sigma^{2}\mathbf{H}\sum_{i=1}^{n}1=n\sigma^{2}\mathbf{H} \end{align*}\]

donde

$(1)$ Aplicamos propiedades de la varianza.

$(2)$ Usamos que $\bf{H}$ es simétrica.

$(3)$ Usamos que $\bf{H}$ es idempotente. Además ocupamos que tanto $\sigma^{2}$ como $\bf{H}$ no dependen del índice de la suma $\square$

2. (1.5 puntos)

Considere el modelo de regresión siguiente: \[ y_i = \beta_0 + \beta_1x_i + \beta_2(\frac{3}{4}x_i^2 -2) + \epsilon_i, \ \ i=1, 2, 3,\]

donde $x_1=0, x_2 = 2, x_3 = -2.$

Defina la matrix diseño $\textbf{X}$ asociada a este modelo. Calcule $\textbf{X}^t\textbf{X}$ y su inversa.

Solución: Para definir la matriz de diseño $\textbf{X}$ debemos ver primero los resultados individuales al evaluar los tres valores de $x_{i}$, con $i=1,2,3$:

\[\begin{align*} x_{1}=0: \ \ \ \ y_{1}&=\beta_0 + \beta_1\cdot (0) + \beta_2\left(\frac{3}{4}\cdot (0)^2 -2\right) + \epsilon_1\\ &=\beta_{0}-2\beta_{2}+\epsilon_{1}\\ x_{2}=2: \ \ \ \ y_{2}&=\beta_0 + \beta_1\cdot (2) + \beta_2\left(\frac{3}{4}\cdot (2)^2 -2\right) + \epsilon_2\\ &=\beta_{0}+2\beta_{1}+\beta_{2}+\epsilon_{2}\\ x_{3}=-2: \ \ \ \ y_{3}&=\beta_0 + \beta_1\cdot (-2) + \beta_2\left(\frac{3}{4}\cdot (-2)^2 -2\right) + \epsilon_3\\ &=\beta_{0}-2\beta_{1}+\beta_{2}+\epsilon_{3} \end{align*}\]

Lo anterior puede verse, de forma equivalente, como una ecuación en términos de matrices

\[ \left(\begin{array}{ccc} y_{1} \\ y_{2} \\ y_{3} \end{array}\right)= \left(\begin{array}{ccc} 1 & 0 & -2\\ 1 & 2 & 1\\ 1 & -2 & 1 \end{array}\right)\cdot \left(\begin{array}{ccc} \beta_{0} \\ \beta_{1} \\ \beta_{2} \end{array}\right)+\left(\begin{array}{ccc} \epsilon_{0} \\ \epsilon_{1} \\ \epsilon_{2} \end{array}\right) \]

y por tanto tenemos que

\[ \textbf{X}=\left(\begin{array}{ccc} 1 & 0 & -2\\ 1 & 2 & 1\\ 1 & -2 & 1 \end{array}\right) \]

Luego, si bien es sencillo hallar $\textbf{X}^{t}$ de manera “manual”, ocuparemos python para efectuar dicho cálculo pues también con esta herramienta calcularemos $\textbf{X}^{t}\textbf{X}$ y $(\textbf{X}^{t}\textbf{X})^{-1}$. De tal manera

# Importamos el paquete necesario
import numpy as np

# Matriz X
x = np.array([[1,0,-2], [1,2,1], [1,-2,1]])

x_t = np.transpose(x) 
print(x_t)

## [[ 1  1  1]
##  [ 0  2 -2]
##  [-2  1  1]]

Se procede después a calcular $\textbf{X}^{t}\textbf{X}$:

x_tx = np.dot(x_t, x)

print(x_tx)

## [[3 0 0]
##  [0 8 0]
##  [0 0 6]]

y $(\textbf{X}^{t}\textbf{X})^{-1}$:

x_tx_inv = np.linalg.inv(x_tx)

print(x_tx_inv)

## [[0.33333333 0.         0.        ]
##  [0.         0.125      0.        ]
##  [0.         0.         0.16666667]]

Por lo tanto

\[ \textbf{X}^{t}\textbf{X}=\left(\begin{array}{ccc} 3 & 0 & 0\\ 0 & 8 & 0\\ 0 & 0 & 6 \end{array}\right) \ \ \ \textrm{y} \ \ \ (\textbf{X}^{t}\textbf{X})^{-1}=\left(\begin{array}{ccc} \frac{1}{3} & 0 & 0\\ 0 & \frac{1}{8} & 0\\ 0 & 0 & \frac{1}{6} \end{array}\right) \]

Dé las expresiones de los estimadores por mínimos cuadrados ordinarios de $\beta_0$, $\beta_1$ y $\beta_2$: $\widehat{\beta}_0$, $\widehat{\beta}_1$ y $\widehat{\beta}_2$. Deberán ser una expresión en términos de $y_1,y_2, y_3$.

Solución: Utilizaremos que $\widehat{\beta}=\left(\textbf{X}^{t}\textbf{X}\right)^{-1}\textbf{X}^{t}\textbf{y}$ es el estimador obtenido por el método de los mínimos cuadrados (por ecuación (102) de las notas de clase). Para calcular dicho vector nuevamente utilizaremos python:

# El siguiente vector reemplazara al vector 'y' en los calculos
y = np.array( [1, 1, 1] )

# Primer calculo: ( (X^t)X )^-1 * X^t
calc1 = np.dot(x_tx_inv, x_t)

# Calculo del estimador btea "gorrito"
beta_est = np.dot(calc1,y)
print(beta_est)

## [1. 0. 0.]

lo que nos indica que

\[ \left(\begin{array}{ccc} \widehat{\beta_{0}} \\ \widehat{\beta_{1}} \\ \widehat{\beta_{2}} \end{array}\right)=\left(\begin{array}{ccc} y_{1} \\ 0 \\ 0 \end{array}\right) \]

Muestre que los estimadores por mínimos cuadrados ordinarios del modelo reducido cuando se supone $\beta_2 = 0$ no se alteran, es decir, que $\widehat{\beta}_0^*=\widehat{\beta}_0$ y $\widehat{\beta}_1^*=\widehat{\beta}_1$, donde $\widehat{\beta}_0^*$ y $\widehat{\beta}_1^*$ son los estimadores por mínimos cuadrados del modelo \[ y_i = \beta_0^* + \beta_1^*x_i + \epsilon_i^*, \ \ i=1, 2, 3.\]

Demostración: Suponiendo que $\beta_{2}=0$ obtenemos el modelo reducido

\[ y_{i}=\beta_{0}^\ast+\beta_{1}^{\ast}x_{i}+\epsilon_{i}^\ast, \ \ \ \ i=1,2,3 \]

y queremos probar que $\widehat{\beta}_0^*=\widehat{\beta}_0$ y $\widehat{\beta}_1^*=\widehat{\beta}_1$. Para ello tendremos que

\[\begin{align*} x_{1}=0: \ \ \ \ y_{1}&=\beta_{0}^\ast+\epsilon_{1}^\ast\\ x_{2}=2: \ \ \ \ y_{2}&=\beta_{0}^\ast+2\beta_{1}^\ast+\epsilon_{2}^\ast\\ x_{3}=-2: \ \ \ \ y_{3}&=\beta_{0}^\ast-2\beta_{1}^\ast + \epsilon_{3}^\ast \end{align*}\]

y por ende tenemos la ecuación matricial

\[ \left(\begin{array}{ccc} y_{1} \\ y_{2} \\ y_{3} \end{array}\right)= \left(\begin{array}{ccc} 1 & 0 \\ 1 & 2 \\ 1 & -2 \end{array}\right)\cdot \left(\begin{array}{cc} \beta_{0}^\ast \\ \beta_{1} ^\ast \end{array}\right)+\left(\begin{array}{ccc} \epsilon_{0}^\ast \\ \epsilon_{1}^\ast \\ \epsilon_{2}^\ast \end{array}\right) \]

y a partir de ahí

\[ \textbf{X}=\left(\begin{array}{ccc} 1 & 0 \\ 1 & 2 \\ 1 & -2 \end{array}\right) \]

Procederemos a obtener $\widehat{\beta}=\left(\textbf{X}^{t}\textbf{X}\right)^{-1}\textbf{X}^{t}\textbf{y}$ de manera totalmente análoga al del inciso anterior, donde nuevamente los cálculos serán hecho en python:

# Matriz de disegno x del modelo reducido
rx = np.array([[1,0], [1,2], [1,-2]])

# X transpuesta
rx_t = np.transpose(rx)

# x transpuesta por x
prod_aux = np.dot(rx_t, rx)

# inversa de (x transpuesta por x)
prod_aux_inv = np.linalg.inv(prod_aux) 

# inversa de (x transpuesta por x) por transpuesta de x
prod_aux_inv_xt = np.dot(prod_aux_inv, rx_t)

# nuevamente usamos el vector auxiliar:
ry = ( [1,1,1] )

rbeta_est = np.dot(prod_aux_inv_xt, ry)
print(rbeta_est)

## [1. 0.]

Con lo cual obtenemos que

\[ \left(\begin{array}{ccc} \widehat{\beta_{0}^{\ast}} \\ \widehat{\beta_{1}^\ast} \end{array}\right)=\left(\begin{array}{ccc} y_{1} \\ 0 \end{array}\right) \]

de donde se comprueba que

\[ \widehat{\beta_{0}^{\ast}}=\widehat{\beta_{0}} \ \ \ \ \ \textrm{y} \ \ \ \ \ \widehat{\beta_{1}^{\ast}}=\widehat{\beta_{1}} \hspace{6cm}\square \]

3. (2 puntos)

La Compañía Kenton Food desea comparar 4 diferentes diseños de empaque de un nuevo cereal. Veinte tiendas, con aproximadamente igual volumen de ventas y perfil de clientes, fueron seleccionadas como unidades experimentales. A cada una de las tiendas se le asignó uno de los empaques de forma aleatoria, de manera que cada empaque fuera asignado a 5 tiendas distintas. Las ventas, en número de casos, fueron observadas durante un período de estudio de 2 semanas:

ventas	12	10	15	19	11	11	17	16	14	15	23	20	18	17	27	33	22	26	28
empaque	1	1	1	1	1	2	2	2	2	2	3	3	3	3	4	4	4	4	4

Un incendio ocurrió en una de las tiendas durante el período de estudio y dado que esto cambia las condiciones de venta con respecto a las otras tiendas se decidió eliminar la medición de esa tienda. El número de ventas de esa tienda se excluye de la tabla anterior.

Asuma que se cumplen todos los supuestos de un problema tipo ANOVA.

Presente un gráfico para describir los datos, por ejemplo, un boxplot por tipo de empaque.
Ajuste un modelo de regresión lineal múltiple adecuado. Indique de acuerdo a los parámetros del modelo la expresión del número de ventas promedio por cada tipo de empaque y dé estimaciones puntuales.
Escriba las hipótesis que se contrastan con la tabla ANOVA, calcule ésta e interprete. Use $\alpha=.05$.
¿Se puede considerar que el diseño del empaque afecta las ventas promedio? Use $\alpha=.05$. Argumente indicando con claridad qué hipótesis se están contrastando en términos de los parámetros del modelo ajustado.
Realicé la prueba de hipótesis simultánea asociada a la igualdad de las ventas promedio entre todos los posibles pares de diferentes empaques. Use $\alpha=.05$. Interprete los resultados.
Suponga que los ejecutivos de la empresa tienen la sospecha de que el diseño de empaque 4 es el que aumenta las ventas en comparación con el resto de empaques. Realice una prueba de hipótesis para argumentar en favor o en contra de esta hipótesis de acuerdo con los datos observados. Use $\alpha=.05$

Solución:

Antes de empezar importaremos los datos

rm(list = ls(all.names = TRUE))
gc()

##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2144638  115    3489024  186  3489024  186
## Vcells 3624841   28    8388608   64  5185992   40

Y<-c(12,10,15,19,11,11,17,16,14,15,23,20,18,17,27,33,22,26,28)
X<-c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,4,4,4,4,4)

Verificamos si la variable X es de tipo factor

is.factor(X)

## [1] FALSE

Como vemos que no es de tipo factor la tranformamos y hacemos un dataframe

X1<-as.factor(X)
levels(X1)

## [1] "1" "2" "3" "4"

Vemos que tiene 4 niveles

Datos<-data.frame(X1,Y)

Mostramos algunas estadisticas de los datos

##  X1          Y     
##  1:5   Min.   :10  
##  2:5   1st Qu.:14  
##  3:4   Median :17  
##  4:5   Mean   :19  
##        3rd Qu.:22  
##        Max.   :33

Podemos ver que las ventas del empaque 4 son mayores que la de losd demás

boxplot(Y ~ X1, data = Datos, col = "white", outline=FALSE)
stripchart(Y ~ X1, data = Datos,method = "jitter",pch = 19,col = 2:4, vertical = TRUE, add = TRUE)

De aqui podemos ver que las medianas de los empaques se reflejan de manera ascendente, es decir, la mediana del empaque 1 es la mas baja, le sigue la del segundo y así consecutivamente

Primero ajustamos el modelo:

fit=lm(Y ~ X1, data = Datos)
summary(fit)

## 
## Call:
## lm(formula = Y ~ X1, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -5.20  -1.95  -0.20   1.50   5.80 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    13.40       1.45    9.23  1.4e-07 ***
## X12             1.20       2.05    0.58    0.568    
## X13             6.10       2.18    2.80    0.013 *  
## X14            13.80       2.05    6.72  6.9e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.2 on 15 degrees of freedom
## Multiple R-squared:  0.788,  Adjusted R-squared:  0.746 
## F-statistic: 18.6 on 3 and 15 DF,  p-value: 2.58e-05

Con este modelo tenemos lo siguiente: \[\mathbb{E}(Y|X_{1}) = \beta_{0} + \beta_{1}X_{12} + \beta_{2}X_{13} + \beta_{3}X_{14}\]

Por lo que tendremos:

$\mathbb{E}(Y|X_{1} = X_{11}) = \beta_{0}$
$\mathbb{E}(Y|X_{1} = X_{12}) = \beta_{0} + \beta_{1}$
$\mathbb{E}(Y|X_{1} = X_{13}) = \beta_{0} + \beta_{2}$
$\mathbb{E}(Y|X_{1} = X_{14}) = \beta_{0} + \beta_{3}$

Así, veamos las estimaciones puntuales:

(b0<-fit$coefficients[1])

## (Intercept) 
##          13

(b1<-fit$coefficients[2])

## X12 
## 1.2

(b2<-fit$coefficients[3])

## X13 
## 6.1

(b3<-fit$coefficients[4])

## X14 
##  14

De lo anterior tenemos de manera puntual

(mu1<-b0)

## (Intercept) 
##          13

(mu1<-b0+b1)

## (Intercept) 
##          15

(mu1<-b0+b2)

## (Intercept) 
##          20

(mu1<-b0+b3)

## (Intercept) 
##          27

Veamos que la prueba de la tabla ANOVA sirve para hacer este test

summary(fit)

## 
## Call:
## lm(formula = Y ~ X1, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -5.20  -1.95  -0.20   1.50   5.80 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    13.40       1.45    9.23  1.4e-07 ***
## X12             1.20       2.05    0.58    0.568    
## X13             6.10       2.18    2.80    0.013 *  
## X14            13.80       2.05    6.72  6.9e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.2 on 15 degrees of freedom
## Multiple R-squared:  0.788,  Adjusted R-squared:  0.746 
## F-statistic: 18.6 on 3 and 15 DF,  p-value: 2.58e-05

La prueba F asociada a la tabla ANOVA indica que se rechaza $H_{0}$ por lo que $b1 \neq 0$ o $b2 \neq 0$ o $b3 \neq 0$ ya que tenemos un p-value menor a la significancia que es 0.05

Primero veamos como son las esperanzas en términos de los parmetros

\[\begin{align*} \mu_{1} &= \beta_{0} \\ \mu_{2} &= \beta_{0} + \beta_{1} \\ \mu_{3} &= \beta_{0} + \beta_{2} \\ \mu_{4} &= \beta_{0} + \beta_{3} \end{align*}\]

Si asumimos que no hay diferencia tenemos:

\[\begin{align*} &\mu_{1} = \mu_{2} = \mu_{3} = \mu_{4} \\ \Rightarrow &\beta_{0} = \beta_{0} + \beta_{1} = \beta_{0} + \beta_{2} = \beta_{0} + \beta_{3} \\ \Rightarrow &\beta_{1} = \beta_{2} = \beta_{3} = 0 \\ \therefore &H_{0} : \beta_{1} = 0 \hspace{.1cm}, \hspace{.1cm} \beta_{2} = 0 \hspace{.2cm} y \hspace{.2cm} \beta_{3} = 0 \hspace{.2cm} vs \hspace{.2cm} H_{a} : \beta_{1} \neq 0 \hspace{.1cm}, \hspace{.1cm} \beta_{2} \neq 0 \hspace{.2cm} y \hspace{.2cm} \beta_{3} \neq 0 \end{align*}\]

library(multcomp)
K=matrix(c(0,1,0,0,0,0,1,0,0,0,0,1), ncol=4, nrow=3, byrow=TRUE)
m=c(0,0,0)
summary(glht(fit, linfct=K, rhs=m),test=Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0      1.2
## 2 == 0      6.1
## 3 == 0     13.8
## 
## Global Test:
##      F DF1 DF2   Pr(>F)
## 1 18.6   3  15 2.58e-05

De aqui obtenemos un p-value igual a $2.58$ x $10^{-5}$, lo que es menor a $\alpha = 0.05$, por lo que rechazamos $H_{0}$ lo que quiere decir que no tenemos informacion suficiente que respalde la hipotesis de que el diseño del empaque las ventas promedio.

Tendremos las siguientes hipótesis

Cuando el empaque uno es igual al empaque dos:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{11}) = \mathbb{E}(Y|X_{1}=X_{12}) \\ \Rightarrow & \beta_{0} = \beta_0 + \beta_{1} \\ \Rightarrow & \beta_{1} = 0 \\ \Rightarrow & H_{0}: \beta_{1} =0 \end{align*}\]

Cuando el empaque uno es igual al empaque tres:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{11}) = \mathbb{E}(Y|X_{1}=X_{13}) \\ \Rightarrow & \beta_{0} = \beta_0 + \beta_{2} \\ \Rightarrow & \beta_{2} = 0 \\ \Rightarrow & H_{0}: \beta_{2} =0 \end{align*}\]

Cuando el empaque uno es igual al empaque cuatro:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{11}) = \mathbb{E}(Y|X_{1}=X_{14}) \\ \Rightarrow & \beta_{0} = \beta_0 + \beta_{3} \\ \Rightarrow & \beta_{3} = 0 \\ \Rightarrow & H_{0}: \beta_{3} =0 \end{align*}\]

Cuando el empaque dos es igual al empaque tres:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{12}) = \mathbb{E}(Y|X_{1}=X_{13}) \\ \Rightarrow & \beta_{0} + \beta_{1} = \beta_0 + \beta_{2} \\ \Rightarrow & \beta_{1} = \beta_{2} \\ \Rightarrow & H_{0}: \beta_{1} - \beta_{2} =0 \end{align*}\]

Cuando el empaque dos es igual al empaque cuatro:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{12}) = \mathbb{E}(Y|X_{1}=X_{14}) \\ \Rightarrow & \beta_{0} + \beta_{1} = \beta_0 + \beta_{3} \\ \Rightarrow & \beta_{1} = \beta_{3} \\ \Rightarrow & H_{0}: \beta_{1} - \beta_{3} =0 \end{align*}\]

Cuando el empaque tres es igual al empaque cuatro:

\[\begin{align*} &\mathbb{E}(Y|X_{1}=X_{13}) = \mathbb{E}(Y|X_{1}=X_{14}) \\ \Rightarrow & \beta_{0} + \beta_{2} = \beta_0 + \beta_{3} \\ \Rightarrow & \beta_{2} = \beta_{3} \\ \Rightarrow & H_{0}: \beta_{2} - \beta_{3} =0 \end{align*}\]

K=matrix(c(0,1,0,0,0,0,1,0,0,0,0,1,0,1,-1,0,0,1,0,-1,0,0,1,-1), ncol=4, nrow=6, byrow=TRUE)
m=c(0,0,0,0,0,0)
summary(glht(fit, linfct=K, rhs=m, alternative="two.sided"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Y ~ X1, data = Datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0     1.20       2.05    0.58    0.935    
## 2 == 0     6.10       2.18    2.80    0.058 .  
## 3 == 0    13.80       2.05    6.72   <0.001 ***
## 4 == 0    -4.90       2.18   -2.25    0.154    
## 5 == 0   -12.60       2.05   -6.13   <0.001 ***
## 6 == 0    -7.70       2.18   -3.53    0.014 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Los p-value nos indican que los empaques uno, dos y tres son iguales entre si y todos son diferentes al empaque 4, así concluimos que el que realmente representa una diferencia es el empaque 4

Para verificar si el empaque 4 es el que aumenta las vetas realizaremos las siguientes pruebas simultaneas:

-Empaque cuatro mejor que empaque uno:

\[\beta_{3} > 0\]

-Empaque cuatro mejor que empaque dos:

\[\beta_{3} - \beta_{1} > 0\]

-Empaque cuatro mejor que empaque tres:

\[\beta_{3} - \beta_{2} > 0\] Por lo que tendremos la siguiente prueba de hipotesis: \[H_{0}: \beta_{3} \leq 0 \hspace{.1cm}, \hspace{.1cm} \beta_{3} - \beta_{1} \leq 0 \hspace{.2cm} y \hspace{.2cm} \beta_{3} - \beta_{2} \leq 0 \hspace{.2cm} vs \hspace{.2cm} H_{a}:\beta_{3} > 0 \hspace{.1cm}, \hspace{.1cm} \beta_{3} - \beta_{1} > 0 \hspace{.2cm} y \hspace{.2cm} \beta_{3} - \beta_{2} > 0\]

K=matrix(c(0,0,0,1,
0,-1,0,1,
0,0,-1,1), ncol=4, nrow=3, byrow=TRUE)
m=c(0,0,0)
summary(glht(fit, linfct=K, rhs=m, alternative="greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Y ~ X1, data = Datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0    13.80       2.05    6.72 <0.001 ***
## 2 <= 0    12.60       2.05    6.13 <0.001 ***
## 3 <= 0     7.70       2.18    3.53  0.004 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

De aqui podemos ver que los p-value son menores que la significancia, por lo tato se rechaza $H_{0}$, con lo que concluimos que en efecto, el empaque cuatro aumenta las ventas.

4. (2 puntos)

Una institución de investigación realiza un estudio para analizar los efectos de un nuevo tratamiento para controlar los niveles altos de ansiedad. Para eso consideran un puntaje (a mayor valor mayores niveles de ansiedad) y definen un conjunto experimental con 120 individuos que en ese puntaje presentan valores similares al inicio del estudio, 60 son hombres y 60 mujeres. En el mercado se sabe que hay otro tratamiento que se usa comúnmente para este fin, de manera que de forma aleatoria han divido a los 120 en tres grupos: 40 a los que no se aplicó ningún tratamiento (control), 40 a los que se aplicó el tratamiento actual (Trat1) y 40 a los que se aplicó el nuevo tratamiento (Trat2); 20 hombres y 20 mujeres en cada grupo. Los datos se presentan en el archivo Ex4B.csv.

Los investigadores sospechan que para el nuevo tratamiento podría existir un efecto diferenciado de acuerdo con el sexo, por lo que consideran conveniente incluir esta variable en el análisis.

(para este ejercicio no se requiere verificar supuestos del modelo, asuma que se cumplen)

Realice un análisis descriptivo de los datos. Observe que hay dos variables categóricas, incluya entonces un boxplot para cada posible combinación de niveles que se pueden observar en esas dos variables categóricas (boxplot(Puntaje~Trat+Sexo, …)). Comente lo que observe.

Solución: Comenzamos por cargar los datos

datos<- read.csv(file = 'Ex4B.csv')

Procedemos a revisar la estructura de los datos.

datos1 <- data.frame(as.factor(datos$Sexo), as.factor(datos$Trat), datos$Puntaje) 
colnames(datos1) <- c('Sexo', 'Trat', 'Puntaje')
str(datos1)

## 'data.frame':    120 obs. of  3 variables:
##  $ Sexo   : Factor w/ 2 levels "Hombre","Mujer": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Trat   : Factor w/ 3 levels "Control","Trat1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Puntaje: num  9.19 11.42 12.63 7.48 11.64 ...

Con esto podemos ver que tenemos un conjuntos con dos variables categóricas, tenemos la variable Sexo con dos niveles y la variable Trat con tres niveles.

Ahora veremos algunas estadísticas de los datos:

summary(datos1)

##      Sexo         Trat       Puntaje    
##  Hombre:60   Control:40   Min.   : 2.9  
##  Mujer :60   Trat1  :40   1st Qu.: 6.3  
##              Trat2  :40   Median : 7.7  
##                           Mean   : 8.1  
##                           3rd Qu.:10.2  
##                           Max.   :14.8

Podemos observar que la muestra está balanceada en cada grupo, además podríamos pensar que la variable Puntaje sigue una distribución simétrica, pues la media y la mediana tienen valores parecidos, para poder corroborar esto, agregaremos el histograma de la variable Puntaje.

hist(datos1$Puntaje, main = 'Histograma de la variable Puntaje',
     xlab = 'Puntaje', ylab = 'Frecuencia', col = 'forestgreen')

Con el cual podemos ver que NO tenemos una distribución simétrica. Procederemos a construir el boxplot con nuestros datos para las categorías que tenemos:

with(datos1, boxplot(Puntaje ~ Trat + Sexo, col = 'yellow'))

En este boxplot podemos observar que para el grupo control los resultados en promedio para hombres y mujeres son parecidos, sin embargo, tenemos un promedio ligeramente mayor para mujeres. En el tratamiento uno podemos observar que si disminuye en promedio la variable puntaje, algo que se podría esperar de un tratamiento ya en el mercado, a su vez podemos observar que en promedio el valor en Mujeres es ligeramente menor.

Por último, para el nuevo tratamiento podemos observar que, en promedio, tiene una menor eficacia que el tratamiento 1, asimismo tiene un efecto notablemente menor en mujeres ya que el puntaje promedio en mujeres es mayor. En general, para el tratamiento 2, tenemos que en promedio hay mejores resultados para los hombres ya que en promedio el puntaje es menor que para mujeres. Con estos datos podríamos pensar que la hipótesis de que el puntaje se ve descrito por el sexo de la persona podría ser cierto.

Considere un modelo de regresión donde en las covariables se incluyan las dos variables categóricas de forma individual y también su interacción. De acuerdo con los parámetros de ese modelo, ajuste la regresión y dé la expresión del puntaje promedio para cada valor de las variables categóricas: $E(puntaje|Trat=k, Sexo=l)$, con $k \in \{Control, Trat1, Trat2 \}$ y $l \in \{Hombre, Mujer\}$; así como estimaciones puntuales.

Solución: Comenzando con la variable Sexo, como tiene dos niveles es necesario agregar una variable $X$, de forma que $X = 1$ si es mujer, $X = 0$ en otro caso. De esta forma incluimos el efecto de esta variable.

datos1$X = (datos1$Sexo == 'Mujer')*1

Ahora para la variable Trat, tiene tres niveles por lo que hay que agregar 2 variables dicotómicas, $Y$,$Z$ tales que $Y = 1$ si la persona recibio Trat 1, $Y = 0$ en otro caso. $Z = 1$ si la persona recibio Trat 2, $Z = 0$ en otro caso.

datos1$Y = (datos1$Trat == 'Trat1')*1
datos1$Z = (datos1$Trat == 'Trat2')*1

Ahora para añadir el modelo de interacción agregaremos dos variables dicotómicas U, V tales que $U = 1$ si recibio Trat 1 y es mujer, $U = 0$ en otro caso. $V = 1$ si recibió Trat 2 y es mujer, $V = 0$ en otro caso.

datos1$U = (datos1$Trat == 'Trat1' & datos1$Sexo == 'Mujer')*1
datos1$V = (datos1$Trat == 'Trat2' & datos1$Sexo == 'Mujer')*1

Por lo que podemos quedarnos con el modelo

\[ Puntaje = \beta_{0} + \beta_{1}X + \beta_{2}Y + \beta_{3}Z + \beta_{4}U + \beta_{5}V + \epsilon, \]

Su interpretación en términos de la esperanza de la variable puntaje está dada por:

\[\begin{align*} &E(Puntaje|Trat = Control, Sexo = Hombre) = \beta_{0},\\ &E(Puntaje|Trat = Control, Sexo = Mujer) = \beta_{0} + \beta_{1},\\ &E(Puntaje|Trat = Trat1, Sexo = Hombre) = \beta_{0} + \beta_{2},\\ &E(Puntaje|Trat = Trat1, Sexo = Mujer) = \beta_{0} + \beta_{1} + \beta_{2} + \beta_{4},\\ &E(Puntaje|Trat = Trat2, Sexo = Hombre) = \beta_{0} + \beta_{3},\\ &E(Puntaje|Trat = Trat2, Sexo = Mujer) = \beta_{0} + \beta_{1} + \beta_{3} + \beta_{5}. \end{align*}\]

Por último ajustamos el modelo de regresión líneal multiple:

fit1 <- lm(Puntaje ~ X + Y + Z + U + V, data = datos1)
summary(fit1)

## 
## Call:
## lm(formula = Puntaje ~ X + Y + Z + U + V, data = datos1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.143 -0.906 -0.187  0.756  4.000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   10.624      0.307   34.66  < 2e-16 ***
## X              0.807      0.433    1.86   0.0653 .  
## Y             -4.490      0.433  -10.36  < 2e-16 ***
## Z             -3.291      0.433   -7.59  9.5e-12 ***
## U             -1.639      0.613   -2.67   0.0086 ** 
## V             -0.310      0.613   -0.51   0.6143    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.4 on 114 degrees of freedom
## Multiple R-squared:  0.736,  Adjusted R-squared:  0.724 
## F-statistic: 63.5 on 5 and 114 DF,  p-value: <2e-16

Con el valor de los coeficientes del modelo, damos las estimaciones puntuales para todas las esperanzas de la variable Puntaje:

\[\begin{align*} &E(Puntaje|Trat = Control, Sexo = Hombre) = 10.624,\\ &E(Puntaje|Trat = Control, Sexo = Mujer) = 11.431, &E(Puntaje|Trat = Trat1, Sexo = Hombre) = 6.134,\\ &E(Puntaje|Trat = Trat1, Sexo = Mujer) = 5.302,\\ &E(Puntaje|Trat = Trat2, Sexo = Hombre) = 7.333,\\ &E(Puntaje|Trat = Trat2, Sexo = Mujer) = 7.83 \end{align*}\]

Escriba las hipótesis que se contrastan con la tabla ANOVA, calcule ésta e interprete. Use $\alpha=.05$.

Solución: La prueba que queremos contrastar es:

\[\begin{align*} &H0 : \beta_{1} = \beta_{2} = \ldots = \beta_{5} = 0 vs\\ &Ha : \beta_{1} \neq 0 \ \ o \ \ \beta_{2} \neq 0 \ \ o \ldots o \ \ \beta_{5} \neq 0, \end{align*}\]

Esta prueba nos dice que queremos contrastar que ninguna de las variables utilizadas sirve para el modelo contra que al menos una de estas variables nos aporta información

Veamos la prueba ANOVA

K0 <- matrix(c(0,1,0,0,0,0,
               0,0,1,0,0,0,
               0,0,0,1,0,0,
               0,0,0,0,1,0,
               0,0,0,0,0,1), ncol = 6, nrow = 5, byrow = T)
summary(multcomp::glht(fit1, linfct = K0, rhs = c(0,0,0,0,0)),
        test = multcomp::Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    0.807
## 2 == 0   -4.490
## 3 == 0   -3.291
## 4 == 0   -1.639
## 5 == 0   -0.310
## 
## Global Test:
##      F DF1 DF2   Pr(>F)
## 1 63.5   5 114 2.43e-31

Con un $p-value$ de $2.43e-31$ podemos decir que contamos con suficiente evidencia para rechazar la hipotesis nula, es decir, alguna de las variables que agregamos aporta información a la variable puntaje.

iV. ¿Se puede considerar que el sexo tiene un efecto en el puntaje, es decir, al menos para un tratamiento existe un efecto diferenciado en el puntaje derivado del sexo de los individos? Use una prueba F con $\alpha=.05$. Argumente. Aquí: $H_0:E(puntaje|Trat=k, Sexo=Hombre)=E(puntaje|Trat=k, Sexo=Mujer) \forall k \in \{Control, Trat1, Trat2 \}$. En caso de no rechazar $H_0$ considere el modelo reducido eliminando la variable Sexo; pero si se rechaza $H_0$, considere una prueba simultánea que ayude a identificar para qué tratamiento se puede considerar que el sexo tiene un efecto, con los resultados de esa prueba reduza el modelo si es posible.

Solución: Traduciendo nuestra hipótesis nula en términos de los coeficientes de la regresión, esta queda expresada como

\[ H_0 : \beta_{1} = 0, \ \ \beta_{1} + \beta_{4} = 0, \ \ \beta_{1} + \beta_{5} = 0 \]

Realizamos ahora la prueba F para esta hipótesis:

K1 <- matrix(c(0,1,0,0,0,0,
               0,1,0,0,1,0,
               0,1,0,0,0,1), ncol = 6, nrow = 3, byrow = T)
summary(multcomp::glht(fit1, linfct = K1, rhs = c(0,0,0)),
        test = multcomp::Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    0.807
## 2 == 0   -0.833
## 3 == 0    0.497
## 
## Global Test:
##      F DF1 DF2 Pr(>F)
## 1 2.82   3 114  0.042

Teniendo una significancia a=0.025 podemos decir que no tenemos suficiente información para rechazar la hipótesis nula. Por lo que podemos concluir que el sexo no tiene efecto en el modelo.Por lo que nos quedamos con el modelo

\[ Puntaje = \beta_{0} + \beta_{2}Y + \beta_{3}Z + e. \]

En caso de que en el inciso anterior se haya reducido el modelo, ajuste de nuevo la regresión y dé la expresión del puntaje promedio para cada valor en las variables categóricas: $E(puntaje|Trat=k, Sexo=l)$, con $k \in \{Control, Trat1, Trat2 \}$ y $l \in \{Hombre, Mujer\}$; así como estimaciones puntuales.

Solución: Como en el inciso anterior se redujo el modelo, ya que no tomó en cuenta la variable sexo. Tenemos que:

\[ Puntaje = \beta_{0} + \beta_{1}Y + \beta_{2}Z + e, \]

y su interpretación en términos de la esperanza de la variable Puntaje ahora está dada por:

\[\begin{align*} &E(Puntaje|Trat = Control, Sexo = Hombre) = \beta_{0}, &E(Puntaje|Trat = Control, Sexo = Mujer) = \beta_{0}, &E(Puntaje|Trat = Trat1, Sexo = Hombre) = \beta_{0} + \beta_{1}, &E(Puntaje|Trat = Trat1, Sexo = Mujer) = \beta_{0} + \beta_{1}, &E(Puntaje|Trat = Trat2, Sexo = Hombre) = \beta_{0} + \beta_{2}, &E(Puntaje|Trat = Trat2, Sexo = Mujer) = \beta_{0} + \beta_{2}. \end{align*}\]

Luego, volvemos a ajustar el modelo de regresión lineal múltiple

fit2 <- lm(Puntaje ~ Y + Z, data = datos1)
summary(fit2)

## 
## Call:
## lm(formula = Puntaje ~ Y + Z, data = datos1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.546 -0.990 -0.153  0.687  3.796 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11.027      0.222    49.7   <2e-16 ***
## Y             -5.309      0.314   -16.9   <2e-16 ***
## Z             -3.445      0.314   -11.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.4 on 117 degrees of freedom
## Multiple R-squared:  0.716,  Adjusted R-squared:  0.711 
## F-statistic:  148 on 2 and 117 DF,  p-value: <2e-16

Obteniendo así las siguientes estimaciones puntuales

E(Puntaje|Trat = Control, Sexo = Hombre) = 11.027,
E(Puntaje|Trat = Control, Sexo = Mujer) = 11.027,
E(Puntaje|Trat = Trat1, Sexo = Hombre) = 5.718,
E(Puntaje|Trat = Trat1, Sexo = Mujer) = 5.718,
E(Puntaje|Trat = Trat2, Sexo = Hombre) = 7.582,
E(Puntaje|Trat = Trat2, Sexo = Mujer) = 7.582.

Realice una prueba de hipótesis para argumentar en favor o en contra de la hipótesis: el nuevo tratamiento tiene el mejor desempeño. Use $\alpha=.05$

Solución: Queremos poner a prueba la hipótesis nula $H_0$ dada por

\[ E(Puntaje|Trat = Trat1)\geq E(Puntaje|Trat = Trat2) \]

la cual se traduce en términos de los coeficientes del modelo de regresión como

\[ H_0 : 0 \geq \beta_{1} - \beta_{2} \]

Para ambos casos ya que llegamos a que no influía el sexo, con lo que solo faltaría realizar la prueba de hipótesis:

K2 <- matrix(c(0,-1,1), ncol = 3, nrow = 1, byrow = T)
summary(multcomp::glht(fit2, linfct = K2, rhs = c(0),
                       alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Puntaje ~ Y + Z, data = datos1)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(>t)    
## 1 <= 0    1.864      0.314    5.94 1.5e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

De esta manera podemos decir que tenemos suficiente evidencia para rechazar H0, es decir, podemos afirmar que en promedio el tratamiento 1 tiene mejores resultados.

Realice una prueba de hipótesis para argumentar en favor o en contra de la hipótesis: el tratamiento actual tiene el mejor desempeño. Use $\alpha=.05$ Nota. Suponga que tiene dos variables categóricas, una con $K$ niveles y la otra con $J$ niveles. Para usar el modelo de regresión con interacciones se requiere incluir $K-1$ y $J-1$ variables binarias asociadas a los efectos principales de los niveles, además de $(K-1)\times(J-1)$ variables binarias asociadas a las interacciones. Las variables de las interacciones se construyen como el producto de las $K-1$ y $J-1$ variables binarias que se introducen en el modelo.

Solución: Queremos poner a prueba la hipótesis nula H0 dada por

E(Puntaje|Trat = Trat2) ???

E(Puntaje|Trat = Trat1)

la cual se traduce en términos de los coeficientes del modelo de regresión como

H0 : 0 ??? ??1 ??? ??2.

Para ambos casos ya que llegamos a que no influía el sexo, con lo que solo faltaría realizar la prueba de hipótesis:

K2 <- matrix(c(0,1,-1), ncol = 3, nrow = 1, byrow = T)
summary(multcomp::glht(fit2, linfct = K2, rhs = c(0),
                       alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Puntaje ~ Y + Z, data = datos1)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)
## 1 <= 0   -1.864      0.314   -5.94      1
## (Adjusted p values reported -- single-step method)

De esta manera podemos decir que NO tenemos suficiente evidencia para rechazar H0, es decir, podemos afirmar que en promedio el tratamiento 1 tiene mejores resultados.

5. (2 puntos)

Suponga que una empresa farmacéutica está ofreciendo al gobierno un nuevo medicamento para tratar a pacientes con la enfermedad Covid-19. El costo del medicamento es considerable y para tomar una buena decisión se han acercado a usted para analizar los datos que ha compartido la empresa farmacéutica. El archivo Ex5.csv contiene la información: $Ant$ es el número total de anticuerpos, $Trat$ es una variable con dos niveles dependiendo si se aplicó o no el nuevo medicamento. Se sabe que tener mayores anticuerpos evita que se desarrolle una versión grave de la enfermedad y la empresa afirma que eso se logra al aplicar el medicamento, pues los pacientes que recibieron el medicamento tienen más anticuerpos que los que sólo recibieron placebo. También se sabe que la generación de anticuerpos es diferente dependiendo de la edad de los individuos y se sospecha que eso también podría afectar la efectividad del medicamento, así que al diseñar el experimento se seleccionaron al azar 100 personas de 300 que presentaban síntomas leves al iniciar el cuadro de la enfermedad a los que se les administró el medicamento, al resto se les dió sólo seguimiento. En todos los pacientes se capturó la edad y se procuró tener pacientes en el rango entre 16 y 60 años en ambos grupos. No se sospecha de otro aspecto que pudiera modificar la evaluación del medicamento.

(para este ejercicio no se requiere verificar supuestos del modelo, asuma que se cumplen)

Realice un análisis descriptivo de los datos considerando tanto la información de la edad como de la administración o no del medicamento.

Solución:

Comenzamos por cargar y ver los datos

datos <- read.csv("Ex5.csv")
head(datos)

además

str(datos)

## 'data.frame':    300 obs. of  3 variables:
##  $ Ant : num  24.3 21.2 22.2 15.3 17.7 ...
##  $ Trat: chr  "Med" "Med" "Med" "Med" ...
##  $ Edad: int  29 51 34 55 57 18 39 55 40 36 ...

levels(datos$Trat)

## NULL

# convertimos la variable Trat a tipo factor:
datos$Trat <- as.factor(datos$Trat) 

levels(datos$Trat)

## [1] "Control" "Med"

Es preciso para realizar el análisis descriptivo de los datos que podamos visualizarlos, para ello

plot(datos$Edad, datos$Ant, xlab="Edad", ylab="Anticuerpos",pch=19, col = c("red", "blue")[datos$Trat])
legend("bottomleft", levels(datos$Trat),
       col = c("red", "blue"), pch = 19, inset = 0.01,  pt.cex=1.5,cex = .9, y.intersp = 1.3 , bty="n")

Tanto para control como para med podemos observar cierta nube de puntos además de observar posibles rectas para cada uno. De tal manera puede verse cierta relación lineal entre la edad y los anticuerpos. Algo que se ve claramente es que el número de anticuerpos disminuye conforme la edad aumenta, independientemente del si el tratamiento fue suministrado o no, aunque se ve mayor decaimiento en aquellas personas mayores que no recibieron el medicamento.

Finalmente, vemos en el gráfico que el número de anticuerpos es menor, sin importar la edad, en aquellas personas a las que no se les aplicó el medicamento.

Por otro lado, también observamos

boxplot(Ant~Trat, data=datos, xlab="Tratamiento", ylab="Anticuerpos")

Vemos que la media respectiva de med es mayor a la media de control, es decir, parece ser que en promedio el número de anticuerpos es mayor en aquellos pacientes que recibieron el medicamento versus aquellos que no lo recibieron. Pero recordemos que lo dicho anteriormente es sólo una idea y no puede ser tomado como algo contundente debido a que en el modelo tenemos la presencia de otra variable.

Ajuste un modelo adecuado para evaluar la efectividad del medicamento ajustando por la edad de los pacientes. Es decir, un modelo que incluya como explicativas las variables edad, la binaria asociada a la administración del medicamento y la interacción obtenida como el producto de éstas dos.

Solución: Asumimos que se cumplen los supuestos del modelo.

Buscamos pues realizar el ajuste del siguiente modelo:

\[ \textrm{Anticuerpos}= \beta_{0}+\beta_{1}( \textrm{Trat}\$Med)+\beta_{2}( \textrm{Edad})+\beta_{3}( \textrm{Trat\$Med})( \textrm{Edad})+ \epsilon \]

En consecuencia tenemos que

\[ \mathbb{E}[y_{i}]=\beta_{0}+\beta_{1}\cdot x_{ig_{1}}+\beta_{2}\cdot x_{2}+\beta_{3}\cdot x_{1g_{1}}x_{2} \]

donde $x_{ig_{1}}=1$ si grupo 1 (el grupo referente al medicamente) y $x_{ig_{1}}=0$ en cualquier otro caso; además $x_{2}$ es la variable que hace referencia a la edad.

Observación: Líneas arriba vimos que la referencia en la variable Trat es control.

Luego ajustamos y vemos cierta información sobre dicho ajuste

fit=lm(Ant ~ Trat*Edad, data=datos)
summary(fit)

## 
## Call:
## lm(formula = Ant ~ Trat * Edad, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.986 -1.913 -0.063  1.884  9.704 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   25.4308     0.6874   37.00   <2e-16 ***
## TratMed        0.9351     1.1764    0.79     0.43    
## Edad          -0.3091     0.0172  -17.96   <2e-16 ***
## TratMed:Edad   0.1607     0.0295    5.45    1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3 on 296 degrees of freedom
## Multiple R-squared:  0.713,  Adjusted R-squared:  0.71 
## F-statistic:  245 on 3 and 296 DF,  p-value: <2e-16

Lo primero que debemos preguntarnos es si nuestro modelo tiene sentido, para ello nos ayudaremos de la prueba asociada a la tabla Anova, donde $p-value=2e-16<0.05$. Por ende, podemos asumir que nuestro modelo si tiene sentido.

Luego, basándonos en la información del summary anterior tenemos las siguientes expresiones de las esperanzas

\[\begin{align*} &\mathbb{E}[y|\textrm{Control}, x_{2}]=\beta_{0}+\beta_{2}x_{2}\\ &\mathbb{E}[y|\textrm{Med}, x_{2}]=\beta_{0}+\beta_{1}+\beta_{2}x_{2}+\beta_{3}x_{2}=\beta_{0}+\beta_{1}+(\beta_{2}+\beta_{3})x_{2} \end{align*}\]

Finalmente realizamos una prueba sobre la igualdad de pendientes donde contrastaremos

\[ H_{0}:\beta_{3}=0 \ \ \ vs \ \ H_{a}:\beta_{3}\neq0 \]

library(multcomp)
K = matrix(c(0,0,0,1), ncol=4, nrow=1, byrow=TRUE) 

m=0

summary(multcomp::glht(fit, linfct=K, rhs=m), test=Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    0.161
## 
## Global Test:
##      F DF1 DF2   Pr(>F)
## 1 29.7   1 296 1.05e-07

donde $p-value=1e-07<0.05$, rechazándose $H_{0}$ y por ende es plausible asumir que $\beta_{3}\neq 0$. Con lo anterior concluimos que el modelo no puede reducirse, dado que sólo estamos trabajando con dos variables, por lo que podemos continuar trabajando con nuestro modelo.

De acuerdo con el modelo ajustado, indique las expresiones asociadas a la relación de la generación promedio de anticuerpos con la edad en a) el grupo control y b) en el grupo que recibe el medicamento.

Solución: En términos matemáticos:

la generación promedio de anticuerpos está expresada como $\mathbb{E}[y]$.
la generación promedio de anticuerpos con la edad en el grupo control está expresada como $\mathbb{E}[y| \textrm{control}, x_{2}]=\beta_{0}+\beta_{2}\cdot x_{2}$.
la generación promedio de anticuerpos con la edad en el grupo que recibe el medicamento está expresada como $\mathbb{E}[y| \textrm{control}, x_{2}]=\beta_{0}+\beta_{1}+(\beta_{2}+\beta_{3})\cdot x_{2}$.

¿Se puede decir que la edad afecta de la misma forma la generación de anticuerpos en el grupo control que en el grupo que recibe el medicamento? Realice una prueba de hipótesis apropiada e interprete.

Solución: Equivalentemente nos preguntamos: ¿el efecto de la variable $x_{2}$ es el mismo en las diferentes clases de $x_{1}$?

Observación: $x_{1}=1$ si hablamos de med y $x_{1}=0$ en cualquier otro caso.

Ahora, planeamos la siguiente hipótesis:

$H_{0}$: la edad afecta de la misma forma la generación de anticuerpos en el grupo control que en el grupo que recibe el medicamento

o equivalentemente

\[ H_{0}: \mathbb{E}(y|Control, x_{2}+1)-\mathbb{E}(y|Control,x_{2})=H_{0}: \mathbb{E}(y|Med, x_{2}+1)-\mathbb{E}(y|Med,x_{2}) \]

finalmente se tiene la equivalencia $H_{0}: \beta_{3}=0$ y por ende la prueba de hipótesis a contrastar es

\[ H_{0}: \beta_{3}=0 \ \ \ \ vs \ \ \ \ H_{a}: \beta_{3}\neq 0 \]

de donde

K = matrix(c(0,0,0,1), ncol=4, nrow=1, byrow=TRUE) 
m=0
summary(multcomp::glht(fit, linfct=K, rhs=m), test=Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    0.161
## 
## Global Test:
##      F DF1 DF2   Pr(>F)
## 1 29.7   1 296 1.05e-07

obtenieno que $p-value=1.049e-07<0.05$ por lo cual se rechaza $H_{0}$, es decir, es plausible asumir que

la edad NO afecta de la misma forma la generación de anticuerpos en el grupo control que en el grupo que recibe el medicamento

Comente sobre el ajuste del modelo incluyendo la interpretación de cada uno de los coeficientes.

Solución: Comenzaremos por ver nuevamente el summary del modelo

summary(fit)

## ## Call: ## lm(formula = Ant ~ Trat * Edad, data = datos) ## ## Residuals: ## Min 1Q Median 3Q Max ## -7.986 -1.913 -0.063 1.884 9.704 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 25.4308 0.6874 37.00 <2e-16 *** ## TratMed 0.9351 1.1764 0.79 0.43 ## Edad -0.3091 0.0172 -17.96 <2e-16 *** ## TratMed:Edad 0.1607 0.0295 5.45 1e-07 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3 on 296 degrees of freedom ## Multiple R-squared: 0.713, Adjusted R-squared: 0.71 ## F-statistic: 245 on 3 and 296 DF, p-value: <2e-16

Como ya se dijo anteriormente, de acuerdo a la prueba asocida a la tabla Anova podemos concluir que nuestro modelo tiene sentido. Luego, como $R^{2}=0.713$ podemos decir que en general nuestro modelo es bueno.

Procedemos después a la interpretación de los coeficientes:

coefficients(fit)

## (Intercept) TratMed Edad TratMed:Edad ## 25.43 0.94 -0.31 0.16

$\beta_{0}=25.43$: Cantidad promedio de anticuerpos para aquellas personas que tienen 16 años y a las cuales no se les ha suministrado el medicamento.

$\beta_{1}=0.94$: Cantidad promedio de anticuerpos que ganan aquellas personas de 16 años a las que se les ha suministrado el medicamento.

$\beta_{2}=-0.31$: Cantidad promedio de anticuerpos que pierden las personas por cada año más de edad.

$\beta_{3}=0.16$: Cantidad promedio de anticuerpos que se evitan perder cada año mas de edad a las personas que se les suministró el medicamento.

Argumente en contra o a favor de la afirmación: “El medicament0 funciona aumentando el número de anticuerpos para todos los pacientes entre 16 y 60 años”. Se puede apoyar de pruebas de hipótesis o intervalos de confianza simultáneos.

Solución: Procederemos a realizar intervalos de confianza simultáneos.

Primero:

# Creamos un "intervalo" para las edades desde los 16 hasta los 60 edad <- seq(from = 16, to = 60, by = 0.5) length(edad)

## [1] 89

procedemos después a programar las bandas para las rectas respectivas de control y Med

# Banda para la recta del Control # E(Y|Control,x2)=b0+b2*edad KA <- cbind(1, 0, edad, 0) # Banda para la recta del Medicamento # E(Y|Med,x2)=(b0+b1)+(b2+b3)*edad KB <- cbind(1, 1, edad, edad) # Juntamos la información en una sola matriz K=rbind(KA, KB)

Después utilizamos la función glht para definir las combinaciones lineales de la matriz k para calcular los intervalos de confianza de éstas mediante confint:

fitE <- glht(fit, linfct = K) # Dado que son bastantes intervalos bajamos el nivel de confianza fitci <- confint(fitE, level = 0.90)

plot(datos$Edad, datos$Ant, pch=19, col = c("#FC7702", "#02E5FC")[datos$Trat], xlab = "Edad", ylab = "Anticuerpos") # Intervalo para el Control lines(edad, coef(fitE)[1:89], col="red") lines(edad, fitci$confint[1:89,"upr"], col="red") lines(edad, fitci$confint[1:89,"lwr"], col="red") # Intervalo para Med lines(edad, coef(fitE)[90:178], col="blue") lines(edad, fitci$confint[90:178,"upr"], col="blue") lines(edad, fitci$confint[90:178,"lwr"], col="blue")

Vemos que las rectas no se intersectan en ningún punto entre los 16 y 60 años. De lo anterior se deduce que el número de anticuerpos aumenta para aquellas personas que se les ha suministrado el medicamento. Por tanto, lo anterior argumenta a favor de la afirmación de los investigadores.

6. (2 puntos)

Considere los datos del archivo Ex6.csv.

Considere un modelo de regresión lineal múltiple con las covariables $X_1$ a $X_6$ sin ninguna interacción. La variable dependiente es $Y$. Verifique los supuestos del modelo. En caso de que alguno no se cumpla, realice transformaciones convenientes hasta que obtenga un modelo que parezca cumplir con los supuestos del modelo de regresión lineal múltiple.

Con las variables transformadas en el inciso anterior, realice una selección de variables. Justifique su respuesta.

Solución:

Antes de empezar, cargamos los datos

rm(list = ls(all.names = TRUE)) gc()

## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 2422154 129 4541850 243 3489024 186 ## Vcells 4125515 32 8388608 64 6842354 52

datos <- read.csv(file = 'Ex6.csv')

Consideramos el modeo de regresion multiple y checamos el nombre de las varaibles y si son de tipo categoricas, en cuyo caso verificamos si están codificadas como factores.

str(datos)

## 'data.frame': 400 obs. of 7 variables: ## $ X3: num 4.17 4.42 2.78 2.09 4.21 ... ## $ X4: num 5.63 5 5.28 4 4.38 ... ## $ X5: num -5.452 -3.748 -1.611 -0.998 -1.535 ... ## $ X6: num -0.536 1.875 0.068 0.486 2.287 ... ## $ X1: chr "A3" "A3" "A3" "A3" ... ## $ X2: chr "B1" "B2" "B3" "B4" ... ## $ Y : num 1.38e+31 1.96e+20 2.22e+24 1.15e+14 9.35e+15 ...

Ahora, checamos algunas estadisticas descriptivas.

summary(datos)

## X3 X4 X5 X6 X1 ## Min. :-2.2 Min. :2.2 Min. :-10.7 Min. :-6.7 Length:400 ## 1st Qu.: 1.9 1st Qu.:4.4 1st Qu.: -5.2 1st Qu.:-1.5 Class :character ## Median : 3.3 Median :5.0 Median : -3.1 Median :-0.1 Mode :character ## Mean : 3.2 Mean :5.1 Mean : -3.2 Mean :-0.1 ## 3rd Qu.: 4.5 3rd Qu.:5.7 3rd Qu.: -1.1 3rd Qu.: 1.4 ## Max. : 8.5 Max. :8.3 Max. : 4.4 Max. : 6.6 ## X2 Y ## Length:400 Min. :0.00e+00 ## Class :character 1st Qu.:3.04e+16 ## Mode :character Median :1.86e+23 ## Mean :3.66e+67 ## 3rd Qu.:5.87e+30 ## Max. :1.46e+70

De aqui, podemos ver que los datos están balanceados con respecto a las variales categoricas, además, tanto la variable $X_{4}$ como $Y$ toman valores positivos. Otra cosa a considerar de la variable $Y$ es que tiene un rango muy amplio, además su distribución es bastante asimetrica ya que la mediana tien un valor del orden $10^{23}$ y la media del orden $10^{67}$, asi que lo que podemos hacer para empezar es transfromar la variable $Y$ a una escala logarítmica:

datos$logY <- log(datos$Y)

Ahora presentamos un resumen grafico del conjunto de datos

ggpairs(datos)

En este grafico podemos ver que las variables continuas sean significativas en el modelo, excepto la variable $X_4$ y un poco la variable $X_6$, esto lo analizaremos más a detalle despues.

Ahora ajustamos el modelo

fit1 <- lm(logY ~ X1 + X2 + X3 + X4 + X5 + X6, data = datos) summary(fit1)

## ## Call: ## lm(formula = logY ~ X1 + X2 + X3 + X4 + X5 + X6, data = datos) ## ## Residuals: ## Min 1Q Median 3Q Max ## -5.363 -1.854 -0.737 0.883 18.361 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -62.47547 0.93162 -67.06 < 2e-16 *** ## X1A2 35.16379 0.44143 79.66 < 2e-16 *** ## X1A3 19.47794 0.44189 44.08 < 2e-16 *** ## X1A4 19.95520 0.44452 44.89 < 2e-16 *** ## X2B2 -5.83115 0.43991 -13.26 < 2e-16 *** ## X2B3 -5.80323 0.43984 -13.19 < 2e-16 *** ## X2B4 -3.68654 0.43837 -8.41 7.9e-16 *** ## X3 -0.12391 0.07762 -1.60 0.11 ## X4 20.31372 0.15545 130.68 < 2e-16 *** ## X5 0.00841 0.05409 0.16 0.88 ## X6 -3.16612 0.07728 -40.97 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3.1 on 389 degrees of freedom ## Multiple R-squared: 0.986, Adjusted R-squared: 0.985 ## F-statistic: 2.69e+03 on 10 and 389 DF, p-value: <2e-16

La prueba ANOVA del modelo nos dice que al menos alguna de las variables $X_{1}$ a $X_{6}$ es útil en el modelo, así que podemos continuar con la verificación de supuestos. Empezaremos con el supuesto de linealidad, para el cual nos apoyamos de la gráfica de los residuales observados.

plot(fit1,1)

Este gráfico nos deja ver que no se cumple el supuesto, para verificar cuales son las variables que causan esto, realizamos el siguiente gráfico:

car::residualPlots(fit1)

## Test stat Pr(>|Test stat|) ## X3 -1.17 0.24 ## X4 56.95 <2e-16 *** ## X5 0.80 0.43 ## X6 -1.35 0.18 ## Tukey test 18.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos notar que la variable que causa el problema es $X_4$.

Ahora podemos verificar el supuesto de homocedasticidad, para esto veamos la siguiente gráfica

plot(fit1,3)

En este grafico parece que no se cumle el supuesto, por lo que ahora realizaremos la prueba de Breusch-Pagan:

lmtest::bptest(fit1)

## ## studentized Breusch-Pagan test ## ## data: fit1 ## BP = 13, df = 10, p-value = 0.2

Ahora la prueba nvc:

car::ncvTest(fit1)

## Non-constant Variance Score Test ## Variance formula: ~ fitted.values ## Chisquare = 16, Df = 1, p = 6e-05

Con el p-value obtenido en la prueba Breusch-Pagan no se rechaza $H_{0}$ con una significancia $\alpha = 0.05$; por lo tanto se cumple el supuesto de homocedasticidad, sin embargo el p-value obtenido en la prueba nvc indica que se rechaza la hipótesis de homocedasticidad, por lo que no es posible asumir el supuesto de homocedassticidad.

Ahora veamos el supuesto de normalidad, para lo que construiremos el Q-Q plot:

plot(fit1,2)

En este gráfico podemos ver que no se cumple el supuesto de normalidad.

Ahora solo queda verificar el supuesto de aleatoriedad. Lo haremos con el autocorrelograma:

Datosfit=broom::augment(fit1) acf(Datosfit$.std.resid, main = 'Autocorrelograma')

De aqui podemos decir que se puede asumir el supuesto de aleatoriedad, pero lo cotejaremos con una prueba de aleatoriedad

randtests::runs.test(Datosfit$.std.resid)

## ## Runs Test ## ## data: Datosfit$.std.resid ## statistic = -0.5, runs = 196, n1 = 200, n2 = 200, n = 400, p-value = ## 0.6 ## alternative hypothesis: nonrandomness

Con el p-value obtenido no se rechaza $H_{0}$, por lo que podemos asumir que se cumple el supuesto de aleatoriedad.

Ahora que sabemos cuales son los supuesto que no se cumplen podemos trabajar sobre estos, comencemos con el supuesto de linealidad.

Para este supuesto ya sabemos que el problema se presenta en la variable $X_{4}$ y ya que esta variable toma valores positivos, intentaremos arreglar el problema con una transformación Box-Tidwell:

car::boxTidwell(logY~X4, ~X1+X2+X3+X5+X6 , data = datos)

## MLE of lambda Score Statistic (z) Pr(>|z|) ## 2 54 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## iterations = 2

Con base en esto, proponemos la tramsformación : \[Z_{4} = X^{2}_{4}\] Por lo que ajustamos nuestro nuevo modelo:

datos$Z4 <- (datos$X4)^2 fit2 <- lm(logY ~ X1 + X2 + X3 + Z4 + X5 + X6, data = datos)

Verifiquemos de nuevo el supuesto de linealidad:

plot(fit2, 1)

Podemos ver que esta vez si se puede asumir el suepuesto de linealidad, ya que en esta gráfica no se obsera ningún patrón.

Ahora analicemos de nuevo el supuesto de homocedasticidad

plot(fit2, 3)

Vemos que con este nuevo modelo tambien parece que se puede asumir el supuesto de homocedasticidad, ahora realizaremos la prueba de Breusch-Pagan para confirmarlo:

lmtest::bptest(fit2)

## ## studentized Breusch-Pagan test ## ## data: fit2 ## BP = 4, df = 10, p-value = 1

Esta vez con los p-values de las dos pruebas no se rechaza la hipótesis de homocedasticidad, por lo que podemos decir que el nuevo modelo cumple con el supuestp de homocedasticidad.

Comprobemos ahora el supuesto de normalidad, construyendo el Q-Q plot:

plot(fit2,2)

Con este gráfico parece cumplirse el supuesto, ahora verifiquemoslo con las pruebas Shapiro-Wilk y Kolmogorov-Smirnov; comencemos con la prueba de Shapiro-Wilk:

Datosfit2=broom::augment(fit2) shapiro.test(Datosfit2$.std.resid)

## ## Shapiro-Wilk normality test ## ## data: Datosfit2$.std.resid ## W = 1, p-value = 0.1

Ahora veamos la prueba de Shapiro-Wilk:

nortest::lillie.test(Datosfit2$.std.resid)

## ## Lilliefors (Kolmogorov-Smirnov) normality test ## ## data: Datosfit2$.std.resid ## D = 0.03, p-value = 0.6

Vemos que en las dos pruebas no se rechaza la hipotesis de que los datos provienen de una distribucion normal, por lo tanto podemos asumir el supuesto de normalidad.

Por ultimo verifiquemos de nuevo el supuesto de aleatoriedad, para esto construimos de nuevo el autocorrelograma

acf(Datosfit2$.std.resid, main = 'Autocorrelograma')

Con base en el autocorrelograma parece que se cumple el supuesto, ahora corroboremos con la prueba de aleatoriedad:

randtests::runs.test(Datosfit2$.std.resid)

## ## Runs Test ## ## data: Datosfit2$.std.resid ## statistic = -0.5, runs = 196, n1 = 200, n2 = 200, n = 400, p-value = ## 0.6 ## alternative hypothesis: nonrandomness

Con este p-value no se rechaza la la hipotesis de aleatoriedad de los datos, así que podemos asumir que se cumple el supuesto de aleatoriedad.

Por lo tanto, podemos decir que el modelo:

\[log(Y) = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \beta_{3}X_{3} + \beta_{4}X_{4} + \beta_{5}X_{5} + \beta_{6}X_{6} + \epsilon\]

Cumple todos los supuestos del modelo de regresion lineal múltiple

Ya que el conjunto de covarables es pequeño, usraemos el ,etodo del mejor subconjunto, además usaremos el o. Adem´as, usaremos el criterio de optimización del estadístico $R^{2}_{adj}$.

subconjuntos <- leaps::regsubsets(logY ~ X1 + X2 + X3 + Z4 + X5 + X6, data = datos, nbest = 2) summary(subconjuntos)

## Subset selection object ## Call: regsubsets.formula(logY ~ X1 + X2 + X3 + Z4 + X5 + X6, data = datos, ## nbest = 2) ## 10 Variables (and intercept) ## Forced in Forced out ## X1A2 FALSE FALSE ## X1A3 FALSE FALSE ## X1A4 FALSE FALSE ## X2B2 FALSE FALSE ## X2B3 FALSE FALSE ## X2B4 FALSE FALSE ## X3 FALSE FALSE ## Z4 FALSE FALSE ## X5 FALSE FALSE ## X6 FALSE FALSE ## 2 subsets of each size up to 8 ## Selection Algorithm: exhaustive ## X1A2 X1A3 X1A4 X2B2 X2B3 X2B4 X3 Z4 X5 X6 ## 1 ( 1 ) " " " " " " " " " " " " " " "*" " " " " ## 1 ( 2 ) "*" " " " " " " " " " " " " " " " " " " ## 2 ( 1 ) "*" " " " " " " " " " " " " "*" " " " " ## 2 ( 2 ) " " " " " " " " " " " " " " "*" " " "*" ## 3 ( 1 ) "*" " " " " " " " " " " " " "*" " " "*" ## 3 ( 2 ) "*" " " "*" " " " " " " " " "*" " " " " ## 4 ( 1 ) "*" "*" "*" " " " " " " " " "*" " " " " ## 4 ( 2 ) "*" " " "*" " " " " " " " " "*" " " "*" ## 5 ( 1 ) "*" "*" "*" " " " " " " " " "*" " " "*" ## 5 ( 2 ) "*" "*" "*" "*" " " " " " " "*" " " " " ## 6 ( 1 ) "*" "*" "*" " " "*" " " " " "*" " " "*" ## 6 ( 2 ) "*" "*" "*" "*" " " " " " " "*" " " "*" ## 7 ( 1 ) "*" "*" "*" "*" "*" " " " " "*" " " "*" ## 7 ( 2 ) "*" "*" "*" " " "*" "*" " " "*" " " "*" ## 8 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" " " "*" ## 8 ( 2 ) "*" "*" "*" "*" "*" " " "*" "*" " " "*"

Una vez hecho esto podemos ver que la variable $Z_{4} = X^{2}_{4}$ está presemte en casi todos los modelos mientras que las variables $X_{3} \hspace{.1cm} y \hspace{.1cm} X_{5}$ no.

Para ayudarnos en encontrar entre estos modelos, al que tiene el estadístico $R^{2}_{adj}$ más alto, construimos la siguiente tabla:

subconjuntos2 <- summary(subconjuntos) combine <- cbind(subconjuntos2$which, subconjuntos2$adjr2) ndim <- dim(subconjuntos2$which) ncols <- ndim[2] dimnames(combine)[[2]][ncols+1] <- "R^2_adj" round(combine, digits = 4)

## (Intercept) X1A2 X1A3 X1A4 X2B2 X2B3 X2B4 X3 Z4 X5 X6 R^2_adj ## 1 1 0 0 0 0 0 0 0 1 0 0 0.67 ## 1 1 1 0 0 0 0 0 0 0 0 0 0.16 ## 2 1 1 0 0 0 0 0 0 1 0 0 0.82 ## 2 1 0 0 0 0 0 0 0 1 0 1 0.76 ## 3 1 1 0 0 0 0 0 0 1 0 1 0.89 ## 3 1 1 0 1 0 0 0 0 1 0 0 0.85 ## 4 1 1 1 1 0 0 0 0 1 0 0 0.94 ## 4 1 1 0 1 0 0 0 0 1 0 1 0.91 ## 5 1 1 1 1 0 0 0 0 1 0 1 0.99 ## 5 1 1 1 1 1 0 0 0 1 0 0 0.94 ## 6 1 1 1 1 0 1 0 0 1 0 1 0.99 ## 6 1 1 1 1 1 0 0 0 1 0 1 0.99 ## 7 1 1 1 1 1 1 0 0 1 0 1 1.00 ## 7 1 1 1 1 0 1 1 0 1 0 1 0.99 ## 8 1 1 1 1 1 1 1 0 1 0 1 1.00 ## 8 1 1 1 1 1 1 0 1 1 0 1 1.00

Así, concluimos que de acuerdo al criterio de optimización del estadístico $R^{2}_{adj}$ el mejor modelo es el que contiene a todas las variables excepto $X_{3} \hspace{.1cm} y \hspace{.1cm} X_{5}$.

Por lo tanto con un $R^{2}_{adj} = 0.9984$, seleccionamos al modelo dado por:

\[log(Y) = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \beta_{4}X_{4} + \beta_{6}X_{6} + \epsilon\]

Examen 2, versión B. Regresión lineal múltiple

19 de diciembre de 2021

1. (1.5 puntos)

2. (1.5 puntos)

3. (2 puntos)

Solución:

4. (2 puntos)

5. (2 puntos)

6. (2 puntos)

Solución: