#
# Regresión Poisson
#
rm(list=ls())
help(warpbreaks)
data <- warpbreaks
#
# Método frecuentista
#
poisson.model <- glm(breaks ~ wool + tension, data, family = poisson(link = "log"))
summary(poisson.model)
#
# Predicción
#
nuevosdatos <- data.frame(wool=as.factor("A"),tension=as.factor("H"))
pred <- predict(poisson.model, newdata = nuevosdatos, type = "response")
pred
#
# Cambiar los datos al formato BUGS.
#
y <- data$breaks
n <- length(y)
mean(y)
woolB <- as.numeric(data$wool)-1
tensiont <- as.numeric(data$tension)-1
tensionM <- tensiont
tensionM[tensionM!=1] <- 0
tensionH <- tensiont
tensionH[tensionH!=2] <- 0
#
# Código OpenBUGS
#
if (!require("R2OpenBUGS")) install.packages("R2OpenBUGS")
library(R2OpenBUGS)
mlg <- function(){
  for (i in 1:n) {
    y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- intercepto+betaB*woolB[i]+betaM*tensionM[i]+betaH*tensionH[i]
  }
  intercepto ~ dnorm(0,0.0001)
  betaB ~ dnorm(0,0.0001)
  betaM ~ dnorm(0,0.0001)
  betaH ~ dnorm(0,0.0001)
}
#
# Definir los datos
#
mlgdata <- list(n = n,y = y,woolB = woolB,tensionM = tensionM,tensionH = tensionH)
#
# Definir valores iniciales
#
mlginits <- function() {
  list(intercepto=3.7,betaB=0,betaM=0,betaH=0)
}
#
# Correr el programa
#
mlgout <- bugs(data=mlgdata,inits=mlginits,parameters.to.save=c("intercepto","betaB","betaM","betaH"),model.file=mlg,
               n.chains = 1,n.iter = 10000)
mlgout
#
# Convergencia (debemos mirar los demás parámetros también)
#
intercepto <- mlgout$sims.list$intercepto
plot(intercepto,type='l',ylab='intercepto')
plot(cumsum(intercepto)/c(1:length(intercepto)),type='l',ylab='E[intercepto]')
acf(intercepto)
#
# Gráficos de las densidades a posteriori
#
plot(density(intercepto),ylab='f',xlab='intercepto',main="Densidad a posteriori del intercepto")
#
# Predicciones
#
betaH <- mlgout$sims.list$betaH
lambda <- exp(intercepto+betaH)
mean(lambda)
ypred <- rpois(length(lambda),lambda)
ynum <- table(ypred)
barplot(ynum/sum(ynum),xlab='y',ylab='P(y)')
quantile(ypred,c(0.025,0.5,0.975))