setwd("C:/Users/Lisa/Dropbox/CSUMB stats class")

library("BRugs") 

#Data
# Column 1: Species
#1 = guinea pig
#2 = rat
#3 = baboon
#4 = human
#5 = horse
#6 = rabbit

# Column 2: H2O_ratio = % H2O/% not H2O
# Column 3: logit_Pfat = ln(% fat/% not fat)

waterfatdata <- read.csv(file = "HW8data.csv", header = TRUE)

# Put data in to OpenBugs format.
data.in <- list(Species = waterfatdata$Species,
                logit_Pfat = waterfatdata$logit_Pfat,
                H2O_ratio = waterfatdata$H2O_ratio,
                N = length(waterfatdata$Species),
                M = 6)
bugsData(data = data.in, fileName = "waterfatdata.txt")


# Model

model{
# Likelihood.
  for(i in 1:N){  # March through each animal (N = total sample size).
    SCode <- Species[i]
    logit_Pfat[i] ~ dnorm(mu[i], tau)
    mu[i] <- Beta0[SCode] + (Beta1[SCode]*H2O_ratio)
  }
# Priors on Beta0 and Beta1.
# This is hierarchical, so the parameters for each species come from an
#   overall distribution with hyper-parameters.
  for(i in 1:M){  # March through each species (M = total number of species).
    Beta0[i] ~ dnorm(Beta0.mu, Beta0.tau)
    Beta1[i] ~ dnorm(Beta1.mu, Beta1.tau)
  }
# Vague hyperpriors.
  Beta0.mu ~ dnorm(0,0.00001)
  Beta0.tau ~ dunif(0,100)
  Beta1.mu ~ dnorm(0,0.00001)
  Beta1.tau ~ dunif(0,100)
# Don't forget about the prior on error variance (not hierarchical).
  tau ~ dunif(0,100)
}

# Initial values.
inits.func <- function (){
  list (Beta0 = rnorm(n = 6,mean = c(0.91,
                                      -0.16,
                                      4.77,
                                      0.01,
                                      0.87,
                                      1.16), sd = 0.6),
        Beta1 = rnorm(n = 6,mean = c(-1.71,
                                      -0.93,
                                      -3.6,
                                      -1.0,
                                      -1.63,
                                      -1.65),sd = 0.4),
        tau = 1/(rnorm(1,10,10)^2))
}

bugsInits(inits = inits.func, numChains = 3, fileName = c("hier1.txt",
                                                          "hier2.txt",
                                                          "hier3.txt"),
          format="E")