# Modelos Lineales Generalizados.

# Se pretende modelizar, para enfermos de Leucemia,
# la relación existente entre el número de semanas transcurridas
# hasta la muerte de un paciente desde su diagnóstico y el
# log_10 del número inicial de leucocitos que presentaba. Para
# ello se han tomado 17 pacientes con Leucemia, cuyos datos están en
# el archivo leucemia.dat

# Modelo:
# 1) La distribución exponencial es utilizada habitualmente para
# describir tiempos de supervivencia:  Y_i ~ Exp(theta_i),
# f_{theta_i}(y_i)=theta_i*exp(-y_i*theta_i), i=1,...,n.
#
# 2) Una posible especificación para E[Y] es 
# E[Y_i]=exp(beta_0 + \beta_1*x_i), i=1,...,n, que asegura que E[Y] 
# es no negativa para todos los valores de los parámetros y todos los valores 
# de x. En este caso la función de enlace es g(x)=log(x).

datos<-read.table("leucemia.dat",head=T)
attach(datos)

plot(y~x)

res<-glm(y~x,family=Gamma(link="log"))
summary(res)

# Adecuación del modelo
1-pchisq(res$deviance,res$df.residual)

# Contraste H_0: beta1=0

res1<-glm(y~1,family=Gamma(link="log"))
anova(res1,res,test="Chisq")


xx<-seq(min(x),max(x),len=1000)

lines(xx,exp(coefficients(res)[1]+coefficients(res)[2]*xx),col=2)



# Regresión de Poisson

# 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. 
# Ajustar un modelo de regresión de Poisson a dichos datos.

datos<-read.table("sida.dat",head=T)
attach(datos)

plot(periodo,muertos)

res<-glm(muertos~log(periodo),family=poisson(link="log"))
summary(res)

1-pchisq(res$deviance,res$df.residual)

pp<-seq(min(periodo),max(periodo),len=1000)
lines(pp,exp(res$coefficients[1]+res$coefficients[2]*log(pp)),col=2)




# Modelos Dosis-Respuesta
#
# Los vectores dosis, tumor y total contienen los datos de un experimento efectuado
# con una serie de grupos de ratas que fueron expuestas a diferentes niveles de 
# concentración de ETU (Ethylene thiourea) observándose el número de ratas que 
# desarrollaron tumores tiroideos (Graham et al. (1975)).


dosis<-c(0,5,25,125,250,500)
tumor<-c(2,2,1,2,16,62)
total<-c(72,75,73,73,69,70)

plot(dosis,tumor/total)

# Modelo logístico
res<-glm(cbind(tumor,total-tumor)~dosis, family=binomial(link="logit"))
summary(res)

1-pchisq(res$deviance,res$df.residual)

x<-seq(0,500,len=1000)
lines(x,1/(1+exp(-(res$coefficients[1]+res$coefficients[2]*x))),col=2)


# Modelo probit
res1<-glm(cbind(tumor,total-tumor)~dosis, family=binomial(link="probit"))
summary(res1)

1-pchisq(res1$deviance,res1$df.residual)


lines(x,pnorm(res1$coefficients[1]+res1$coefficients[2]*x),col=3)

# Modelo log-log complementario
res2<-glm(cbind(tumor,total-tumor)~dosis, family=binomial(link="cloglog"))
summary(res2)

1-pchisq(res2$deviance,res2$df.residual)


lines(x,1-exp(-exp(res2$coefficients[1]+res2$coefficients[2]*x)),col=4)