# plot cholestrole values over time for the paper

postscript(file="cholst2.ps", horizontal = F, 
width=8.5, height=6.5)
par(mfrow=c(1,1))

out <- read.table(file="cholst.dat", col.names=c("newid", "id", "cholst",
 "sex", "age",  "year"))

plot(out$year, out$cholst, xlab="Time in years", ylab="Cholesterol level",
cex=0.8)

b0 <-  138.18 
b1 <- -9.6393 
b2 <- 2.0509  
b3 <- 6.8003  
b4 <- 1.7995  
b5 <- -0.1145 

x0 <- 10.8

title("Cholesterol levels over time", cex=0.8)
lines(out$year[out$newid==2], out$cholst[out$newid==2], lty=2)
int <- 100.50
slp <- 2.7414 
sex <- out$sex[out$newid==2]
age <- out$age[out$newid==2]
year <- out$year[out$newid==2]
y <- b0 + int + b1*sex + b2*age + (b3 + b4*sex + b5*age + slp)*year
lines(year, y, lty=1)
text(x0, y[length(y)], "ID=2", cex=0.7)

lines(out$year[out$newid==74], out$cholst[out$newid==74], lty=2)
int <- 46.9844 
slp <- 1.3579  
sex <- out$sex[out$newid==74]
age <- out$age[out$newid==74]
year <- out$year[out$newid==74]
y <- b0 + int + b1*sex + b2*age + (b3 + b4*sex + b5*age + slp)*year
lines(year, y, lty=1)
text(x0, y[length(y)], "ID=74", cex=0.7)

lines(out$year[out$newid==171], out$cholst[out$newid==171], lty=2)	
int <- -51.5764 
slp <-  -0.6812 
sex <- out$sex[out$newid==171]
age <- out$age[out$newid==171]
year <- out$year[out$newid==171]
y <- b0 + int + b1*sex + b2*age + (b3 + b4*sex + b5*age + slp)*year
lines(year, y, lty=1)
text(x0, y[length(y)], "ID=171", cex=0.7)

dev.off()