#############################################################################################
### Curso modelos lineales generales, generalizados y mixtos en ecologia IE-UNAM          ###
### Modelos lineales generalizados con distribucion binomial                              ###
### Roberto Munguia-Steyer                                    2011/11/24                  ###
#############################################################################################

R

library(MASS)
library(bbmle)
library(car) #util para estimar los logits empiricos
library(ISwR)


### 
## Distribucion binomial, num de exitos del total de intentos independientes

# a)Cual es la probabilidad de que una familia con 5 vastagos, 4 de ellos sean hijas? 

# 1) si la probabilidad de que un bebe nazca ninho o ninha es la mísma

pr0 <- dbinom(x=4, size=5, prob=0.5)
pr0

#equivalente a p(x) = n!/(X!(n-X)! p^x (1-p)^n-x

coefbin <- choose(5,4) # Coeficiente binomial, no importa el orden en que nacieron el hijo y las hijas
#con tal de que sean 4 hijas y un hijo

pr <- coefbin*(0.5^4)*((1-0.5)^1)
pr;pr0

# 2) La probabilidad de que que un bebe sea ninha en Mexico ha sido reportada de 0.487
#http://www.indexmundi.com/mexico/sex_ratio.html

psmex= 0.487
dbinom(x=4, size=5, prob=psmex)

# Valor esperado de la probabilidad, calculemos las probabilidades de tener de 0 a 5 hijas

phijas <- dbinom(x=0:5, size=5, prob=psmex)
phijas
names(phijas) <- 0:5
#Grafica
plot(0:5, phijas, type="h",lwd=2.5, xlab="Numero de hijas", col= "red", ylab= "Probabilidad")

# La suma de probabilidades de los posibles resultados exhaustivos y complementarios = 1

sum(phijas)

#Ahora simulemos que tenemos 2000 familias con cinco vastagos y queremos saber el valor esperado de hijas para cada familia
nfam =2000
hsim <- rbinom(nfam, size=5, prob=psmex)

#Frecuencia de hijas por familia con cinco descendientes
table(hsim)

#Frecuencia relativa de hijas por familia de n=5

frs <- table(hsim)/nfam
frs

#Valor esperado de numero de hijas por familia teorico E(y) = Np

ve <- sum(0:5* phijas)
ve

#Valor esperado de hijas por familia de cinco vastagos: datos simulados
mean(hsim)
var(hsim)

#Valor esperado: E(y) = Np, Varianza: var(y) = Np(1-p)

# Por lo general nosotros desconocemos los parametros de la distribucion probabilistica a partir de la cual se generaron nuestros datos
# y realizamos una inferencia estadistica para estimar estos parametros

#Numero de hijas en las 2000 familias
hsim 

#Los hijos son el complemento ya que todas las familias simuladas tuvieron 5 descendientes

msim <- 5 -hsim # numero de hijos varones por familia
unique(hsim+msim)

#Modelo lineal generalizado de distribucion binomial
m1 <- glm(cbind(hsim,msim)~1, binomial)
coef(m1)

# Probabilidad de tener una hija, estimada a partir de nuestra base de datos simulada
# plogis equivalente a utilizar una funcion llamada: logistica <-  function(x) exp(x)/(1+ exp(x))
# la funcion convierte los valores estimados del predictor lineal a valores de probabilidad

logistica <-  function(x) exp(x)/(1+ exp(x))

plogis(coef(m1))
logistica(coef(m1))

# Intervalos de confianza al 95%

# en logits: log(p)/1-p logit(y) ~ a
confint(m1) 

# en valores de probabilidad
logistica(confint(m1))

#Y si nuestro tamanho de muestra fuera diez veces menor

nfam2=200
hsim2 <- rbinom(nfam2, size=5, prob=psmex)

m2 <- glm(cbind(hsim2, 5-hsim2)~1,binomial)

logistica(coef(m1));logistica(coef(m2))

logistica(confint(m1));logistica(confint(m2))

# La incertidumbre asociada a la estimacion de los parametros decrece con el aumento
# de la informacion presente en el modelo (n mayor)

## Regresion logistica, base de datos besouro. Los datos describen el efecto de diferentes concentraciones de insectisida
## en la mortalidad de los escarabajos. Estimemos la probabilidad de la mortalidad dependiendo de la concentracion del insecticida.
## Funciones mle2 y glm

besor<-read.table("besouro.csv",sep=",", header=T) 
head(besor)
str(besor)

#graficando las proporciones de individuos muertos en escala logit, util para sugerir valores iniciales en la funcion mle2
bes.logit<-logit(besor$affected/besor$exposed)
plot(besor$conc,bes.logit)

##Ajuste de los modelos
#con glm
mrl1<-glm(cbind(affected,exposed-affected)~conc,data=besor,binomial(link="logit")) #funcion de ligacion logit es la default
summary(mrl1)

#con mle2, necesario asignar valores iniciales a los parameteros
mrl2<-mle2(affected~dbinom(logistica(a + b*conc), size = exposed), start=list(a =-3, b =0.5), data= besor)
summary(mrl2)

#Perfiles de verosimiliud e intervalos de confianza
perfil<-profile(mrl1)
confint(perfil)

#Graficando los datos (proporcion de escarabajos muertos en funcion de la conc. del insecticida) y la curva logística
par(las=1)
propor<-besor$affected/besor$exposed
plot(besor$conc,propor, cex=besor$exposed/15, ylab="Probabilidad de muerte", xlab="Conc. de insecticida" )
a.est<-coef(mrl1)[1]
b.est<-coef(mrl2)[2]
curve(logistica(a.est + b.est*x),add=T)

#sobredispersion, manera rapida y burda de detectarla Devianza res/gl residuales, correccion de los errores estandar
#los valores estimados de los parametros no cambian
summary(mrl1)

mrl3 <- update(mrl1,family= quasibinomial)

summary(mrl3) #Veamos la sobredispersion estimada y la variacion en los errores estandar
summary(mrl1) #Los errores estandar estimados sin sobredispersion eran mayores y por tanto mayor la prob. de cometer error del tipo 1

#Validacion del modelo
par(mfrow=c(2,2))
plot(mrl3)


# Ejercicio 11.1 
# Evalua el riesgo de obtener malaria en funcion de la edad y del nivel de anticuerpos en el cuerpo (en escala log). 
# La base de datos malaria se encuentra en el paquete (ISwR)
# 1. Efectua un glm con dist binomial y funcion liga logit, teniendo presentes como variables predictoras: edad, lanticuerpos y la interaccion de ambas.
# 2. Realiza seleccion de modelos para determinar cual es el modelo minimo adecuado, ¿que variable contiene?
# 3. Grafica los valores predichos empleando la funcion predict con el argumento type="response", por que empleamos este ultimo argumento.
# 4. Coteja que los valores predichos  del punto 3 coincidan con los valores predichos que obtenemos a partir de los coeficientes 
#    y la funcion que creamos llamada logistica que convierte los valores predichos de logits a valores estimados de probabilidad.
# 5. ¿Por que en este caso no evaluamos la sobredispersion?



data(malaria)
names(malaria) <- c("sujeto", "edad", "ac", "enfermedad")
head(malaria)

#Creando una la variable log(ac)
malaria$lac <- log(malaria$ac)

mbin <- glm (enfermedad ~ edad*lac, binomial, data=malaria)
summary(mbin)

mbin2 <- glm(enfermedad ~ lac, binomial, data=malaria)

#Graficando 

ie <- coef(mbin2)[1] 
pe <- coef(mbin2)[2] 

par(las=1)
plot(malaria$lac, malaria$enfermedad, ylab= "Probabilidad de enfermarse", xlab= "Conc. de anticuerpos (log)")
curve(logistica(ie + pe*x), add=T) #o la funcion plogis

#Funcion predict
par(las=1)
plot(malaria$lac, malaria$enfermedad, ylab= "Probabilidad de enfermarse", xlab= "Conc. de anticuerpos (log)")
range(malaria$lac)

#Generemos 100 valores de lac dentro del rango presente en nuestra base
seqx <- seq(min(malaria$lac), max(malaria$lac), length=100)
seqx

#Transformemos los valores de logits a valores de probabilidad con la funcion predict y el argumento type="response"
valpred <- predict(mbin2,list(lac= seqx), type="response")
valpred

#Grafica2
par(las=1)
plot(malaria$lac, malaria$enfermedad, ylab= "Probabilidad de enfermarse", xlab= "Conc. de anticuerpos (log)")
lines(seqx, valpred)