Taller 4 Modelos

---
title: "Taller 4"
author: "Jorge Moreno, Aarón Calderón"
output:
  pdf_document: default
  html_notebook: default
---

# Ejercicio 1: Dataset: cigarrillos - Suma de cuadrados extra y multicolinealidad

La dataset “cigarette” contiene información de 25 cigarrillos:

**Name:** Marca del cigarillo

**Tar:** Contenido de Alquitrán en miligramos

**Nicotine:** Contenido de Nicotina en miligramos

**Weight:** Peso del cigarrillo en gramos

**CO:** Contenido de Monóxido de Carbono en miligramos

Queremos modelar el contenido de Monóxido de Carbono como una función de las variables explicativas.

## Pasos introductorios:

Importamos el set de datos:

```{r}
#cigarette <- read.csv("C:/Users/USUARIO/Downloads/cigarette.txt", sep="")
cigarette <- read.csv("C:/Users/aleja/Desktop/cigarette.txt", sep="")
attach(cigarette)

```

## Exploración gráfica y estadística descriptiva:

La exploración gráfica y la estadística descriptiva del set de datos es el primer paso de la modelación. Con
ello podemos conocer las variables y hacernos una idea de los tipos de relaciones que existen.

El histograma de la variable de respuesta Y nos ayuda a saber si el supuesto de normalidad es adecuado o si
deberíamos considerar alguna transformación de Y.

```{r}
hist(CO)
```

En este caso el histograma tiene una forma de campana por lo que no se sospecharía de violación grave del
supuesto de normalidad. Recuerde sin embargo que el supuesto de la normalidad de Y es condicional en los
valores de X, por lo que el histograma de los valores observados de y no es evidencia de cumplimiento del
supuesto. Se deben siempre analizar los residuos.

Para graficar la relación de todas las variables de la dataset se puede usar

```{r}
pairs(cigarette[, -1], panel = function(x, y) {
  points(x, y)
  lines(lowess(x, y), col = "red")
})
```

### Pregunta 1. Estime el siguiente modelo y obtenga el resumen del modelo y la tabla ANOVA con la suma de cuadrados extra. ¿Observa algún resultado extraño en la estimación del modelo? (No centre las variables)

$$CO=\beta_0+\beta_1Tar+\beta_2Nicotine+\beta_3Weight+\beta_4Weight^2+\epsilon$$

```{r}
modelo1=lm(CO~Tar+Nicotine+Weight+I((Weight)^2),data=cigarette)

summary(modelo1)
```
En primera instancia, podemos apreciar que el $\beta_2$ corrrespondiente a Nicotina es un valor negativo, o en otras palabras con pendiente negativa, sin embargo, si apreciamos las graficas correspondientes a esta, vemos que contradice dicho resultado.

Adicionalmente, la prueba T da un valor p relativamente alto, lo que no nos permitiría rechazar la hipótesis de $\beta_{2}=0$, pero en el gráfico se ve que sí existe una correlación entre Nicontina y Alquitrán.


```{r}
anova(modelo1)
```
Notamos que la suma cuadrática de regresión de Weight da 0, esto ocurre si para cada valor $\hat y_i=\bar y$ lo cual conllevaría que el modelo lineal de Weight con la variable de respuesta es la recta constante $\bar y$, por los gráficos notamos que no es así.

Se repiten los resultados con la prueba F para Nicotine, indicando que no existe evidencia estadística para rechazar $\beta_{2}$.

Encontramos que con el summary, según la prueba T, para $\beta_{3}$, correspondiente al peso, existe evidencia para rechazar la hipótesis que dice que dice que $\beta_{3}=0$, sin embargo, la prueba F contradice estos resultados.

### Pregunta 2. Estime el siguiente modelo sin Tar y obtenga el resumen del modelo y la tabla ANOVA con la suma de cuadrados extra. ¿Que observa comparando este modelo al modelo de la pregunta 1? ¿Cómo explicaría el cambio? (No centre las variables)

$$CO=\beta_0+\beta_2Nicotine+\beta_3Weight+\beta_4Weight^2+\epsilon$$

```{r}
modelo2=lm(CO~Nicotine+Weight+I((Weight)^2),data=cigarette)

summary(modelo2)
anova(modelo2)
```

Del nuevo modelo, ahora se obtuvo el $\beta_2$ positivo correspondiente a Nicotina, además, en este modelo solo para Nicotina tenemos evidencia de rechazar la hipotesis nula, la cual indica que el $\beta=0$. Se obtuvo un menor coeficiente de determinación, por lo que explica un poco menos la variabilidad de la variable CO. Esto se debe a que al quitar la variable Tar, la mayor parte de la variabilidad fue absorbida por Nicotina.



### Pregunta 3. Obtenga los factores de inflación de la varianza para detectar problemas de multicolinealidad en la matriz X del modelo 1. (Puede usar la función vif de la librería car

```{r}
library(car)
vif(modelo1)
```

Por las reglas de thumb, tenemos pruebas evidentes que existe multicolinealidad, ya que al menos uno de los VIF es mayor a 10 y el promedio de los VIF es considerablemente mayor a 1. Obviamente, por los VIF similares, concluimos que las variables Weight y Weight^2 están correlacionadas.


### Pregunta 4. Corrija el problema de multicolinealidad, removiendo Nicotine y centrando weight.

$$CO=\beta_0+\beta_1Tar+\beta_3Weight+\beta_4Weight^2+\epsilon$$

```{r}
modelo4=lm(CO~Tar+Weight+I((Weight-mean(Weight))^2),data=cigarette)

summary(modelo4)
vif(modelo4)
```

Notamos como reducieron considerablemente los VIF, ya no superan a 10 y su promedio ya no es considerablemente mayor a 1.

# Ejercicio 2: FEV: Interacciones

```{r}
#fev <- read.csv("C:/Users/USUARIO/Downloads/fev.txt", sep="")
fev<- read.csv("C:/Users/aleja/Desktop/fev.txt", sep="")

fev_subset=fev[,1:3]
attach(fev_subset)

pairs(fev_subset, panel = function(x, y) {
  points(x, y)
  lines(lowess(x, y), col = "red")
})

```

En base a talleres anteriores sabemos que un buen modelo para FEV en base a la edad y a la altura es un
modelo aditivo con efecto cuadrático de altura y estimado por mínimos cuadrados ponderados.

```{r}
modelo0 = lm(FEV ~ age + height + I((height - mean(height))^2), data = fev_subset)
modelo0_ws = lm(abs(residuals(modelo0)) ~ fitted.values(modelo0))
ws = 1/fitted.values(modelo0_ws)^2

modelo_1 = lm(FEV ~ age + height + I((height - mean(height))^2), data = fev_subset,
weights = ws)

```

### Pregunta 5. Grafice los residuos ponderados del modelo aditivo anterior vs el efecto interactivo de age y height para detectar si es necesario incluir un efecto de interacción.

```{r}
res_pond=sqrt(ws)*modelo_1$residuals

plot((age*height),res_pond)
```

Podemos ver que están distribuidos alrededor del 0, sin patrón aparente. Adicionalmente, los supuestos tienen una gran variabilidad, por lo que es una buena idea incluir el efecto de interacción.

### Pregunta 6. Estime un modelo incluyendo el efecto de interacción entre age y height.

$$FEV=\beta_0+\beta_1age+\beta_2height+\beta_3age\hspace{0.05cm}height+\epsilon$$

El modelo puede estimarse por:

```{r}
modelo_2=lm(FEV~age+height+age:height,data=fev_subset)
```


### Pregunta 7. Analice los residuos del modelo anterior.


```{r}
plot(modelo_2$residuals,ylim=c(-3,3))
abline(2.5,0)
abline(-2.5,0)
```

Como puede verse, los residuos están dsitribuidos de forma aleatoria alrededor del cero. Adicionalmente, no parece haber valores aberrantes.

### Pregunta 8. ¿Es el término de interacción importante para explicar la respuesta?. Responda a esta pregunta utilizando una prueba F parcial.

Lo que queremos probar es que:

$$H_0: \beta_3=0 \hspace{0.5cm} H_1:\beta_3\neq0$$

```{r}
modelo_red=lm(FEV~age+height,data=fev_subset)
modelo_com=modelo_2

anova(modelo_red,modelo_com)
```

Por el valor p de la prueba F parcial podemos rechazar la hipótesis nula, por lo que el término de interacción si es importante para explicar la respuesta.


### Pregunta 9. En base al último modelo, ¿A cuánto se estima que cambie el promedio de FEV por un aumento de una pulgada en la altura de los niños?

$$FEV=\beta_0+\beta_1age+\beta_2height+\beta_3age\hspace{0.05cm}height+\epsilon$$

$$
\begin{aligned}
E\left[FEV|age=age,height=1+height\right]&-E\left[FEV|age=age,height=height\right]\\
\beta_0+\beta_1age+\beta_2(1+height)+\beta_3age\hspace{0.05cm}(1+height)&-\beta_0-\beta_1age-\beta_2height-\beta_3age\hspace{0.05cm}height\\
&\beta_2+\beta_3age
\end{aligned}
$$

Por ende, el cambio estimado en promedio de FEV por aumento de una pulgada de altura en los niños es: $\beta_2+\beta_3age$