---
title: "Práctica 5"
output: html_document
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Datos binomiales
## Ejercicio 1
En el archivo chdage.dat aparecen los datos de edad (age) y de presencia o ausencia de una enfermedad coronaria (chd) de 100 sujetos seleccionados para participar en un estudio. Nuestro interés se centra en estudiar la relación entre la edad y la presencia de la enfermedad coronaria. 
 

Importamos los datos:

```{r chdage, message=FALSE}
chdage<-read.table("./datos/chdage.dat",head=T)
```


Ajustamos a estos datos un modelo de regresión logística, 
$$
logit(p_i)=\log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_1\, age_i
$$

y damos las estimaciones de los parámetros, así como intervalos de confianza al 95% para los mismos.

```{r chdage ajuste logit, message=FALSE}
mlg1<-glm(chd~age, family=binomial(link="logit"),data=chdage)
summary(mlg1)
confint(mlg1)
```

Contrastamos la bondad de ajuste del modelo:
$$
\begin{array}{l}
H_0: \mbox{el modelo ajusta bien a los datos}\\
H_1: \mbox{el modelo no ajusta bien a los datos}
\end{array}
$$


```{r chdage bondad ajuste logit, message=FALSE}
1-pchisq(deviance(mlg1),df.residual(mlg1))
```


Contrastemos la edad como factor de riesgo para la presencia de la enfermedad coronaria mediante una comparación de devianzas:

```{r chdage significación age logit, message=FALSE}
mlg0<-glm(chd~1, family=binomial(link="logit"),data=chdage)
anova(mlg0, mlg1, test="Chisq")
```

¿Qué probabilidad predice el modelo de presentar la enfermedad para una
persona de 28 años? ¿Y para uno de 76?

```{r chdage prdicciones logit, message=FALSE}
predict(mlg1,data.frame(age=28))
predict(mlg1,data.frame(age=28),type="response")
predict(mlg1,data.frame(age=76),type="response")
```

Repitamos el análisis con el modelo probit
$$
F^{-1}(p_i)=\beta_0+\beta_1\, age
$$
con $F$ la función de distribución de una $N(0,1)$.

```{r chdage probit, message=FALSE}
#ajuste del modelo
mlg2<-glm(chd~age, family=binomial(link="probit"),data=chdage)
summary(mlg2)
#intervalos de confianza
confint(mlg2)
#contraste de bondad de ajuste
1-pchisq(deviance(mlg2),98)
#significación de la edad
mlg0<-glm(chd~1, family=binomial(link="probit"),data=chdage)
anova(mlg0, mlg2,test="Chisq")
#predicciones
predict(mlg2,data.frame(age=28),type="response")
predict(mlg2,data.frame(age=76),type="response")
```


## Ejercicio 2
En enero de 1986, el transbordador espacial Challenger explotó poco después del lanzamiento. Se inició una investigación sobre la causa del accidente y la atención se centró en las juntas de goma (orings)  de los cohetes aceleradores. A temperaturas más bajas, la goma se vuelve más frágil y es un sellador menos eficaz. En el momento del lanzamiento, la temperatura era de 31°F. ¿Podría haberse predicho el fallo de las juntas? En las 23 misiones anteriores del transbordador de las que existen datos, se registraron algunas evidencias de daños por soplado y erosión en algunas juntas. Cada transbordador tenía dos aceleradores, cada uno con tres juntas. Para cada misión, conocemos el número de juntas (de un total de seis) que mostraban algún daño y la temperatura de lanzamiento. Los datos se encuentran en el archivo orings.dat

Vamos a ajustar modelos de regresión logística y regresión probit a estos datos, obtener las ecuaciones que nos dan estos modelos, comprobar la bondad del ajuste y la significación de la temperatura a la hora de predecir la probabilidad de rotura de las juntas de goma. Por último los utilizaremos para hacer predicciones.



Importamos los datos:

```{r orings, message=FALSE}
orings<-read.table("./datos/orings.dat",head=T, sep="\t")
```


En este ejercicio, la variable respuesta es el número de juntas dañadas (damage) y la variable predictora la temperatura. Se podría pensar que si estimamos la probabilidad de fallo en cada misión, se podría ajustar un modelo lineal normal a los datos:

```{r orings lineal, message=FALSE}
plot (damage/6 ~ temp, data= orings, xlim=c(25,85), ylim = c(0,1), xlab="Temperature", ylab="Prob of damage")
lmod <- lm(damage/6 ~ temp, orings)
abline(lmod)
```

Es obvio que el modelo lineal normal no es adecuado. De hecho, la variable respuesta es discreta y no parece apropiado ajustar una distribución normal. Puesto que cada cohete tiene seis juntas de goma, y la variable respuesta es el número de juntas dañadas, la distribución más adecuada es la binomial con $n=6$.

Como para la distribución binomial tenemos varias funciones de enlace, probemos primero con la función logit (***regresión logística***)
$$
logit(p_i)=\log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_1 temp_i
$$
```{r orings logística, message=FALSE}
logitmod <- glm(cbind(damage,6-damage) ~ temp, family=binomial(link="logit"),data=orings)
summary(logitmod)
```
También podemos ajustar el modelo tomando como función de enlace la inversa de la función de distribución $N(0,1)$ (***regresión probit***)
$$
F^{-1}(p_i)=\beta_0+\beta_1\, temp_i
$$
con $F$ la función de distribución de una $N(0,1)$. Ajustemos este modelo:
```{r orings probit, message=FALSE}
probitmod <- glm(cbind(damage,6-damage) ~ temp, family=binomial(link="probit"),data=orings)
summary(probitmod)
```
En este caso en particular, es posible "dibujar" los ajustes de ambos modelos:


```{r oring logit vs probit, message=FALSE}
plot (damage/6 ~ temp, orings, xlim=c(25,85), ylim=c(0,1), xlab="Temperature", ylab="Prob of damage")
x <- seq(25,85,1)
ilogit<-function(x){exp(x)/(1+exp(x))}
lines(x,ilogit(11.6630 - 0.2162*x))
lines(x, pnorm(5.5915 - 0.1058*x), lty=2)
```

Calculemos los intervalos de confianza:

```{r intervalos oring, message=FALSE}
confint(logitmod)
confint(probitmod)
```

Utilizaremos la devianza residual para contrastar la bondad de ajuste de los modelos:
$$
\begin{array}{l}
H_0: \mbox{el modelo ajusta bien a los datos}\\
H_1: \mbox{el modelo no ajusta bien a los datos}
\end{array}
$$
Los p-valores para este contraste en ambos modelos son
```{r bondad de ajuste oring, message=FALSE}
#regresión logística
deviance(logitmod)
1-pchisq(deviance(logitmod),21)
#regresión probit
deviance(probitmod)
1-pchisq(deviance(probitmod),21)
```

Utilicemos ambos modelos para predecir la probabilidad de rotura de una junta a una temperatura de 31 ºF:

```{r predicciones oring, message=FALSE}
predict(logitmod,data.frame(temp=31),type="response")
predict(probitmod,data.frame(temp=31),type="response")
```


## Ejercicio 3
El archivo icu.dat contiene los datos de 200 individuos que formaron parte deun estudio mucho mayor sobre la supervivencia de los pacientes en una Unidad de Cuidados Intensivos (UCI). La descripción de cada una de las variables aparece en la hoja anexa.
El objetivo del estudio era la predicción de la probabilidad de supervivencia de un paciente tras ser ingresado en la UCI.

Importamos los datos:

```{r icu , message=FALSE}
icu<-read.table("./datos/icu.dat",head=T)
colnames(icu)
```


Vamos a ajustar un modelo de regresión logística de la variable sta sobre las variables age, can, cpr(Cardiopulmonary Resuscitation), inf, y race. 
$$
logit(p_i)=\log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_1\, age+\beta_2\, can +\beta_3\, cpr+ \beta_4\, inf +\beta_5\,race2+\beta_6\,race3
$$

Obtenemos también los intervalos de confianza al 95% para dichos parámetros.

```{r icu logit, message=FALSE}
mlg3<-glm(sta~age+can+cpr+inf+factor(race), family=binomial(link="logit"),data=icu)
summary(mlg3)
confint(mlg3)
```

Contrastamos la bondad de ajuste del modelo 
$$
\begin{array}{l}
H_0: \mbox{el modelo ajusta bien a los datos}\\
H_1: \mbox{el modelo no ajusta bien a los datos}
\end{array}
$$


```{r bondad ajuste icu logit, message=FALSE}
1-pchisq(deviance(mlg3),df.residual(mlg3))
```

```{r}
mlg3b<-glm(sta~age+can+cpr+inf, family=binomial(link="logit"),data=icu)
anova(mlg3b,mlg3, test="Chisq")
```

Aplicamos un método backward de selección de variables:

```{r backward icu logit, message=FALSE}
step(mlg3, direction='backward', criterion='AIC')
```

Aplicamos un método forward de selección de variables:

```{r forward icu logit, message=FALSE}
mlg0<-glm(sta~1, family=binomial(link="logit"),data=icu)
step(mlg0,scope=list(upper=mlg3,lower=mlg0), direction='forward', criterion='AIC')
```

Ajustemos el modelo con el que nos quedamos tras el procedimiento de selección:

```{r icu logit 2,message=FALSE}
mlg4<-glm(sta~age+cpr+inf, family=binomial(link="logit"),data=icu)
summary(mlg4)
```

Utilicemos este modelo para realizar predicciones:

```{r icu logit 2 predicciones,message=FALSE}
predict(mlg3,data.frame(age=47,race=1,can=0,cpr=1,inf=1),type="response")
predict(mlg4,data.frame(age=47,race=1,can=0,cpr=1,inf=1),type="response")
```