#############################################################################################
### Curso modelos lineales generales, generalizados y mixtos en ecologia IE-UNAM          ###
### Modelos lineales generalizados con distribucion Poisson y binomial negativa           ###
### Funciones de verosimilitud y metodos de optimizacion                                  ###
### Roberto Munguia-Steyer                                    2011/11/25                  ###
#############################################################################################


R

# Cargando paquetes
# install.packages("aod")

library(bbmle)
library(car)
library(MASS)
library(aod)
library(plotrix)

# Modelos sin variables predictoras

# Imaginemos que el valor esperado de insectos atrapados en una planta carnivora del genero Drosera son 3 insectos por dia (lambda= 3)
# a) ¿Cual es el es la probabilidad de que una planta haya atraapado solo 1 insecto?  b)  ¿y cinco?

# a)  
dpois(1, lambda=3)

# b)
dpois(5, lambda=3)


# Podemos calcular la probabilidad a mano, empleando la formula de la funcion de probabilidad de masa
# p(y|lambda) = ( e^(-lambda) * lambda^y )/ y!

# a) 

exp(-3)*(3^1)/factorial(1)


# b) 

exp(-3)*(3^5)/factorial(5)

# Calculemos la probabilidad de que la planta carnivora contenga de 0 a 15 insectos y grafiquemos las probabilidades

pr <- dpois(0:20, lambda=3)
barplot(pr)

# Los posibles resultados son exhaustivos y complementarios (suma de las probabilidades respectivas = 1)
sum(pr)

# Realicemos una simulación y consideremos que estudiamos 500 plantas carnivoras, 
# Grabamos y contabilizamos el numero de insectos que fueron capturados en cada una durante un lapso de 24 hrs.


plantcar <- rpois(500, lambda=3)

#Numero de plantas con 0 a n:insectos
freqin <- table(plantcar)
freqin

par(las)
barplot(freqin, xlab= "Frecuencia", ylab= "Numero de insectos en las plantas")

#Cual es el valor esperado
mean(plantcar)

#y la varianza
var(plantcar)

# Con la distribucion  Poisson el valor esperado y su varianza es la misma.


## Modelo de distribución Poisson sin y con variables predictoras
## Numero de frutos presentes en los arboles 

# Fijando el valor del generador de numeros pseudoaleatorios

set.seed(1234)

#Parametros

nplan<-1000
lambda=5
frutos<-rpois(nplan,lambda)
table(frutos)

#¿Cual es la probabilidad de que una planta sexualmente madura se encuentre sin frutos?

dpois(0,lambda)

# Comparelo con la proporcion de plantas sin frutos obtenidas a partir de nuestros datos simulados:

sum(frutos==0)/length(frutos)

# Vamos a comparar las distribucion de probabilidad teorica con las proporciones empiricas para todos los valores:

mfrut<-max(frutos)
fa<-factor(frutos, levels=0:mfrut)
prob.obs<-table(fa)/nplan
par(las=1)
plot(0:mfrut,prob.obs, xlab="Numero de frutos", ylab="Probabilidad", type="h", lwd=5)

#Agregamos los valores teoricos esperados:

prob.tr<-dpois(0:mfrut, lambda)
points(0:mfrut,prob.tr, pch=21, col="red")

#Definimos la funcion de maxima verosimilitud del parametro lambda:

x<-frutos
poisNLL<-function(lambda){
-sum(dpois(x, lambda, log=TRUE))
}

#Aplicamos la funcion sobre una serie de valores:

xvec<-seq(4.85,5.3, length=1000)
LLest<-sapply(xvec,poisNLL)

#Determinamos el valor de maxima verosimilitud o equivalentemente
# el valor minimo de la log(verosimilitud) negativa

#Por fuerza bruta
minLLest<-min(LLest)
lambdaLL.fb<-xvec[ LLest == min(LLest)]

#Por metodos de optimizacion usando la funcion mle2
lambdaLL.nm<-mle2(poisNLL, start=list(lambda=4))


# Comparemos las estimaciones de lambda obtenidas por fuerza bruta,
# por metodos numericos y analiticos 

lambdaLL.fb
lambdaLL.nm
mean(frutos)

#Graficamos la funcion de verosimilitud

par(las=1)
mfrutos<-mean(frutos)
LLest2<-LLest-min(LLest)
plot(xvec,LLest2, typ="l", xlab="frutos", ylab="loglik")
abline(v=mfrutos, col="blue", lwd = 3)
abline(v=coef(lambdaLL.nm),col ="darkgray")
abline(v=lambdaLL.fb, col="red")

# Si consideramos que la produccion de frutos depende de un nutriente, por ejemplo la cantidad de fosforo
# en el suelo, ojo las cantidades del nutriente son puramente especulativas


set.seed(1234)
phos<-runif(100,0,10)
a= 1
b= 0.3
x<-phos

# Los valores esperados de frutos en relacion al nivel de fosforo seguiran una relacion exponencial
# log(y) = a + b*x equivalente a y = exp(a + b*x)

ydet<-exp(a+b*x)
plot(x,ydet)

# Agreguemos el componente aleatorio al modelo:

fec<-rpois(100,ydet)

# Grafiquemos nuestros datos que simulan un proceso Poisson con un valores esperados de frutos dependiendo de la concentracion de fosforo
# en el suel o, que simulam um processo Poisson com valor esperado que é uma função do nível de fósforo:

par(las=1)
plot(phos, fec, xlab="Fosforo mg/Kg", ylab="Numero de frutos")

# Definimos la funcion de verosimilitud para este modelo:

poisglmNLL = function(a,b) {
   ypred= exp(a+b*x)
  -sum(dpois(fec,lambda=ypred, log=TRUE))
}


# Realizamos la optimizacion de la funcion de verosimilitud con el paquete mle2:

mod.pois<-mle2(poisglmNLL, start= list(a=2.5,b= 0.33))

# Veamos los parametros estimados:

summary(mod.pois)

# Los perfiles de verosimilitud y los intervalos de confianza:

plot(profile(mod.pois))
confint(mod.pois)

# Comparemos los resultados obtenidos con la funcion generica glm

mod.pois2 <- glm(fec~phos, poisson(link="log")) #No es necesario especificar la funcion de ligacion log, es la default.
coef(mod.pois);coef(mod.pois2)

# Validacion del modelo
par(mfrow=c(2,2))
plot(mod.pois2)

# A pesar de haber simulado los datos con una distribucion Poisson encontramos mayor varianza que la esperada

summary(mod.pois2) #Devianza residual/gl residuales

# Corrijamos los valores estimados para los errores estandar

mod.pois3 <- update(mod.pois2, family= quasipoisson)
summary(mod.pois3) #La sobredispersion no es muy alta, de 1.24

# Cuando realizamos un analisis de devianza tanto en modelos Poisson como binomiales y estos presentan sobredispersion
# la prueba estadística sera de F en vez de la  ji cuadrada anova(mod.pois3, test="Chi")

anova(mod.pois3, test="F") # En este caso los valores de probabilidad no diferirian sustancialmente por la magnitud del efecto


## Distribucion binomial negativa

# Otra posibilidad para lidiar con la sobredispersion es emplear la distribucion binomial negativa
# que adicionalmente presenta el parametro k, el cual mientras tenga un valor menor significara que existe una mayor
# sobredispersion, ks tendiendo al infinito no distinguiran a los modelos poisson y los de la binomial negativa


#Funcion de verosimilitud
negbinNLL<- function(a,b,k){
	ypred<-exp(a+b*x)
	-sum(dnbinom(fec, mu=ypred, size=k, log = TRUE))
        }

#Obtenemos un valor inicial aproximado de k:

med<-mean(x)
vari<-var(x)
k.init <-med^2/(vari-med)
k.init

mod.negbin <- mle2(negbinNLL, start=list(a=2.5, b= 0.33, k=k.init))

# Veamos los parametros estimados del modelo:

summary(mod.negbin)

# Equivalentemente podemos utilizar la funcion glm.nb que se encuentra en el paquete MASS:

mod.negbin2 <- glm.nb(fec ~ phos)
summary(mod.negbin2)

# Comparemos los coeficientes de ambos modelos

coef(mod.negbin); coef(mod.negbin2)

# Obtengamos los coeficientes estimados del modelo con dist binom neg
a.est2 <- coef(mod.negbin)[1]
b.est2 <- coef(mod.negbin)[2]

# Resulta necesario utilizar la dist. binomial negativa en vez de la dist Poisson:

AICtab(mod.pois,mod.negbin, delta=T, sort=T, weights = TRUE)

# Grafiquemos la curva de los valores esperados:

par(las=1)
plot(phos,fec, xlab="Fosforo mg/Kg", ylab="Numero de frutos" )
a.est<-coef(mod.pois)[1]
b.est<-coef(mod.pois)[2]
curve(exp(a + b*x),add=TRUE, col="red")
curve(exp(a.est + b.est*x), add=TRUE, col="blue", lty=2) # valores predichos Poisson
curve(exp(a.est2 + b.est2*x), add = TRUE, col = "purple", lwd=3, lty=3) # valores predichos binomial negativa
legend("topleft", c("Parametro simulado","Parametro estimado Pois", "Parametro estimado Binomial Negativa"),col=c("red","blue", "purple"), lty=c(1,2,3), lwd=c(1,1,3))


# Ejercicio 12.1 Existe una especie de hongo macrofago boletaceo que crece en diferentes tipo de vegetacion
# Realizamos conteos de los hongos presentes en 20 cuadros de 10 X 5 m, en los bosques de encino, de pino y oyamel
# Queremos saber si la abundancia del hongo difiere entre los tres tipos de vegetacion
# Ajusta y valida el modelo
# En caso de existir sobredispersion, a) utiliza la correccion para los modelos Poisson y determina de cuanto es la sobredispersion 
# b) ajusta un modelo con distribucion binomial negativa 
# De existir diferencias en la abundancia, calcula los valores esperados de hongos para cada tipo de vegetacion con 
# sus intervalos de confianza al 95 % y en base a ello indica si la abundancia difiere entre be, bp y bo con todas las combinaciones posibles.