# 2012.11.05

setwd("/Users/jameshanley/Dropbox/Courses/bios602/MultilevelModels/GaltonFamilyData")

ds=read.csv("galtonheights197.txt",sep=",",header=T,as.is=rep(T,23))
head(ds); tail(ds)

n=length(ds$Family)
height=c()
family.no=c()
i.male=c()
father=c(); mother=c()
f.no=0
for (i in 1:n){
	if( ds[i,2] != "* ") { 
		f.no = f.no + 1 ;  print(f.no) 
		fa = as.numeric(ds$Father[i] ) + 60;
		mo = as.numeric(ds$Mother[i] ) + 60;
	    for (j in 4:22){
		   x=ds[i,j] ; x=as.numeric(x)+60; print(x)
		   if(!is.na(x)) { family.no=c(family.no,f.no); i.male=c(i.male,1*(j<=13) ); 
		   	               father=c(father, fa );
		   	               mother=c(mother, mo ); 
		   	               height=c(height,x) ;
		   	               }
		 }
	}
}

family.no = family.no - 1*( family.no >162 )

DS=cbind(family.no,father,mother,i.male,height)


head(DS); tail(DS)
write.csv(DS,"heightsTallFile.csv")

############## can start here for ############################################  R JAGS

ds = read.csv("heightsTallFile.csv",header=T)

family.no = ds$family.no
father    = ds$father
mother    = ds$mother
i.male    = ds$i.male
height    = ds$height



FamilyHeightModel=readLines(,30)
model
{
  #data
  for( i in 1:899) {
       height[i] ~ dnorm(mu[i],inv.var.within)
       mu[i]    <- mu.female + extra.if.male * i.male[i]  + alpha[ family.no[i]  ] 
  }     
  
  #priors
  mu.female     ~ dnorm(0, 0.001)
  extra.if.male ~ dnorm(0, 0.001)

  # random.effects
  for( f in 1:197 ) {
      alpha[f] ~ dnorm(0, inv.var.family)
  }
 
  #priors
   sigma.family  ~ dunif(0,1000)
   sigma.within   ~ dunif(0,1000)

  #structural
   var.family      <- sigma.family   * sigma.family
   var.within      <- sigma.within   * sigma.within
   
   inv.var.family <- 1/var.family
   inv.var.within <- 1/var.within

   icc             <- var.family / (var.family + var.within)
}

writeLines(FamilyHeightModel,'FamilyHeightModel.txt')


#inits
INITS=list( mu.female =  63,
      extra.if.male   =   5,
      alpha = rnorm(197), 
      sigma.family =  1, sigma.within =  1)
 
library('rjags')

jags <- jags.model("FamilyHeightModel.txt",
                   data = list(i.male=i.male, family.no=family.no, height=height),
                   inits=INITS,
                   n.chains = 1,
                   n.adapt = 4000)
 
update(jags, 10000)
 
results=jags.samples(jags, c('alpha','mu.female','extra.if.male', 'var.family', 'var.within', 'icc'), 10000)
str(results)

summary(results$alpha[1,,1])
summary(results$alpha[197,,1])

summary(results$mu.female[1,,1])
summary(results$extra.if.male[1,,1]) 
summary(results$var.family[1,,1])
summary(results$var.within[1,,1])
summary( results$icc[1,,1] ) 

hist( results$icc[1,,1] )
hist( results$mu.female[1,,1] )
hist( results$extra.if.male[1,,1] )