#
# Selección de modelos.
#
rm(list=ls())
#
# Ejemplo: regresión Poisson.
#
help(warpbreaks)
data <- warpbreaks
#
# Método frecuentista.
#
poisson.model <- glm(breaks ~ wool + tension, data, family = poisson(link = "log"))
summary(poisson.model)
poissonw <- glm(breaks ~ wool, data, family = poisson(link = "log"))
summary(poissonw)
poissont <- glm(breaks ~ tension, data, family = poisson(link = "log"))
summary(poissont)
#
# Con MCMCpack.  Observamos que es importante usar a prioris propias.  Si no, la verosimilitud marginal no está definida.
#
if (!require("MCMCpack")) install.packages("MCMCpack")
library(MCMCpack)
bpoisson <- MCMCpoisson(breaks ~ wool + tension, data,burnin=10000*30,mcmc=10000*30,thin=30,
                        B0=0.0001,marginal.likelihood="Laplace")
bpoisson
bpoissonw <- MCMCpoisson(breaks ~ wool, data,burnin=10000*30,mcmc=10000*30,thin=30,
                         B0=0.0001,marginal.likelihood="Laplace")
bpoissonw 
bpoissont <- MCMCpoisson(breaks ~ wool, data,burnin=10000*30,mcmc=10000*30,thin=30,
                         B0=0.0001,marginal.likelihood="Laplace")
bpoissont
#
# Usamos OpenBUGS para calcular el DIC.
#
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
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
#
# Modelo con sólo wool
#
mlgw <- function(){
  for (i in 1:n) {
    y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- intercepto+betaB*woolB[i]
  }
  intercepto ~ dnorm(0,0.0001)
  betaB ~ dnorm(0,0.0001)
}

mlgdataw <- list(n = n,y = y,woolB = woolB)
mlginitsw <- function() {
  list(intercepto=3.7,betaB=0)
}
mlgoutw <- bugs(data=mlgdataw,inits=mlginitsw,parameters.to.save=c("intercepto","betaB"),model.file=mlgw,
                n.chains = 1,n.iter = 10000)
mlgoutw
#
# Modelo con sólo tensión
#
mlgt <- function(){
  for (i in 1:n) {
    y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- intercepto+betaM*tensionM[i]+betaH*tensionH[i]
  }
  intercepto ~ dnorm(0,0.0001)
  betaM ~ dnorm(0,0.0001)
  betaH ~ dnorm(0,0.0001)
}
#
# Definir los datos
#
mlgdatat <- list(n = n,y = y,tensionM = tensionM,tensionH = tensionH)
#
# Definir valores iniciales.
#
mlginitst <- function() {
  list(intercepto=3.7,betaM=0,betaH=0)
}
#
# Correr el programa.
#
mlgoutt <- bugs(data=mlgdatat,inits=mlginitst,parameters.to.save=c("intercepto","betaM","betaH"),model.file=mlgt,
                n.chains = 1,n.iter = 10000)
mlgoutt