# source( "CLLEO_QR_ReadingScores_IX_exposedXexcelhealth.Rscript" )

# R code Script to generate quantile regression plots as in Magzamen, Amato, et al., 2015,  Environmental Research, 117, 108-119. 
## script by Michael Amato, http://mikeamatowi.wix.com/profile
## Script is posted ot:   http://psych.wisc.edu/moore/R_Analysis_of_Variance_Handouts_html.html

# Uses the quantreg package, created by Roger Koenker, William B. McKinley Professor of Economics and Professor of Statistics, University of Illinois at Urbana-Champaign

#	http://cran.r-project.org/web/packages/quantreg/quantreg.pdf , retrieved 7/31/2014

# TESTS FOR INTERACTIONS

#################################################################

## Read in the data, require the package
	d1 = read.csv("CLLEO_qr.csv", header=T)
	require(quantreg)


## Specify the model
## Although specifying that the model fit an intercept with "+1" to the immediate right of the "~" is unnecessary, it helps when counting terms to specify which parameters to plot
	m4 = formula( RDG_KCE_SCALE_SCORE ~ 1 + exposed + child_sex + black + excelhealth + lths + ell + DPIpov + exposed:excelhealth )
	
	
## Specify some human-friendly titles for the plots
## Clearly it is EXTREMELY important that the order matches the order in the above formula
	nicetitles.m4 = c("Lead Exposed", "Parent-Rated Excellent Health", "Exposed X Parent-Rated Health") 
	parameters.m4 = c(2, 5, 9) 
	mylimits.m4 = c(-50,50)	 
	 

## Create an object of class rqprocess
## Either z1 or z2 may be used for plotting
## The quantreg package in R differs from SAS in that it plots discrete points at values of tau, rather than interpolating a line
## As the number of points approaches infinity, the plot starts to look like a line, but takes longer to run
	z1 = rq(m4, tau = 1:49/50, data=d1)
	z3 = rq(m4, tau = 1:69/70, data=d1)	

## Create a plot with confidence bands [STEP1]
## set the "outside margin area" so we have space to write
## note that this global parameter will remain set until changed, or R is restarted
## only needs to be set once per session
	par(oma=c(3,2,3,.5))	# order: bottom, left, top, right


## Create a plot with confidence bands [STEP2]
## Use the following command for more info:  help(plot.summary.rqs)
	plot(summary(z1, se="boot"), ylim=mylimits.m4,
		mfrow=c(1,3), 					# specify 1 row of plots, 3 columns
		parm = parameters.m4, main=nicetitles.m4,	
		level = .95)	


## Write main title and axis text
	overalltitle = "Interaction on Reading: Lead Exposure X Parent-Rated Health"
	overallxlabel = "x-axes: quantiles of the reading scores distribution"
	overallylabel = "coefficient for predictor term"

	mtext(overallxlabel, side=1, line=1, outer=TRUE, cex=1, font=2)
		# "side=1" means draw on bottom side of plotting area
		# cex refers to the size of the font
	mtext(overallylabel, side=2, line=1, outer=TRUE, cex=1, font=2)	
		# "side=2" means draw on left side of plotting area
	mtext(overalltitle, side=3, line=1, outer=TRUE, cex=1.5, font=2)	
		# "side=3" means draw on top side of plotting area