Aclaración

Cabe resaltar que el uso de este material de estudio es didactico, y los datos usados no representan ningun estudio real, se busca discreción, ya que el objetivo de este trabajo es mostrar el uso de r para la regresion lineal multiple.

Ejemplo 1

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 observadasdurante un período de estudio de 2 semanas:

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

ventas	11
empaques	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.

I.

I. Presente un gráfico para describir los datos, por ejemplo, un boxplot por tipo de empaque.

Solución:

#---- librerías ----
library(ggplot2)
library(dplyr)
library(multcomp)
#----- Datos ----
ventas = c(27,33,23,26,28,11,17,16,14,15,23,20,18,17,12,10,15,18,11)
empaques = c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,4,4,4,4,4)

empaques = factor(empaques)

datos = data.frame(ventas, empaques)

#---- Grafica ----
datos %>% ggplot(aes(x = empaques,
                     y = ventas,
                     col = empaques)) +
  geom_boxplot()+
  theme_bw()+
  labs(title = "Boxplot de los datos de la compañia kenton Food",
       x = "Tipo de empaque",
       y = "Numero de ventas",
       col = "")

El empaque tipo 1 es el que tiene una media más grande que los demás, mientras que el tipo 4 tiene la media de ventas menor al resto; en los empaques tipo 1 y 2 podemos observar que tienen una varianza similar, pero el tipo 1 tiene un outlier.

II.

II. 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.

Solución:

Observemos que el modelo del problema es el siguiente \[ E(Y;X) = \beta_0 + \beta_1x_1 + \beta_2x_2+\beta_3x_3 \] donde \(x_i\) es 1 si es el empaque \(i+1\) y 0 en otro caso.

Ahora veamos caso por caso \[ E(Y;X = empaque1 ) = \beta_0 \] \[ E(Y;X = empaque2 ) = \beta_0 + \beta_1 \] \[ E(Y;X = empaque3 ) = \beta_0 + \beta_2 \] \[ E(Y;X = empaque4) = \beta_0 + \beta_3 \]

#---- Modelo de regresión multiple ----

fit = lm(ventas ~ empaques,data = datos)
coeficientes = coef(fit)
medias = c(coeficientes[1],coeficientes[1] + coeficientes[2:4])

names(coeficientes) = c("B_0",
                       "B_1",
                       "B_2",
                       "B_3")

names(medias) = c("E(y;x = empaque1)",
                  "E(y;x = empaque2)",
                  "E(y;x = empaque3)",
                  "E(y;x = empaque4)")
medias

## E(y;x = empaque1) E(y;x = empaque2) E(y;x = empaque3) E(y;x = empaque4) 
##              27.4              14.6              19.5              13.2

coeficientes

##   B_0   B_1   B_2   B_3 
##  27.4 -12.8  -7.9 -14.2

III.

III. Escriba las hipótesis que se contrastan con la tabla ANOVA, calcule ésta e interprete. Use \(\alpha = 0.05\). Solución:

\[ H_0: \beta_1 = 0 \wedge \beta_2 = 0 \wedge \beta_3 = 0 \thinspace vs \thinspace H_a : \beta_i \neq 0 \thinspace para \thinspace algun \thinspace i = {1,2,3} \]

summary(fit)

## 
## Call:
## lm(formula = ventas ~ empaques, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -4.40  -1.85  -0.40   1.60   5.60 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   27.400      1.357  20.185 2.76e-12 ***
## empaques2    -12.800      1.920  -6.668 7.51e-06 ***
## empaques3     -7.900      2.036  -3.880  0.00148 ** 
## empaques4    -14.200      1.920  -7.397 2.23e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.035 on 15 degrees of freedom
## Multiple R-squared:  0.8168, Adjusted R-squared:  0.7802 
## F-statistic: 22.29 on 3 and 15 DF,  p-value: 8.794e-06

Nos fijamos p-value: 8.794e-06 < 0.05 por lo que es plausible rechazar \(H_0\) y concluir el tipo de empaque afecta al promedio de ventas.

IV.

IV. ¿Se puede considerar que el diseño del empaque afecta las ventas promedio? Use \(\alpha= 0.05.\) Argumente indicando con claridad qué hipótesis se están contrastando en términos de los parámetros del modelo ajustado.

Solución:

En este caso vamos a utilizar otro código ya que es la misma prueba pero ahora utilizaremos la librería de R llamada multcomp

Donde la hipótesis es la siguiente:

\[ K = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix} \]

Esta matriz K cada renglón representa la combinación lineal de las \(\beta\)´s del modelo.

\[ m = \begin{pmatrix} 0\\ 0\\ 0\\ \end{pmatrix} \]

Ahora este vector representa la igualdad de cada una de las combinaciones lineales de K, ya con estas dos estructuras podemos hacer la siguiente prueba de hipótesis:

\[ H_0: K\beta = m\thinspace vs \thinspace H_a: K \beta \neq m\]

k = matrix(c(0,1,0,0,
             0,0,1,0,
             0,0,0,1), nrow = 3, byrow = T)
m = c(0,0,0)

summary(glht(fit,linfct = k,rhs = m), test = Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    -12.8
## 2 == 0     -7.9
## 3 == 0    -14.2
## 
## Global Test:
##       F DF1 DF2    Pr(>F)
## 1 22.29   3  15 8.794e-06

Notemos que 8.794e-06 < 0.05 entonces podemos concluir lo mismo que el inciso anterior.

V.

V. Realice 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 = 0.05.\) Interprete los resultados.

Solución:

La hipótesis tiene que partir de: \[ Ventas\thinspace i := Ventas \thinspace promedio \thinspace del \thinspace empaque \thinspace i \] \[ Ventas\thinspace 1 = Ventas \thinspace2 \] \[ Ventas\thinspace 1 = Ventas \thinspace3 \] \[ Ventas\thinspace 1 = Ventas \thinspace4 \] \[ Ventas\thinspace 2 = Ventas \thinspace3 \] \[ Ventas\thinspace 2 = Ventas \thinspace4 \] \[ Ventas\thinspace 3 = Ventas \thinspace4 \] Que en nuestro modelo se puede traducir como: \[ \beta_0 = \beta_0+\beta_1 \Leftrightarrow \beta_1=0 \] \[ \beta_0 = \beta_0+\beta_2 \Leftrightarrow \beta_2=0\] \[ \beta_0 = \beta_0+\beta_3 \Leftrightarrow \beta_3 = 0\] \[ \beta_0 + \beta_1 = \beta_0+\beta_2 \Leftrightarrow \beta_1 - \beta_2 =0\] \[ \beta_0 + \beta_1 = \beta_0+\beta_3 \Leftrightarrow \beta_1 -\beta_3 = 0 \] \[ \beta_0 + \beta_2 = \beta_0+\beta_3 \Leftrightarrow \beta_2-\beta_3 = 0 \]

Asi que la matriz K de las combinaciones lineales es la siguiente:

\[ K = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \\ 0 & 1 & -1 & 0\\ 0 & 1 & 0 &-1 \\ 0 &0&1&-1 \end{pmatrix} \]

Mientras que el vector en el que se busca la igualdad es el siguiente:

\[ m = \begin{pmatrix} 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \end{pmatrix} \] Entonces la prueba de hipótesis es la siguiente:

\[ H_0 : K\beta = m \thinspace vs \thinspace H_a: K\beta \neq m \]

k = 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)
k = matrix(data = k,nrow = 6,ncol = 4,byrow = T)

m = rep(0,6)

summary(glht(fit,linfct = k,rhs = m))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = ventas ~ empaques, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  -12.800      1.920  -6.668  < 0.001 ***
## 2 == 0   -7.900      2.036  -3.880  0.00736 ** 
## 3 == 0  -14.200      1.920  -7.397  < 0.001 ***
## 4 == 0   -4.900      2.036  -2.406  0.11848    
## 5 == 0    1.400      1.920   0.729  0.88373    
## 6 == 0    6.300      2.036   3.094  0.03334 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

De la prueba simultánea podemos concluir que es muy plausible decir que el empaque 1 tiene ventas diferentes a los demás, mientras que las ventas promedio del empaque 2 y 3, y el empaque 2 y 4 según la prueba no se rechaza que sean iguales, pero como vimos en el boxplot las ventas promedio del empaque 2 esta entre las ventas promedio del empaque 3 y 4, por lo que con una signicancia del 0.05 podemos decir que las ventas de 3 y 4 no son iguales.

VI.

VI. Suponga que los ejecutivos de la empresa tienen la sospecha de que el diseño de empaque 1 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= 0.05\)

Solución:

Contruimos la hipótesis alternativa (la que queremos comprobar), utilizando la notación del inciso pasado:

\[ Ventas \thinspace 1> Ventas \thinspace2 \] \[ Ventas \thinspace 1> Ventas \thinspace3\] \[ Ventas \thinspace 1> Ventas \thinspace4 \] Y en el modelo tenemos que:

\[ \beta_0 > \beta_0 + \beta_1 \Leftrightarrow \beta_1<0 \] \[ \beta_0 > \beta_0 + \beta_2 \Leftrightarrow \beta_2<0 \] \[ \beta_0 > \beta_0 + \beta_3 \Leftrightarrow \beta_3<0 \] Por lo que en la hipótesis alternativa

\[ H_a: \beta_i < 0 \thinspace para \thinspace algun \thinspace i = 1,2,3 \] Mientras que la hipótesis nula global \[ H_0: \beta_i \geq 0 \thinspace \forall i \in \{ 1,2,3 \} \] Contruimos la matriz K para las combinaciones lineales y el vector m para la comparación de las combinaciones

\[ K = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \\ \end{pmatrix} \]

\[ m = \begin{pmatrix} 0\\ 0\\ 0\\ \end{pmatrix} \]

k= c(0,1,0,0,
     0,0,1,0,
     0,0,0,1)

k = matrix(data = k, nrow = 3,ncol = 4,byrow = TRUE)

m = c(0,0,0)

summary(glht(fit, linfct=k, rhs=m, alternative="less"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = ventas ~ empaques, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(<t)    
## 1 >= 0  -12.800      1.920  -6.668 < 0.001 ***
## 2 >= 0   -7.900      2.036  -3.880 0.00198 ** 
## 3 >= 0  -14.200      1.920  -7.397 < 0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Asi que con una significancia de 0.05 notemos que rechazamos la hipótesis nula global, asi que es plausible decir que el empaque 1 en promedio tuvo mas ventas que los demás empaques.

Ejemplo 2

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 individuos 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 Ex4A.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)

I.

I. Realice un análisis descriptivo de los datos. Dado que las dos covariables son categóricas, incluya 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:

#---- Preambulo ----
rm(list = ls(all.names = T))
setwd("~/Facultad/semestre_7/Modelos no Parametricos/TareaExamen2")

#librerías
library(ggplot2)
library(dplyr)
library(multcomp)
library(readr)

#---- importacion de los datos ----

datos = read_csv("datos/Ex4A.csv")
datos$Sexo = factor(datos$Sexo)
datos$Trat = factor(datos$Trat)

#---- grafica boxplot ----

datos %>% mutate(nueva = paste(Sexo, Trat, sep = "|")) %>% ggplot(aes(x = nueva,
                                                                      y = Puntaje,
                                                                      col = Sexo))+
  geom_boxplot() +
  theme_bw() +
  labs(title ="Grafico de caja",
       subtitle = "Puntaje según Sexo y tratamiento",
      x = "Sexo|Tratamiento",
      y= "Puntaje de ansiedad",
      caption = "el valor control indica que no tiene tratamiento\nmientras que Trat1 hace referencia al tratamiento 1 y analogamente para Trat2") + 
  theme(axis.text.x = element_text(angle = 90))

En la gráfica podemos observar que en efecto hay diferencia al tratamiento de acuerdo al sexo, el tratamiento 1 es mas efectivo en mujeres, mientras que el 2 le funcionó a los hombres, y el control es similar en ambos sexos.

II.

II. Considerando un modelo de regresión que incluye las dos variables categóricas de forma individual y también su interacción, dé la expresión del puntaje promedio para cada valor de las varia bles categóricas, es decir: E(puntaje; T rat = k, Sexo = l), con k ∈ {Control, T rat1, T rat2} y l ∈ {Hombre, Mujer}; así como la estimación puntual correspondiente.

Solución:

\[ E(Y;X) = \beta_0 + \beta_1x_1+\beta_2x_2+ \beta_3x_3 + \beta_4x_1x_2 + \beta_5x_1x_3\] donde \(x_1\) es una indicadora si es mujer, mientras que \(x_2\) y \(x_3\) son indicadoras del tratamiento 1 o el tratamiento 2 respectivamente.

Para Hombres:

\[ E(y;x = (Hombre, no\thinspace tratamiento ) ) = \beta_0 \] \[ E(y;x = (Hombre, tratamiento \thinspace 1) ) = \beta_0 + \beta_2 \] \[ E(y;x = (Hombre, tratamiento 2) ) = \beta_0 + \beta_3 \] Para Mujeres: \[ E(y;x = (Mujer, no tratamiento ) ) = \beta_0 + \beta_1 \] \[ E(y;x = (Mujer, tratamiento 1) ) = \beta_0 + \beta_1 + \beta_2 + \beta_4 \] \[ E(y;x = (Mujer, tratamiento 2) ) = \beta_0 + \beta_1 + \beta_3 + \beta_5\]

#---- modelo de regresión ----

fit = lm(Puntaje ~ Sexo + Trat+Trat*Sexo, data = datos)

coeficientes = coef(fit)

esperanza = c(coeficientes[1],#Hombre No Tratamiento
              coeficientes[1] + coeficientes[3:4],#Hombre Tratamiento 1 y Tratamiento 2
              coeficientes[1] + coeficientes[2], #Mujer No tramiento 
              coeficientes[1] + coeficientes[2]+ coeficientes[3:4]+coeficientes[5:6]) #Mujer tratamientos

names(coeficientes) = c("b0","b1","b2","b3", "b4", "b5")

names(esperanza) = c("E(y;x = (Hombre, no tratamiento ))",
                     "E(y;x = (Hombre, tratamiento  1))",
                     "E(y;x = (Hombre, tratamiento  2))",
                    "E(y;x = (Mujer, no tratamiento ))",
                    "E(y;x = (Mujer, tratamiento  1))",
                    "E(y;x = (Mujer, tratamiento  2))")

coeficientes

##         b0         b1         b2         b3         b4         b5 
## 10.2832476 -0.5230817 -2.8857619 -4.0702771  1.3757855  3.5913503

esperanza

## E(y;x = (Hombre, no tratamiento ))  E(y;x = (Hombre, tratamiento  1)) 
##                          10.283248                           7.397486 
##  E(y;x = (Hombre, tratamiento  2))  E(y;x = (Mujer, no tratamiento )) 
##                           6.212970                           9.760166 
##   E(y;x = (Mujer, tratamiento  1))   E(y;x = (Mujer, tratamiento  2)) 
##                           8.250189                           9.281239

III.

III. Escriba las hipótesis que se contrastan con la tabla ANOVA, calcule ésta e interprete. Use \(\alpha = 0.05\).

Solución:

\[ H_0: \beta_1 = 0 \wedge \beta_2 = 0 \wedge \beta_3 = 0 \wedge\beta_4 = 0 \wedge \beta_5 = 0\thinspace vs \thinspace H_a : \beta_i \neq 0 \thinspace para \thinspace algun \thinspace i = {1,2,3,4,5} \]

summary(fit)

## 
## Call:
## lm(formula = Puntaje ~ Sexo + Trat + Trat * Sexo, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3785 -1.1800 -0.0518  1.2159  4.3400 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          10.2832     0.3948  26.044  < 2e-16 ***
## SexoMujer            -0.5231     0.5584  -0.937   0.3509    
## TratTrat1            -2.8858     0.5584  -5.168 1.02e-06 ***
## TratTrat2            -4.0703     0.5584  -7.289 4.39e-11 ***
## SexoMujer:TratTrat1   1.3758     0.7897   1.742   0.0842 .  
## SexoMujer:TratTrat2   3.5914     0.7897   4.548 1.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.766 on 114 degrees of freedom
## Multiple R-squared:  0.4007, Adjusted R-squared:  0.3744 
## F-statistic: 15.24 on 5 and 114 DF,  p-value: 1.873e-11

Notemos que p-value: 1.873e-11 < 0.05 por lo que es plausible rechazar \(H_0\) por lo que alguna de las betas aportan algo al promedio.

IV.

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 = 0.025\) Interprete.

Solución:

Haciendo un análisis tenemos que hacer una prueba de hipótesis tal que \[ E(Y;X = (Hombre, tratamiento))=\beta_0+\beta_2x_2+\beta_3x_3 =\] \[E(Y; X=(Mujer, tratamiento)) = \beta_0+\beta_1+ \beta_2x_2+\beta_3x_3 + \beta_4x_2 + \beta_5x_5\] Por lo que la prueba tiene que verse como \(\beta_1 = 0\), \(\beta_4 = 0\) y \(\beta_5 = 0\) es decir:

\[ k = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ \end{pmatrix} \]

\[ m= \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \]

\[ H_0 : k\beta = m \thinspace \thinspace vs \thinspace \thinspace H_a:k\beta \neq m \] Para esta prueba utilizamos la librería multcomp de R.

k = c(0,1,0,0,0,0,
      0,0,0,0,1,0,
      0,0,0,0,0,1)

  k = matrix(data = k, nrow = 3,ncol = 6, byrow = T)
 m = c(0, 0 ,0)
 summary(glht(fit, linfct = k , rhs = m), test = Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0  -0.5231
## 2 == 0   1.3758
## 3 == 0   3.5914
## 
## Global Test:
##       F DF1 DF2    Pr(>F)
## 1 11.13   3 114 1.828e-06

Notemos que Pr(>F) = 1.828e-06 < 0.025 por lo que es plausible rechazar la hipótesis nula, es decir que el sexo tiene un efecto en la puntuación de los niveles de ansiedad con significancia 0.025.

Reducción del modelo.

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 reduzca el modelo si es posible.

Hagamos una prueba simultánea para reconocer donde se diferencia:

Igualamos los promedios de los tratamientos:

Que el promedio de ser Hombre y no recibir tratamiento sea igual al promedio de ser Mujer y no recibir tratamiento:

\[ \beta_0 = \beta_0 + \beta_1 \Leftrightarrow \beta_1 = 0\]

Que el promedio de ser Hombre y recibir el tratamiento actual sea igual al promedio de ser Mujer y recibir el tratamiento actual:

\[ \beta_0 + \beta_2 = \beta_0 + \beta_1 + \beta_2 + \beta_4 \Leftrightarrow \beta_1 + \beta_4 = 0\] Que el promedio de ser Hombre y recibir el tratamiento nuevo sea igual al promedio de ser Mujer y recibir el tratamiento nuevo:

\[ \beta_0 + \beta_3 = \beta_0+\beta_1 + \beta_3 + \beta_5 \Leftrightarrow \beta_1 +\beta_5 = 0 \] Ya con esto armamos la matriz k para las combinaciones lineales:

\[ k = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 1& 0 & 0 & 1 & 0\\ 0 & 1 & 0 & 0 & 0 & 1\\ \end{pmatrix} \]

mientras que

\[ m= \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \]

\[ H_0 : k\beta = m \thinspace vs \thinspace H_a: k\beta \neq m \]

k = c(0,1,0,0,0,0,
      0,1,0,0,1,0,
      0,1,0,0,0,1)

k = matrix(data = k, nrow = 3,ncol = 6, byrow = T)
m = c(0, 0 ,0)
 summary(glht(fit, linfct = k , rhs = m))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Puntaje ~ Sexo + Trat + Trat * Sexo, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  -0.5231     0.5584  -0.937    0.724    
## 2 == 0   0.8527     0.5584   1.527    0.339    
## 3 == 0   3.0683     0.5584   5.495   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Ya con esta prueba veamos que es plausible decir que los promedios entre hombres y mujeres al no aplicar tratamiento y al aplicar el tratamiento 1 no hay distinción por sexo.

Asi que modificamos nuesto dataset para poder crear un modelo que nos ayude a modelar mejor este problema:

datos2 = datos %>% mutate(x1 = case_when(
  Sexo == "Hombre" ~ 0,
  Sexo == "Mujer" ~ 1
), x2 = case_when(
  Trat == "Control" ~ 0,
  Trat == "Trat1" ~ 1,
  Trat == "Trat2" ~0
), x3 = case_when(
  Trat == "Control"~ 0,
  Trat == "Trat1" ~ 0,
  Trat == "Trat2" ~ 1)) %>% mutate( x1x2 = x1*x2, x1x3 = x1*x3) %>% mutate(x1 = factor(x1), x2 = factor(x2), x3 = factor(x3), x1x2 = factor(x1x2),x1x3 = factor(x1x3))

fit3 = lm(Puntaje ~ x2 + x3 + x1x3, data = datos2)
summary(fit3)

## 
## Call:
## lm(formula = Puntaje ~ x2 + x3 + x1x3, data = datos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6400 -1.0701  0.0033  1.0982  4.1249 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.0217     0.2806  35.709  < 2e-16 ***
## x21          -2.1979     0.3969  -5.538 1.94e-07 ***
## x31          -3.8087     0.4861  -7.835 2.47e-12 ***
## x1x31         3.0683     0.5613   5.466 2.67e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.775 on 116 degrees of freedom
## Multiple R-squared:  0.3838, Adjusted R-squared:  0.3678 
## F-statistic: 24.08 on 3 and 116 DF,  p-value: 3.461e-12

V

V. 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; T rat = k, Sexo = l), con k ∈ {Control, T rat1, T rat2} y l ∈ {Hombre, Mujer}; así como estimaciones puntuales.

Solución:

\[ E(Y; X =(Tratamiento, Sexo)) = \beta_0 + \beta_2x_2 +\beta_3x_3 +\beta_5x_1x_3\] \[ E(Y;X = (Hombre, no \thinspace tratamiento)) = \beta_0 \] \[ E(Y; X = (Hombre, tratamiento \thinspace 1))= \beta_0 + \beta_2 \] \[ E(Y; X = (Hombre, tratamiento \thinspace 2))= \beta_0 + \beta_3 \] \[ E(Y; X = (Mujer, no \thinspace tratamiento))= \beta_0 \] \[ E(Y; X = (Mujer, tratamiento \thinspace 1))= \beta_0 + \beta_2 \] \[ E(Y; X = (Mujer, tratamiento \thinspace 2))= \beta_0 + \beta_3 + \beta_5 \]

Notemos que en este nuevo modelo el sexo no tiene relevancia en el promedio del puntaje de ansiedad en los grupos que no se les aplicó y tampoco en la poblacion del tratamiento 1; Que se puede ver como las siguientes expresiónes:

\[ E(No\thinspace tratamiento) = \beta_0\] \[ E(Tratamiento\thinspace1) = \beta_0 + \beta_2\] \[ E(Y; X = (Hombre, tratamiento \thinspace 2))= \beta_0 + \beta_3 \] \[ E(Y; X = (Mujer, tratamiento \thinspace 2))= \beta_0 + \beta_3 + \beta_5 \]

coeficientes = coef(fit3)

names(coeficientes) = c("b0", "b2", "b3", "b5")

Esperanzas = c(coeficientes[1],
               coeficientes[1] + coeficientes[2],
               coeficientes[1] + coeficientes[3],
               coeficientes[1] + coeficientes[3] + coeficientes[4])

names(Esperanzas) = c("E(no tratamiento)",
                      "E(Tratamiento 1)",
                      "E(Hombre, Tratamiento 2)",
                      "E(Mujer, Tratamiento 2)")

Esperanzas

##        E(no tratamiento)         E(Tratamiento 1) E(Hombre, Tratamiento 2) 
##                10.021707                 7.823838                 6.212970 
##  E(Mujer, Tratamiento 2) 
##                 9.281239

VI

VI. 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:

Usando las expresiónes anteriores podemos sacar las siguientes combinaciones lineales de los parámetros para la prueba que se nos pide:

El promedio del puntaje de los Hombres con el nuevo tratamiento es menor al promedio del puntaje de los hombres sin ningun tratamiento: \[ \beta_0 + \beta_3 < \beta_0 \Leftrightarrow \beta_3 < 0 \] El promedio del puntaje de los Hombres con el nuevo tratamiento es menor al promedio de puntaje de los Hombres que se les aplicó el tratamiento actual: \[ \beta_0 + \beta_3 < \beta_0 +\beta_2 \Leftrightarrow \beta_3 - \beta_2 < 0 \] El promedio del puntaje de las mujeres con el nuevo tratamiento es menor al promedio del puntaje de las mujeres que no se le aplicó ningun tratamiento. \[\beta_0 + \beta_3 + \beta_5 < \beta_0 \Leftrightarrow \beta_3+\beta_5 < 0 \] El promedio del puntaje de las mujeres con el nuevo tratamiento es menor al promedio del puntaje de las mujeres que se le aplicó el tratamiento actual: \[ \beta_0 + \beta_3 +\beta_5 < \beta_0 + \beta_2 \Leftrightarrow \beta_3 +\beta_5 - \beta_2<0 \] Viendolo con la matriz k (matriz de combinaciones lineales de la prueba): \[ k = \begin{pmatrix} 0 & 0 &1 & 0 \\ 0 & -1 &1 & 0 \\ 0 & 0 &1 & 1 \\ 0 & -1 &1 & 1 \end{pmatrix} \]

Mientras que:

\[ m= \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \]

Las hipótesis alternativas son : \[ \beta_3 < 0 \] \[ \beta_3 - \beta_2 < 0 \] \[\beta_3+\beta_5 < 0 \] \[ \beta_3 +\beta_5 - \beta_2<0 \] Y, las nulas son las siguientes: \[ \beta_3 \geq 0 \] \[ \beta_3 - \beta_2 \geq 0 \] \[ \beta_3+\beta_5\geq0 \] \[ \beta_3 +\beta_5 - \beta_2\geq0 \]

k = c(0,0,1,0,
      0,-1,1,0,
      0,0,1,1,
      0,-1,1,1)
(k = matrix(data = k, nrow = 4, ncol = 4, byrow = T))

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    1    0
## [2,]    0   -1    1    0
## [3,]    0    0    1    1
## [4,]    0   -1    1    1

m = c(0,0,0,0)

summary(glht(fit3, linfct = k, rhs = m, alternative = "less"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Puntaje ~ x2 + x3 + x1x3, data = datos2)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(<t)    
## 1 >= 0  -3.8087     0.4861  -7.835 < 0.001 ***
## 2 >= 0  -1.6109     0.4861  -3.314 0.00239 ** 
## 3 >= 0  -0.7405     0.4861  -1.523 0.18733    
## 4 >= 0   1.4574     0.4861   2.998 1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

De acuerdo con lo visto en la prueba simultánea, es plausible que el tratamiento nuevo sea mejor que el actual y que no utilizar un tratamiento, pero solo para hombres, mientras que para mujeres no tenemos evidencia para concluir que es el mejor.

VII.

VII. 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 en hombres, aunque el tratamiento actual lo tiene en mujeres. Use \(\alpha = .05\):

Solución:

Primero fijémonos en las combinaciones lineales

Nuevo tratamiento en los hombres es menor a no tener tratamiento siendo hombre en promedio

\[ \beta_0 + \beta_3 < \beta_0 \Leftrightarrow \beta_3 < 0\] Nuevo tratammiento en los hombre es menor al aplicar el tratamiento actual en hombres en promedio: \[ \beta_0 + \beta_3 < \beta_0 + \beta_2 \Leftrightarrow \beta_3 - \beta_2 < 0 \] El tratamiento actual es menor en mujeres que el no aplicar tratamiento en mujeres en promedio: \[ \beta_0 + \beta_2 < \beta_0 \Leftrightarrow \beta_2 < 0\] El tratamiento actual es menor en mujeres que el tratamiento nuevo en mujeres en promedio: \[\beta_0 + \beta_2 < \beta_0 + \beta_3 + \beta_5 \Leftrightarrow \beta_2 - \beta_3 - \beta_5 < 0\] Ya con esto podemos crear la matriz k de las combinaciones lineales de los \(\beta's\)

\[ k = \begin{pmatrix} 0 & 0 &1 & 0 \\ 0 & -1 &1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 &-1 & -1 \end{pmatrix} \]

Por lo que la m es la siguiente:

\[ m= \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \]

Nuestras hipótesis alternativas son las siguientes:

\[ \beta_3 < 0 \] \[ \beta_3 - \beta_2 < 0 \] \[ \beta_2 < 0 \] \[ \beta_2 - \beta_3 - \beta_5 < 0 \] Por lo que, nuestras hiptesis nulas son las siguientes

\[ \beta_3 \geq 0 \] \[ \beta_3 - \beta_2 \geq0 \] \[ \beta_2 \geq 0 \] \[ \beta_2 - \beta_3 - \beta_5 \geq 0 \]

k = c(0,0,1,0,
      0,-1,1,0,
      0,1,0,0,
      0,1,-1,-1)
(k = matrix(data = k, nrow = 4, ncol = 4, byrow = T))

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    1    0
## [2,]    0   -1    1    0
## [3,]    0    1    0    0
## [4,]    0    1   -1   -1

m = c(0,0,0,0)

summary(glht(fit3, linfct = k, rhs = m, alternative = "less"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Puntaje ~ x2 + x3 + x1x3, data = datos2)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(<t)    
## 1 >= 0  -3.8087     0.4861  -7.835 < 0.001 ***
## 2 >= 0  -1.6109     0.4861  -3.314 0.00230 ** 
## 3 >= 0  -2.1979     0.3969  -5.538 < 0.001 ***
## 4 >= 0  -1.4574     0.4861  -2.998 0.00632 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Con la prueba podemos concluir que en hombres el tratamiento nuevo es mejor, mientras que en mujeres el actual tiene mejor desempeño.

Ejemplo 3

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)

I.

I. 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:

# ---- previa ----
rm(list = ls(all.names = T))
setwd("~/Facultad/semestre_7/Modelos no Parametricos/TareaExamen2")

#librerías
library(ggplot2)
library(dplyr)
library(multcomp)
library(readr)

# ---- importamos los datos ----

datos = read_csv("datos/Ex5.csv")
datos = datos %>% mutate(Trat = factor(Trat))

# ---- análisis de los datos ----

summary(datos)

##       Ant              Trat          Edad      
##  Min.   : 1.851   Control:200   Min.   :16.00  
##  1st Qu.:11.569   Med    :100   1st Qu.:27.00  
##  Median :14.682                 Median :38.00  
##  Mean   :14.587                 Mean   :37.67  
##  3rd Qu.:17.959                 3rd Qu.:49.00  
##  Max.   :25.075                 Max.   :60.00

summary(datos %>% filter(Trat == "Control"))

##       Ant              Trat          Edad      
##  Min.   : 1.851   Control:200   Min.   :16.00  
##  1st Qu.:10.369   Med    :  0   1st Qu.:27.00  
##  Median :13.280                 Median :39.00  
##  Mean   :13.629                 Mean   :37.87  
##  3rd Qu.:16.741                 3rd Qu.:49.00  
##  Max.   :25.075                 Max.   :60.00

summary(datos %>% filter(Trat == "Med"))

##       Ant              Trat          Edad      
##  Min.   : 8.987   Control:  0   Min.   :16.00  
##  1st Qu.:14.216   Med    :100   1st Qu.:26.75  
##  Median :16.742                 Median :37.00  
##  Mean   :16.503                 Mean   :37.26  
##  3rd Qu.:18.772                 3rd Qu.:50.00  
##  Max.   :23.415                 Max.   :59.00

Podemos notar que hay más personas sin medicamento que personas con tratamiento.
En el summary filtrando la información de los pacientes que tuvieron medicamento podemos ver que la variación de anticuerpos es menor que en el grupo control, esto podria ser un indicio de que hay un efecto en nivel de anticuerpos según el tratamiento.
El intervalo de edades en los dos grupos de personas es casi igual.

datos %>% ggplot(aes(x = Edad, y = Ant, col = Trat)) +
  geom_point() +
  theme_bw() +
  labs(title = "Diagrama de dispersion",
       subtitle = "Anticuerpos según edad y su tratamiento",
       y = "Nivel de anticuerpos",
       col = "Tipo de tratamiento",
       caption = "El tramiento control indica que no se le aplicó el medicamento \n mientras que Med indica que se le aplicó") + 
  scale_colour_manual(values = c("#090088",
                                 "#E4007C"))

Entre mas edad menos anticuerpos.
Es mas notable la diferencia de anticuerpos entre los dos grupos, conforme se aumenta la edad.

II.

II. 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 estas dos.

Solución:

# ---- modelo  de regesion lineal ----

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

summary(fit)

## 
## Call:
## lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6211 -1.9539  0.0277  1.6018  7.0063 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.34298    0.57573  42.282  < 2e-16 ***
## TratMed      -2.25730    0.96763  -2.333   0.0203 *  
## Edad         -0.28290    0.01440 -19.645  < 2e-16 ***
## TratMed:Edad  0.13307    0.02437   5.460 1.01e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.61 on 296 degrees of freedom
## Multiple R-squared:  0.6394, Adjusted R-squared:  0.6357 
## F-statistic: 174.9 on 3 and 296 DF,  p-value: < 2.2e-16

Donde el modelo es el siguiente \[ E(Y; X =(Tratamiento, Edad)) = \beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2 \] donde \(x_1\) es una indicadora que si el tramiento fue aplicado (con valor de 1 y 0 si no fue aplicado) y \(x_2\) es la edad.

III.

III. De acuerdo con el modelo ajustado, indique las expresiónes 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:

Para los que no se les dio el medicamento:

\[ E(Y; X = (Control, Edad)) = \beta_0 + \beta_2(Edad) \] Para los que se les dio el medicamento:

\[ E(Y; X = (Med, Edad)) = \beta_0 + \beta_1 + \beta_2(Edad) + \beta_3(Edad) \] \[ E(Y; X = (Med, Edad)) = (\beta_0 + \beta_1 )+ (\beta_2 + \beta_3)(Edad) \]

b0 = coef(fit)[1]
b1 = coef(fit)[2]
b2 = coef(fit)[3]
b3 = coef(fit)[4]
control <- function(edad) {
  b0 + b2*(edad)
}
med <- function(edad) {
  b0+b1+(b2+b3)*edad
}
datos %>% ggplot(aes(x = Edad, y = Ant, col = Trat)) +
  geom_point() +
  theme_bw() +
  labs(title = "Diagrama de dispersion",
       subtitle = "Anticueropos según edad y su tratamiento",
       y = "Nivel de anticuerpos",
       col = "Tipo de tratamiento",
       caption = "El tramiento control indica que no se le aplicó el medicamento \n mientras que Med indica que se le aplicó") + 
  scale_colour_manual(values = c("#090088",
                                 "#E4007C")) + geom_function(fun = control, col ="#090088" ) + 
geom_function(fun = med, col = "#E4007C")

Podemos decir que la pendiente es mas marcada en el grupo de las personas que no se les proporcionó medicamento

IV

IV. ¿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:

Utilizando lo del inciso anterior y con lo que nos pide se debe a hacer una prueba de hipótesis sobre \(\beta_3\).

Utilizamos la librería multcomp para esta prueba:

\[ K = \begin{pmatrix} 0 & 0& 0 & 1 \end{pmatrix} \]

\[ m = \begin{pmatrix} 0 \end{pmatrix} \]

\[ H_0: K\beta = m\thinspace vs \thinspace H_a: K \beta \ne m\]

summary(fit)

## 
## Call:
## lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6211 -1.9539  0.0277  1.6018  7.0063 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.34298    0.57573  42.282  < 2e-16 ***
## TratMed      -2.25730    0.96763  -2.333   0.0203 *  
## Edad         -0.28290    0.01440 -19.645  < 2e-16 ***
## TratMed:Edad  0.13307    0.02437   5.460 1.01e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.61 on 296 degrees of freedom
## Multiple R-squared:  0.6394, Adjusted R-squared:  0.6357 
## F-statistic: 174.9 on 3 and 296 DF,  p-value: < 2.2e-16

k = matrix(data = c(0,0,0,1),
           nrow = 1,
           ncol = 4,
           byrow = T)
m = c(0)

summary(glht(fit, linfct = k, rhs = m), test = Ftest())

## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0   0.1331
## 
## Global Test:
##       F DF1 DF2    Pr(>F)
## 1 29.81   1 296 1.008e-07

Notemos que el 1.008e-07 < 0.05 es decir que es plausible rechazar la hipótesis nula, y con un nivel de significancia de 0.05 la edad afecta de diferente manera en las dos poblaciones.

V.

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

Solución:

Como no se ajustó, el modelo sigue siendo igual asi que el análisis sigue siendo el mismo, como tambien la interpretación.

VI.

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

Solución:

Para contestar esta pregunta vamos a hacer 7 pruebas de hipótesis y tambien vamos a verificar intervalos de confianza simultáneos.

Pruebas de hipótesis.

Todas las pruebas van a tener la misma estructura asi que vamos a observar como van a ser las combinaciones lineales

\[ H_{01}: E(Y;X=(Medicamento,\thinspace edad)) \leq E(Y; X= (No \thinspace\thinspace \thinspace tratamiento, edad)) \thinspace \thinspace vs \] \[H_{a1} : E(Y;X=(Medicamento,\thinspace edad)) > E(Y; X= (No \thinspace\thinspace \thinspace tratamiento, edad)) \] Desarrollamos la expresión para \(E(Y;X=(Medicamento,\thinspace edad)) > E(Y; X= (No \thinspace\thinspace \thinspace tratamiento, edad))\):

\[ (\beta_0 + \beta_1) +(\beta_2 + \beta_3)(edad) > \beta_0 + \beta_2(edad) \Rightarrow \beta_1 + \beta_3(edad)>0\] Por lo que la combinación lineal es la siguiente

\[ K = \begin{pmatrix} 0 & 1 & 0 & edad\end{pmatrix}\] \[ m = (0) \]

prueba 1: (edad 30)

k = matrix(c(0,1,0,30),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value   Pr(>t)    
## 1 <= 0    1.735      0.368   4.715 1.87e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 30.

prueba 2: (edad 35)

k = matrix(c(0,1,0,35),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(>t)    
## 1 <= 0   2.4003     0.3254   7.377 8.2e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 35.

prueba 3: (edad 40)

k = matrix(c(0,1,0,40),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0   3.0656     0.3256   9.415 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 40.

prueba 4: (edad 45)

k = matrix(c(0,1,0,45),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0   3.7310     0.3686   10.12 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 45.

prueba 5: (edad 50)

k = matrix(c(0,1,0,50),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0   4.3964     0.4421   9.944 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 50.

prueba 6: (edad 55)

k = matrix(c(0,1,0,55),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0   5.0617     0.5336   9.486 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos con el medicamento en las personas de edad 55.

prueba 7: (edad 60)

k = matrix(c(0,1,0,60),nrow = 1, byrow = T)

m = c(0)

summary(glht(fit, linfct = k, rhs = m, alternative = "greater"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ant ~ Trat + Edad + Trat * Edad, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)    
## 1 <= 0   5.7271     0.6353   9.014 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la hipótesis nula, es plausible decir que se eleva el promedio de nivel de anticuerpos el medicamento en las personas de edad 60.

Intervalos de confianza.

En los intervalos de confianza, la confianza se reducirá al 90%

edades = seq(from = 30, to = 60, by = 0.5)
medicamento = cbind(1,1,edades,edades)
no = cbind(1,0,edades,0)
# los primeros 61 son los de medicamento
K = rbind(medicamento,no)
fitE = glht(fit,linfct = K)
fitci = confint(fitE, level = 0.90)
confianza = data.frame(fitci$confint)
confianza$EDAD = c(edades, edades)
confianza$Tipo = c(replicate(n = 61, "Medicamente"),
                   replicate(n = 61, "No medicamento"))



ggplot()  + 
  geom_line(data = confianza[1:61,],
            mapping = aes(x = EDAD, y = lwr),
            col = "#E4007C",
            linetype = 5) +
  geom_ribbon(data = confianza[1:61,],
            mapping = aes(x = EDAD, 
                          ymin = lwr, 
                          ymax = upr,
                          fill = Tipo),
            alpha = 0.2)+
  geom_line(data = confianza[1:61,],
            mapping = aes(x = EDAD, y = upr),
            col = "#E4007C",
             linetype = 5) +
  geom_line(data = confianza[62:122,],
            mapping = aes(x = EDAD, y = lwr),
            col ="#090088",
             linetype = 5)+
  geom_line(data = confianza[62:122,],
            mapping = aes(x = EDAD, y = upr),
            col = "#090088",
             linetype = 5) + 
  geom_ribbon(data = confianza[62:122,],
              mapping = aes(x = EDAD, 
                          ymin = lwr, 
                          ymax = upr,
                          fill = Tipo),
              alpha = 0.2)+
  theme_bw() + 
  labs(title = "Intervalos de confianza", 
       x = "Edad",
       subtitle = "Anticuerpos según edad y su tratamiento",
       y = "Nivel de anticuerpos") + 
  geom_function(fun = control, col ="#090088" ) +
geom_function(fun = med, col = "#E4007C")

Podemos afirmar con una confianza del 90% que para las personas entre las edades de 30 y 60 años que tomaron el medicamento es mayor el nivel promedio de anticuerpos que para las personas que no se medicaron.

Ejemplo de preguntas de regresión lineal multiple

Fuentes Tapia Ivan Baruch

Aclaración

Ejemplo 1

I.

II.

III.

IV.

V.

VI.

Ejemplo 2

I.

II.

III.

IV.

Reducción del modelo.

V

VI

VII.

Ejemplo 3

I.

II.

III.

IV

V.

VI.

Pruebas de hipótesis.

Intervalos de confianza.