## raw growth data should be saved in tab separated
## with column names: c("id","growth","experim.1","experim.2", etc)
## datapub generated in read.data.r

## lme4 package needs to be installed.
## install.packages("lme4")

library(lme4)

## reshape package needs to be installed.
## install.packages("reshape")
library(reshape)

## read raw growth data

rawgrowth = read.table(file="rawgrowth.txt",sep='\t', header=T,
  quote = "\"", stringsAsFactors = F)

## as factors

rawgrowth$experim = as.factor(rawgrowth$experim)
rawgrowth$meas.by = as.factor(rawgrowth$meas.by)
rawgrowth$serum = as.factor(rawgrowth$serum)

## define reference
## 1 = MALE is reference for sex
## ASW is reference for population

rawgrowth$sex = as.factor(rawgrowth$sex)
rawgrowth$pop = as.factor(rawgrowth$pop)
rawgrowth$pop = relevel(rawgrowth$pop,"ASW")

## fit the mixed effects model

fitlmer = lmer(growth~pop+sex+serum+(1|meas.by)+(1|IID),data=rawgrowth)

## compute fixed effects

coefi = fixef(fitlmer)
xx = with(rawgrowth,model.matrix(~1+pop+sex+serum))
xx[,"serumFBS"] = 1
rawgrowth$xbeta= apply(xx %*% diag(coefi), 1, sum)

## reformat data
## this dataset has multiple entries per IID

rawgrowth= melt(rawgrowth,id=c("IID","pop","sex","experim","meas.by","serum"),measure.vars="xbeta")

## generate growth with only one entry per IID

growth = cast(rawgrowth,IID+pop+sex~.,mean, na.rm=T)
##vartemp = cast(rawgrowth,IID+pop+sex~.,var, na.rm=T) ## to check var=0
names(growth)[names(growth)=="(all)"] = "xbeta"

## compute igrowth0 (random effect for IID)

igrowth0 = ranef(fitlmer)$IID
igrowth0 =  cbind(rownames(igrowth0),igrowth0)
names(igrowth0) = c("IID","igrowth0")

growth = merge(growth, igrowth0, by.x="IID", by.y="IID")

## compute igrowth = igrowth0 + x*betas

growth$igrowth = growth$igrowth0 + growth$xbeta

growth[c("igrowth0","xbeta")] = list()

## write data

write.table(growth,file="igrowth.csv",sep=",",quote=F,col=T,row=F)