y2=read.table('C:/Users/rsmith/jan16/research/BAA/AgeCurves/BostonRunners.txt')
#
# Each row of y2 corresponds to one performance by one of the 806 runners in the database
#
# Columns are:
#
# Runner number (1 through 806)
# Year of race
# Age on day of race
# Sex (1=female, 2=male)
# Finish time in minutes
# Indicator of whether time was known (=1) or estimated (=0; 2013 finishers only)

# Following is a sample analysis using the lme4 package
#
# The age-time curve is modeled through splines with 6 degrees of freedom

library(splines)
library(lme4)
# form spline basis functions
df=6
agemin=min(y2[,3])
agemax=max(y2[,3])
agesp=ns(agemin:agemax,df)
x1=matrix(nrow=length(y2[,1]),ncol=df)
for(i in 1:length(y2[,1])){
j=y2[i,3]-agemin+1
x1[i,]=agesp[j,]
}
y3=log(y2[,5])
year=as.factor(y2[,2])
ind2=as.factor(y2[,1])
sex=as.factor(y2[,4])

# basic fitting routine - individual effect and year effect as random effects
# females first then males
agesp0=cbind(rep(1,62),agesp)
lm2=lmer(formula=y3~x1+(1|ind2)+(1|year),subset=(sex==1))
out1=matrix(nrow=agemax-agemin+1,ncol=4)
out1[,1]=agemin:agemax
out1[,2]=agesp0%*%fixef(lm2)[1:(df+1)]
lmvar=agesp0 %*% as.matrix(vcov(lm2)[1:(df+1),1:(df+1)]) %*% t(agesp0)
lmsd=sqrt(diag(lmvar))
out1[,3]=out1[,2]-2*lmsd
out1[,4]=out1[,2]+2*lmsd
lm2=lmer(formula=y3~x1+(1|ind2)+(1|year),subset=(sex==2))
out2=matrix(nrow=agemax-agemin+1,ncol=4)
out2[,1]=agemin:agemax
out2[,2]=agesp0%*%fixef(lm2)[1:(df+1)]
lmvar=agesp0 %*% as.matrix(vcov(lm2)[1:(df+1),1:(df+1)]) %*% t(agesp0)
lmsd=sqrt(diag(lmvar))
out2[,3]=out2[,2]-2*lmsd
out2[,4]=out2[,2]+2*lmsd
out1[,2:4]=exp(out1[,2:4])
out2[,2:4]=exp(out2[,2:4])

# draw figure
par(mfrow=c(1,1),cex=1.3)
plot(out1[,1],out1[,2],ylim=range(c(out1[,3:4],out2[,3:4])),type='n',xlab='Age',
ylab='Mean Time (Minutes)')
lines(out1[,1],out1[,2],lty=1,lw=5,col='blue')
lines(out1[,1],out1[,3],lty=2,lw=3,col='blue')
lines(out1[,1],out1[,4],lty=2,lw=3,col='blue')
lines(out2[,1],out2[,2],lty=1,lw=5,col='red')
lines(out2[,1],out2[,3],lty=2,lw=3,col='red')
lines(out2[,1],out2[,4],lty=2,lw=3,col='red')