#
# ANOVA bayesiana: y[i,j] = mu + theta[j] + epsilon[i,j] donde supongamos que theta[1] = 0
#
# i = 1, ..., n_j datos en grupo j
# j = 1, ..., k grupos
#
rm(list=ls())
#
data(PlantGrowth)
help(PlantGrowth)
PlantGrowth$group # group es un factor en R.  
PlantGrowth$weight # niveles de crecimiento de plantas bajo un control y dos tratamientos.
#
# Mirar los datos
#
boxplot(weight ~ group, data=PlantGrowth)
#
# Modelo frecuentista
# El efecto del grupo control es fijado en 0 por el modelo
#
summary(lm(weight ~ group, data=PlantGrowth)) # No hay diferencias significativas en media de los demás grupos y el control
#
# Haciendolo a mano.
#
y <- PlantGrowth$weight
X <- matrix(rep(0,90),ncol=3)
X[,1] <- 1
X[11:20,2] <- 1
X[21:30,3] <- 1
X
n <- dim(X)[1]
XTX <- t(X)%*%X
XTy <- t(X)%*%y
freqpars <- solve(XTX)%*%XTy
#
# Haciendolo con MCMCpack con una a priori poca informativa.  
#
if (!require("MCMCpack")) install.packages("MCMCpack")
library(MCMCpack)
resbayes <- MCMCregress(weight ~ group, data = PlantGrowth)
summary(resbayes)
intercepto <- resbayes[,1] 
gtrt1 <- resbayes[,2]
gtrt2 <- resbayes[,3]
#
# Convergencia
#
ts.plot(intercepto)
plot(cumsum(intercepto)/c(1:length(intercepto)),type="l",ylab="E[intercepto]",xlab="iters")
acf(intercepto)
#
# A posterioris
#
plot(density(intercepto),xlab="intercepto",ylab="f",main="")
#
# Haciendolo con R2OpenBUGS.  
#
# Primero tenemos que definir las variables.  
#
weight <- PlantGrowth$weight
xtr1 <- X[,2] # un vector que tiene el valor 1 para los miembros del grupo de tratamiento 1 y 0 para los demás
xtr2 <- X[,3] # un vector que tiene el valor 1 para los miembros del grupo de tratamiento 2 y 0 para los demás
#
# El modelo es weight[i,j] = beta0 + betatr1* xtr1[i]+betatr2*xtr2[i]+epsilon[i,j]
#
# Definir el modelo para OpenBUGS
#
modelobayes <- function(){
  for( i in 1 : n ) {
    weight[i] ~ dnorm(mu[i], phi)
    mu[i] <- beta0 + betatr1 * xtr1[i] + betatr2*xtr2[i]
  }
  beta0 ~ dnorm(0.0, 1.0E-6)
  betatr1 ~ dnorm(0.0, 1.0E-6)
  betatr2 ~ dnorm(0.0, 1.0E-6)
  phi ~ dgamma(0.001, 0.001)
}
#
# Definir los datos como list
#
datos <- list(n=length(weight),weight=weight,xtr1=xtr1,xtr2=xtr2)
#
# Definir los valores iniciales como list en una función
#
mean(weight)
iniciales <- function(){
  list(beta0=5,betatr1=0,betatr2=0,phi=1)
}
#
# Correr el programa con bugs
#
if (!require("R2OpenBUGS")) install.packages("R2OpenBUGS")
library(R2OpenBUGS)
resbayesob <- bugs(data = datos, inits = iniciales, parameters.to.save = c("beta0", "betatr1", "betatr2","phi"), 
                   model.file = modelobayes, n.chains = 1, n.iter = 20000)
#
# Resultados
#
resbayesob
#
# Análisis de convergencia
#
beta0 <- resbayesob$sims.list$beta0
betatr1 <- resbayesob$sims.list$betatr1
betatr2 <- resbayesob$sims.list$betatr2
phi <- resbayesob$sims.list$phi
ts.plot(beta0)
plot(cumsum(beta0)/c(1:length(beta0)),type='l')
acf(beta0)
#
# Mirando diferencias entre efectos medios en cada grupo
#
plot(density(beta0)) # distribución a posteriori del efecto media en el grupo de control.
plot(density(betatr1)) # diferencia entre la media en el grupo tr1 y el grupo de control
plot(density(betatr2)) # diferencia entre la media en el grupo tr2 y el grupo de control
plot(density(betatr2-betatr1)) # diferencia entre la media en el grupo tr2 y el grupo tr1
quantile(betatr2-betatr1,c(0.025,0.5,0.975))
sum(betatr2-betatr1>0)/length(betatr1) # hay una probabilidad de 99% de que la media en tr2 es más que la media en tr1.