Modelos taller 3 completo

---
title: "Taller 3_Modelos I"
author: "Jorge Moreno"
output: pdf_document
---


# Ejercicio 2: Dataset: FEV

La dataset FEV incluye observaciones de 654 niños a los que se les mide el Volumen Espiratorio Forzado (FEV: forced expiratory volume). En este ejercicio estamos interesados en modelar la relación de FEV y las siguientes variables.

**age**: la edad del niño en años

**height**: la altura del niño en pulgadas

Importamos el set de datos:

```{r}
fev <- read.csv("C:/Users/aleja/Desktop/fev.txt", sep="")
fev_subset=fev[,1:3]
```

### Pregunta 1. Realice una exploración gráfica de la relación entre FEV y la edad y la altura del niño. ¿Se ajustaría un modelo aditivo de primer orden o algún orden superior?

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

Puede verse que existe una curvatura en la relación entre FEV y edad, y entre FEV y altura, por lo que, en este caso, se ajustaría mejor un modelo de orden superior.

### Pregunta 2. Ajuste el siguiente modelo aditivo de primer orden, analice los residuos del modelo y detecte la curvatura y heterocedasticidad en los residuos:

$$
FEV_i = \beta_0 + \beta_1 age_i + \beta_2 height_i + \epsilon_i
$$

```{r}
modelo1=lm(FEV~age+height,data=fev_subset)
par(mfrow=c(2,2))
plot(modelo1)


par(mfrow=c(1,1))
plot(modelo1$residuals,ylim=c(-3,3))+abline(2.5,0)+abline(-2.5,0)
```
En el gráfico de los valores ajustados contra los residuos, podemos ver que la varianza incremente conforme los valores ajustados aumentan, por lo que existen problemas de heterocedasticidad. También podemos ver una ligera curvatura.

Adicionalmente, en el gráfico Q-Q, podemos ver que tenemos una situación de simétrica con colas pesadas, por lo que no hay normalidad.

No hay valores aberrantes.

### Pregunta 3. Ajuste el siguiente modelo aditivo de segundo orden, analice los residuos del modelo:

$$
 FEV_i = \beta_0 + \beta_1 age_i + \beta_2 height_i + \beta_3 age^2_i + \beta_4 height^2_i + \epsilon_i
$$

```{r}
modelo2=lm(FEV~age+height+I((age-mean(age))^2)+I((height-mean(height))^2),data=fev_subset)
par(mfrow=c(2,2))
plot(modelo2)

par(mfrow=c(1,1))
plot(modelo2$residuals,ylim=c(-3,3))+abline(2.5,0)+abline(-2.5,0)
```

El problema de simetría con colas pesadas persiste, por lo que no se ha solucionado la no normalidad.

Sin embargo, ya no está presente la curvatura en el gráfico de valores ajustados contra residuos. Lamentablemente, todavía es posible observar heterocedasticidad.

### Pregunta 4. Realice una prueba de hipótesis para los efectos cuadráticos de age y height mejoran el ajuste del modelo. Escriba el contraste de hipótesis, el estadístico de prueba, el valor p de la prueba y su conclusión.

$$
 H_0: \beta_3=\beta_4=0
\quad  \quad  \quad \quad  \quad  \quad \quad  \quad  \quad  H_1:
 \beta_3 \neq 0  \vee \beta_4 \neq 0
$$


```{r}
SCEpq=sum(modelo1$residuals^2)
SCEp=sum(modelo2$residuals^2)
q=2
n=length(fev_subset$FEV)
p=5
(F0=((SCEpq-SCEp)/(q))/(SCEp/(n-p)))
(pf(F0,q,n-p,lower.tail = F))
```
Una forma alternativa puede ser:

```{r}
anova(modelo1,modelo2)
```

Si suponemos una significancia del 5%, existe evidencia estadística para rechazar la hipótesis nula, esto es, que $\beta_{3}=\beta_{4}=0$.

### Pregunta 5. Compare el modelo 1 y 2 en base al coeficiente de determinación ajustado. ¿Qué modelo prefiere?

```{r}
sum1=summary(modelo1)
sum2=summary(modelo2)

R2=cbind(sum1$r.squared, sum2$r.squared)
R2_adj=cbind(sum1$adj.r.squared, sum2$adj.r.squared)

colnames(R2)=c("modelo1","modelo2")
colnames(R2_adj)=c("modelo1","modelo2")


R2
R2_adj
```
Como el modelo 2 tiene un mejor coeficiente de determinación ajustado, se prefiere este último, ya que explica mejor la variabilidad.

### Pregunta 6. Estimemos el modelo 2 usando los mínimos cuadrados ponderados. Analice los residuos $r^{(w)}$

```{r}
modelo3=lm(abs(modelo2$residuals)~modelo2$fitted.values)
w=1/(modelo3$fitted.values^2)
modelo4=lm(FEV~age+height+I((age-mean(age))^2)+I((height-mean(height))^2),data=fev_subset,weights=w)
par(mfrow=c(2,2))
plot(modelo4)
par(mfrow=c(1,1))


plot(modelo4$residuals,ylim=c(-3,3))+abline(2.5,0)+abline(-2.5,0)
```


Vemos que no hay ningún patrón presente entre los residuos ponderados y los valores ajustados. Los supuestos de heterocedasticidad y de normalidad no se cumplen todavía, pues la varianza aumenta cuando crecen los valores ajustados, y todavía hay colas pesadas.

### Pregunta 7. Utilice las funciones plot3d() y persp3d() de la librería "rgl" para graficar la superficie de respuesta del último modelo.

```{r}
#install.packages("rgl")
library(rgl)
summary(fev_subset$age)
x=c(3:19)
summary(fev_subset$height)
y=c(46:74)
xy=data.frame(expand.grid(x,y))
colnames(xy)=c("age","height")
pxy=data.frame(predict(modelo2,xy))
wxy=1/(predict(modelo3,pxy)^2)
z=predict(modelo4,xy,weights = w_xy)

plot3d(fev_subset$ag,fev_subset$height,fev_subset$FEV,xlab="age",ylab="height",zlab="FEV")
persp3d(x,y,z,add=TRUE,color="blue")
```


### Pregunta 8. ¿Cuánto sería el FEV de un niño de 5 años y 53 pulgadas? Proporcione un intervalo de predicción e interprételo.

```{r}
predict(object=modelo2,newdata=data.frame(height=53,age=5),interval=c("prediction"),level=0.95,df=n-5)
```
El intervalo nos dice que, con un 95% de predicción, si se tiene un niño de 5 años y de 53 pulgadas, podemos predecir que el FEV del niño está entre 0.54 y 2.12.