### Estimating COMPONENTS Of VARIANCE, and INTRA-CLASS CORRELATIONS

### Example  2 measurements of each of  8 subjects by 4 raters 

setwd("/Users/jameshanley/Documents/Courses/bios602/MultilevelModels/BigEars")

earlength = c(67,65,65,64, 74,74,74,72, 67,68,66,65, 65,65,65,65, 65,62,62,61, 59,56,55,53, 60,62,60,59, 66,65,65,63,
              67,66,66,66, 74,73,71,73, 68,67,68,67, 64,65,65,64, 61,62,60,61, 59,56,55,53, 60,62,60,58, 66,65,65,65)

summary(earlength)

subject =rep(rep(1:8,each=4),2) ; subject
observer=rep(1:4,16) ; observer

dump(c("subject", "observer", "earlength"),"earsData.txt")

results=read.coda("CODAchain1.txt","CODAindex.txt") 

summary(results)

summary( aov(earlength~as.factor(subject)+as.factor(observer) ) )

summary( aov(earlength~as.factor(subject)+as.factor(observer) + as.factor(subject)*as.factor(observer)  ) ) 




####### WinBUGS Directly

model
{
  for( i in 1:64) {
       earlength[i] ~ dnorm(mu[i],inv.var.error)
       mu[i] <- mu.overall + alpha[ subject[i] ] + beta[ observer[i] ]
  }     
  mu.overall ~ dnorm(0, 0.001)
 
  for( s in 1:8 ) {
      alpha[s] ~ dnorm(0, inv.var.subjects)
  }
 for( o in 1:4 ) {
      beta[o] ~ dnorm(0, inv.var.observers)
 }
   sigma.subjects  ~ dunif(0,1000)
   sigma.observers ~ dunif(0,1000)
   sigma.error     ~ dunif(0,1000)

   var.subjects  <- sigma.subjects   * sigma.subjects
   var.observers <- sigma.observers  * sigma.observers
   var.error     <- sigma.error      * sigma.error

 icc   <- var.subjects / (var.subjects + var.observers + var.error)

   inv.var.subjects  <- 1/var.subjects  
   inv.var.observers <- 1/var.observers
   inv.var.error     <- 1/var.error
}

#data
list(earlength = c(67,65,65,64, 74,74,74,72, 67,68,66,65, 65,65,65,65, 65,62,62,61, 59,56,55,53, 60,62,60,59, 66,65,65,63,
              67,66,66,66, 74,73,71,73, 68,67,68,67, 64,65,65,64, 61,62,60,61, 59,56,55,53, 60,62,60,58, 66,65,65,65),
     subject=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,
               1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8),
     observer=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,
                1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4))

#inits
list(alpha = c(-5,-5,-5,-5, 5,5,5,5), beta= c(-5,-5, 5,5), sigma.subjects =  5, 
sigma.observers = 5 , sigma.error = 5 , mu.overall=  60)