# sample R code for mixed effects regression analysis

# first read source data files and recode as an R dataset
#
# Add path descriptions to file names if needed
y1=read.csv('runners1.csv',colClasses="character")
n1=length(y1[,1])/3
y2=matrix(nrow=n1,ncol=6)
name3=rep(NA,2000)
ind1=0
name1='uyx'
for(i in 1:n1){
name2=y1[3*i-2,3]
y2[i,6]=1
if(name2!=name1){
ind1=ind1+1
name3[ind1]=name2
name1=name2}
y2[i,1]=ind1
y2[i,2]=as.numeric(y1[3*i-2,1])
y2[i,3]=as.numeric(y1[3*i-2,4])
y2[i,4]=1
if(y1[3*i-2,5]=='M')y2[i,4]=2
tim=y1[3*i,5]
tim1=y1[3*i,6]
if(tim1!='')tim=tim1
hr=as.numeric(substr(tim,1,1))
min=as.numeric(substr(tim,3,4))
sec=as.numeric(substr(tim,6,7))
y2[i,5]=60*hr+min+sec/60
}
name3=name3[1:ind1]

# add 2013 DNFs
BM13b=read.csv('BM2013_DNF_List.csv',colClasses="character")
Projr=read.csv('ProjectedResults.csv',colClasses="character")
for(i1 in 1:ind1){
n3=nchar(name3[i1])
ind3=rep(0,n3)
for(k in 1:n3){if(substr(name3[i1],k,k)==',')ind3[k]=1}
k=min(which(ind3==1))
name4=paste(substr(name3[i1],k+2,n3),substr(name3[i1],1,k-1))
i2=which(BM13b[,2]==name4)
if(length(i2)>0){
i2=min(i2)
bibnum=BM13b[i2,1]
i3=which(Projr[,2]==bibnum)
if(length(i3)>0){
t3=60*as.numeric(Projr[i3,7])+as.numeric(Projr[i3,8])+as.numeric(Projr[i3,9])/60
i4=max(which(y2[,1]==i1))
if(i4<n1){y2=rbind(y2[1:i4,],y2[i4,],y2[(i4+1):n1,])}
if(i4==n1){y2=rbind(y2[1:n1,],y2[n1,])}
n1=n1+1
y2[i4+1,2]=2013
y2[i4+1,3]=as.numeric(Projr[i3,3])
y2[i4+1,5]=t3
y2[i4+1,6]=0
}
}
}

# y2 is the 7219 x 6 data file for future analysis
#
# 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

# First load splines and lme4 packages - not that lme4 must be installed before it can be loaded

library(splines)
library(lme4)

# form spline basis functions
df=6
agemin=min(y2[,3])
agemax=max(y2[,3])
agesp=ns(agemin:agemax,df)
age30=30-agemin+1
age0=agesp[age30,]
for(j in 1:(agemax-agemin+1)){agesp[j,]=agesp[j,]-age0}
sv1=svd(agesp)
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])

# basic fitting routine - individual effect and year effect as random effects
lm2=lmer(formula=y3~x1+(1|ind2)+(1|year))

# form matrix of output and plot
out1=matrix(nrow=agemax-agemin+1,ncol=4)
out1[,1]=agemin:agemax
out1[,2]=agesp%*%fixef(lm2)[1+1:df]
lmvar=agesp %*% as.matrix(vcov(lm2)[1+1:df,1+1:df]) %*% t(agesp)
lmsd=sqrt(diag(lmvar))
out1[,3]=out1[,2]-2*lmsd
out1[,4]=out1[,2]+2*lmsd
plot(out1[,1],out1[,2],ylim=range(out1[,3:4]),type='n',xlab='Age',ylab='Effect of Age')
lines(out1[,1],out1[,2],lty=1)
lines(out1[,1],out1[,3],lty=2)
lines(out1[,1],out1[,4],lty=2)