# 2012.11.05

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

### Example" weight gains, twin A twin B (1st pair)... twin A twin B (12th pair)
### Bouchard C. NEJM 

# in R

x=  c(13.3,11, 11.1,12.9, 8.2,9.6, 6.1,11.4, 7.9,7.7, 7.1,6.9,  6.9,8, 6.5,8, 6.7,8.8,  7.3,4.3, 6.5,5.3, 5.4,7.1)
	  
m=matrix(x,ncol=12,byrow=F)

## METHOD OF MOMENTS 

var(x) ;mean(x)

family.means=apply(m,2,mean) ; family.means ; 
var(family.means)

between.ss = 2*11*var(family.means) ; between.ss

mean.square.between = between.ss/11 ; mean.square.between

within.family.ss = 1*apply(m,2,var) ;  within.family.ss
pooled.num = sum(within.family.ss);
est.of.within.family.sigma.squared = pooled.num/12
est.of.within.family.sigma.squared

est.of.between.family.sigma.squared = (mean.square.between - est.of.within.family.sigma.squared)/2
est.of.between.family.sigma.squared

est.of.icc = est.of.between.family.sigma.squared/
            (est.of.between.family.sigma.squared + est.of.within.family.sigma.squared)

round(est.of.icc,2)


####################################################### WinBUGS 

model
{ 
   for( i in 1:12 ) {
	mu[i ] ~ dnorm(population.mu, tau.b)
      for( j in 1:2 ) {
         Y[i , j] ~ dnorm(mu[i ], tau.w)
	}
    }		
  # priors
   population.mu ~ dnorm(0, 0.001)
   sigma.b ~ dunif(0,100)
   sigma.w ~ dunif(0,100)

   tau.b <- 1/ (sigma.b * sigma.b)
   tau.w <- 1/ (sigma.w * sigma.w)

   var.b <- 1/tau.b
   var.w <- 1/tau.w
   icc <- ( var.b ) / ( var.b + var.w ) 
}

#data
list(Y = structure(
.Data =  c(13.3,11, 11.1,12.9, 8.2,9.6, 6.1,11.4, 7.9,7.7, 7.1,6.9,  6.9,8, 6.5,8, 6.7,8.8, 7.3,4.3, 6.5,5.3, 5.4,7.1),
.Dim = c(12,2)))

#inits
list(sigma.b=1, sigma.w=1)

list(sigma.b=3 sigma.w=1)

list(sigma.b=1, sigma.w=3)



################################################################### R JAGS

#######  # transpose data for JAGS (opposite to WinBUGS)

MODEL=readLines(,20) 
model
{ 
   for( i in 1:12 ) {
	mu[i ] ~ dnorm(population.mu, tau.b)
      for( j in 1:2 ) {
         Y[j,i] ~ dnorm(mu[i ], tau.w)
	}
    }		
  # priors
   population.mu ~ dnorm(0, 0.001)
   sigma.b ~ dunif(0,100)
   sigma.w ~ dunif(0,100)

   tau.b <- 1/ (sigma.b * sigma.b)
   tau.w <- 1/ (sigma.w * sigma.w)

   var.b <- 1/tau.b
   var.w <- 1/tau.w
   icc <- ( var.b ) / ( var.b + var.w ) 
}

writeLines(MODEL,'ICCModel.txt')


# transpose data for JAGS (opposite to WinBUGS)

d=matrix(c(13.3,11, 11.1,12.9, 8.2,9.6, 6.1,11.4, 7.9,7.7, 7.1,6.9,  
          6.9,8, 6.5,8, 6.7,8.8, 7.3,4.3, 6.5,5.3, 5.4,7.1),
          2,12)

paste(t(d),collapse=",")

DATA=list(Y = structure(
.Data =  c(13.3, 11.1, 8.2,  6.1, 7.9, 7.1, 6.9, 6.5, 6.7, 7.3, 6.5, 5.4,
          11,   12.9, 9.6, 11.4, 7.7, 6.9, 8,   8,   8.8, 4.3, 5.3, 7.1),
.Dim = c(2,12)))


library('rjags')

jags <- jags.model('ICCModel.txt',
                   data = DATA,
                   inits=list(sigma.b=3, sigma.w=1),
                   n.chains = 1,
                   n.adapt = 5000)
 
update(jags, 10000)


results=jags.samples(jags, 
  c('var.b','var.w','icc'), 10000)
str(results)

summary(results$icc[1,,1])
summary(results$var.b[1,,1])
summary(results$var.w[1,,1]) 

hist(results$icc[1,,1],breaks=seq(0,1,0.02)) 

round( ( sort(results$icc[1,,1]) ) [10000*c(0.025,0.975)] , 2)