x  = seq(1,100)/100
y1 = 0.1*sin(2*pi*x)+0.25+(x-0.25)/3
y2 = y1
for(i in 26:100){
   y2[i] = y1[i]+(x[i]-0.25)*2/3-(x[i]-0.25)^2
}

pdf("fig41.pdf",width=5,height=5.2)
par(mar=c(2,2.1,1.2,2))

plot(x,y1,type="l",lty=1,lwd=2,xlab="",
     ylab="Disease Risk Factor",mgp=c(1,1,0),
     xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n",cex=1.4,
     cex.axis=1.4,cex.lab=1.4)
lines(x[26:75],y2[26:75],lty=2,lwd=3,cex=1.4)
lines(c(x[75],x[75]),c(y1[75],y2[75]),lty=1,lwd=4,cex=1.4,col="black")
axis(1,at=c(0.75),labels=c("Current Time"),cex.axis=1.4)
legend("topleft",legend=c("Regular Longitudinal Pattern",
       "Longitudinal Pattern of a Patient"),
       lty=1:2,lwd=3,cex=1.2)
box(lty=1, lwd=4)
polygon(c(x[25:75],rev(x[25:75])),c(y1[25:75],rev(y2[25:75])),
     border=par("bg"),col="skyblue")

graphics.off()