# 2012.11.05 ### 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) ######################################################################################### R JAGS additive.model=readLines(,28) 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 } writeLines(additive.model,"additiveModel.txt") 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) library('rjags') jags <- jags.model("additiveModel.txt", data = DATA, inits=list( mu.overall=60, sigma.subjects=4, sigma.observers=2, sigma.error=1), n.chains = 1, n.adapt = 5000) update(jags, 10000) results=jags.samples(jags, c('mu.overall','var.subjects','var.observers', 'var.error', 'icc'), 10000) str(results) summary( results$icc[1,,1] ) summary(results$mu.overall[1,,1]) summary(results$var.subjects[1,,1]) summary(results$var.observers[1,,1]) summary(results$var.error[1,,1]) summary(results$icc[1,,1]) apply(results$icc,1,summary) hist(results$icc[1,,1],breaks=seq(0,1,0.02)) round( ( sort(results$icc[1,,1]) ) [10000*c(0.025,0.975)] , 2)