rm(list=ls())

#### Bone density example for QUANTILES

library(quantreg)


#### Read data

df = read.table("bone.txt",header=TRUE)
attach(df)
n = nrow(df)
age = df$age
BoneD = df$spnbmd


### NONPARAMETRIC QUANTILE Regression.

### function lprq(...m...) returns a grid xx of m points
###   and fitted values fv on m points

## h=2

pdf("BoneQuant_h=2.pdf")
plot(age,BoneD,pch=21,cex=1.8,bg="lightpink",ylab="",
            xlab="",yaxt="n",cex.axis=1.7)
out = lprq(age,BoneD,h=2,m=50)        ## get fit (fv) on grid of size m
lines(out$xx,out$fv,lwd=5,col="purple")
out = lprq(age,BoneD,h=2,m=50,tau=.1) 
lines(out$xx,out$fv,col="red",lwd=5)
out = lprq(age,BoneD,h=2,m=50,tau=.9) 
lines(out$xx,out$fv,col="red",lwd=5)
title(xlab="Age",line=1.8,cex.lab=2)
title(xlab="h=2",line=3.2,cex.lab=2)
title(ylab="Change in BMD",line=2.3,cex.lab=2)
dev.off()


##  h=4

pdf("BoneQuant_h=4.pdf")
plot(age,BoneD,pch=21,cex=1.8,bg="lightpink",ylab="",
            xlab="",yaxt="n",cex.axis=1.7)
out = lprq(age,BoneD,h=4,m=50)        ## get fit (fv) on grid of size m
lines(out$xx,out$fv,lwd=5,col="purple")
out = lprq(age,BoneD,h=4,m=50,tau=.1) 
lines(out$xx,out$fv,col="red",lwd=5)
out = lprq(age,BoneD,h=4,m=50,tau=.9) 
lines(out$xx,out$fv,col="red",lwd=5)
title(xlab="Age",line=1.8,cex.lab=2)
title(xlab="h=4",line=3.4,cex.lab=2)
title(ylab="Change in BMD",line=2.3,cex.lab=2)
title(xlab="h=4",line=3,cex.lab=2)dev.off()

## NONLINEAR  h=2

pdf("BoneQuant_NonL.pdf")
h = 2
plot(age,BoneD,pch=21,cex=1.8,bg="lightpink",ylab="",  
           xlab="",yaxt="n",cex.axis=1.7)
out = lprq(age,BoneD,h=h,m=50,tau=.5)    
lines(out$xx,out$fv,lwd=5,col="purple")  
out = lprq(age,BoneD,h=h,m=50,tau=.1)    
lines(out$xx,out$fv,lwd=5,col="red") 
out = lprq(age,BoneD,h=h,m=50,tau=.9)    
lines(out$xx,out$fv,lwd=5,col="red")     
title(xlab="Age",line=1.8,cex.lab=2)
title(xlab="h=4",line=3.4,cex.lab=2)
title(ylab="Change in BMD",line=2.3,cex.lab=2)
dev.off()