---
title: "Práctica 5"
output: html_document
date: "`r Sys.Date()`"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# Datos de Poisson
## Ejercicio 1
El número de muertes por Sida en Australia en períodos de tres meses entre 1983 y 1986 aparecen en el archivo sida.dat.
Se pretende ajustar un modelo de regresión de Poisson a dichos datos, determinar los coeficientes e intervalos de confianza al 95% para los mismos así como la bondad del ajuste. Como conclusión tenemos que ver si la variable periodo tiene relación significativa con el número de muertes.
Importamos los datos:
```{r sida, message=FALSE}
sida<-read.table("./datos/sida.dat",head=T)
colnames(sida)
plot(sida$periodo,sida$muertos)
abline(lm(muertos~periodo,data=sida))
```
Ajustamos a estos datos modelo de Poisson, con función de enlace el logaritmo:
$$
\log(\lambda_i)=\beta_0+\beta_1\, i
$$
y damos las estimaciones de los parámetros, así como intervalos de confianza al 95% para los mismos.
```{r sida ajuste log, message=FALSE}
mlg1<-glm(muertos~periodo,family=poisson(link="log"),data=sida)
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 sida bondad ajuste log, message=FALSE}
1-pchisq(deviance(mlg1),df.residual(mlg1))
```
El ajuste no parece muy bueno, cambiamos la función de enlace por la raíz cuadrada:
$$
\sqrt{\lambda_i}=\beta_0+\beta_1\, periodo_i
$$
y damos las estimaciones de los parámetros, así como intervalos de confianza al 95% para los mismos.
```{r sida ajuste sqrt, message=FALSE}
mlg2<-glm(muertos~periodo,family=poisson(link="sqrt"),data=sida)
summary(mlg2)
```
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 sida bondad ajuste sqrt, message=FALSE}
1-pchisq(deviance(mlg2),df.residual(mlg2))
```
¿Qué probabilidad predice el modelo de presentar la enfermedad para el periodo 15?
```{r sida prdicciones sqrt, message=FALSE}
predict(mlg2,data.frame(periodo=16),type="response")
```
## Ejercicio 2
Para 30 islas de las Galápagos, disponemos de un recuento del número de especies de plantas encontradas en cada isla y del número de especies que son endémicas de dicha isla. También contamos con cinco variables geográficas para cada isla y dos climáticas. Se trata de determinar qué variables geográficas influyen en el número de tortugas.
Vamos a ajustar un modelo de regresión de Poisson, comprobar la bondad del ajuste y la significación de las distintas variables geográficas y climáticas. Por último lo utilizaremos para hacer predicciones.
Importamos los datos:
```{r gala, message=FALSE}
gala<-read.table("./datos/gala.dat",head=T, sep="\t")
colnames(gala)
```
En este ejercicio, la variable respuesta es el número de especies de plantas. Ajustamos un modelo de regresión de Poisson de esta variable frente las variables geográficas y climáticas.
$$
logit(\lambda_i)=\beta_0+\beta_1 Area +\beta_2 Elevation + \beta_3 Nearest + \beta_4 Scruz + \beta_5 Adjacent+\beta_6 Temperatura+\beta_7Humidity
$$
```{r gala log, message=FALSE}
mlg3 <- glm(Species ~ Area+ Elevation+ Nearest + Scruz+ Adjacent+Temperature+Humidity, family=poisson(link="log"),data=gala)
summary(mlg3)
```
Calculemos los intervalos de confianza:
```{r intervalos gala, message=FALSE}
confint(mlg3)
```
Utilizaremos la devianza residual para contrastar la bondad de ajuste:
$$
\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 de ajuste gala, message=FALSE}
1-pchisq(deviance(mlg3),df.residual(mlg3))
```
Aplicamos un método backward de selección de variables:
```{r backward gala, message=FALSE}
step(mlg3, direction='backward', criterion='AIC')
```
Aplicamos un método forward de selección de variables:
```{r forward gala, message=FALSE}
mlg0<-glm(Species~1, family=poisson(link="log"),data=gala)
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 gala 2,message=FALSE}
mlg4<-glm(formula = Species ~ Elevation + Adjacent + Area + Scruz + Nearest + Temperature, family = poisson(link = "log"), data = gala)
summary(mlg4)
```
Se trata no obstante, de un modelo poco conveniente para realizar predicciones.
## Ejercicio 3
Los datos del archivo cancer.dat muestran el número de supervivientes en un estudio sobre cáncer de pulmón. En los datos aparecen los intervalos en que se hizo el seguimiento, el tipo de histología, el estadío de la enfermedad y el período de seguimiento (en meses). La variable de interés el el recuento del número de muertes. Estos recuentos son tratados como variables de Poisson.
Vamos a ajustar un modelo log lineal a estos datos, comprobar la bondad del ajuste y la significación de los distintos factores.
Importamos los datos:
```{r cancer, message=FALSE}
cancer<-read.table("./datos/cancer.dat",head=T)
colnames(cancer)
```
Ajustamos el modelo log lineal siguiente
$$
\log(\mu_{ijk})=\beta_0+ \beta_1 risktime_{ijk} + \tau_i +\delta_j+\gamma_k+ (\tau\delta)_{ij}+(\tau\gamma)_{ik}+(\delta\gamma)_{jk}
$$
```{r cancer ajuste , message=FALSE}
mlg5 <- glm(count ~ risktime+ factor(histology) + factor(stage) + factor(time) + factor(histology):factor(stage) + factor(histology):factor(time) + factor(stage):factor(time), family = poisson(link = log), data=cancer)
summary(mlg5)
anova(mlg5,test="Chisq")
```
Utilizaremos la devianza residual para contrastar la bondad de ajuste:
$$
\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 de ajuste cancer, message=FALSE}
1-pchisq(deviance(mlg5),df.residual(mlg5))
```
Aplicamos un método backward de selección de variables:
```{r backward cancer, message=FALSE}
step(mlg5, direction='backward', criterion='AIC')
```
## Ejercicio 4
En el archivo clot.dat aparecen los datos del tiempo de coagulación (time, en segundos), para plasma normal diluido con diferentes concentraciones de plasma sin protombina (conc). La coagulación fue inducida por dos lotes de tromboplastina (lot).
Se pretende ajustar un modelo de regresión de Gamma, usando la inversa como función de enlace y el logaritmo de la concentración como variable predictora. Una vez ajustado, hay que determinar las estimaciones para los coeficientes del modelo, especificando si hay diferencia significativa para ambos lotes.
```{r clot,message=FALSE}
clot<-read.table("./datos/clot.dat",head=T)
colnames(clot)
```
La variable respuesta es el tiempo de coagulación y se asume con distribución Gamma: $Y_i\sim \Gamma(\phi, \theta_i)$
$$
\mu_i=E[Y_i]=\phi\theta_i\ ,\quad Var[Y_i]=\phi\theta_i^2=\frac{\mu_i^2}{\phi}
$$
Como podemos ver, el parámetro $\phi$, sin ser el principal juega un papel importante en el modelo y habrá que tenerlo en cuenta.
Ajustamos un modelo de regresión Gamma con función de enlace inversa:
$$
\frac{1}{\mu_i}=\beta_0+\beta_1x_i+\tau lot.two_i + \delta (x_i*lot.two_i)\quad i=1,\dots,n
$$
```{r ajuste clot, message=FALSE}
mlg6<-glm(time~log(conc)+lot+log(conc):lot,family=Gamma(link="inverse"),data=clot)
summary(mlg6)
confint(mlg6)
```
Contrastamos la significación del lote, tanto en su efecto principal como en su interacción con el logaritmo de la concentración:
```{r significación lot clot, message=FALSE}
mlg7<-glm(time~log(conc),family=Gamma,data=clot)
summary(mlg7)
anova(mlg7,mlg6,test="Chisq")
```
Contrastemos la bondad de ajuste del modelo. En este caso, la devianza hay que dividirla por la estimación del parámetro de dispersión para que se pueda comparar con la chi cuadrado teórica.
```{r bondad de ajuste clot, message=FALSE}
1-pchisq(0.029401/0.002129707,14)
```
Por último, vamos a ver cuánto tardaría en coagularse un plasma del lote one con concentración 27.
```{r sida predicciones sqrt, message=FALSE}
predict(mlg6,data.frame(conc=27, lot="one"),type="response")
```