# 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.

## Generando la base de datos
set.seed(1234)

habitat <- gl(3,20, labels= c("be", "bp", "bo"))

#log(8)
a <- log(8)

# log(12)
b1 <- log(11) - log(8)

b2 <- log(12) - log(8)

mm <- model.matrix(~habitat)
betavec<- c(a,b1,b2)

pl <- mm[,] %*% betavec

mu <- exp(pl)
k <- 5

y <- rnbinom(60,mu=mu, size= k)

boletus <- data.frame(y, habitat)
names(boletus) <- c("abundancia", "habitat")

help(write.table)
write.table(boletus, "boletus.csv", quote=FALSE, sep=",", row.names=FALSE, col.names=c("abundancia", "habitat"))

## Ajustando el modelo
boletus <- read.table("boletus.csv", sep=",", header=TRUE)


m4 <- glm(abundancia~habitat, poisson, data=boletus)
exp(coef(m4))

par(mfrow=c(2,2))
plot(m4)

summary(m4)

m5 <- update(m4, family=quasipoisson)
summary(m5)

par(mfrow=c(2,2))
plot(m5)

anova(m5, test="Chi")
anova(m5, test="F")

library(effects)

#Difiere bp y bo?
efectos <- allEffects(m5)
summary(efectos)

#con la funcion predict
m6 <- glm(abundancia~habitat-1, quasipoisson, data= boletus)
summary(m6)
exp(coef(m6))

m7 <- glm.nb(y~habitat, boletus)
anova(m7)

m8 <- glm.nb(abundancia~habitat-1, boletus)
confint(m8)
m <- exp(confint(m8))
m


# Valores predichos e intervalos de confianza al 95 %
ii <- m[,1]
is <- m[,2]
vp <- exp(coef(m8))

#Graficando
par(las=1)
plotCI(x= 1:3,y=vp, ui=is, li=ii, names=c("be","bp","bo"), ylab="Abundancia Boletus")