FILENAME goptions URL "http://192.38.117.59/~linearpredictors/datafiles/goptions.sas";
%include goptions;

FILENAME url URL "http://192.38.117.59/~linearpredictors/datafiles/readVitaminD.sas";
%include url;

/*
Women from Ireland and Poland
*/ 

DATA irlpolwomen;
	SET vitamind (WHERE = (country in (4 6) and category=2));
RUN;

/*
Get the range of BMI for defining x-axis when plotting
*/

PROC MEANS DATA = irlpolwomen MIN MAX;
	VAR bmi;
	OUTPUT OUT = min_max MIN = min_bmi MAX = max_bmi;
RUN;

DATA min_max;
	SET min_max;
	CALL symput('min_bmi',min_bmi);
	CALL symput('max_bmi',max_bmi);
RUN;

/*
Alternative 1
*/

PROC SORT DATA = irlpolwomen;
	BY bmi; 
RUN;

DATA irlpolwomen;
	SET irlpolwomen;
	BY bmi;
	IF first.bmi THEN CALL symput('min_bmi',bmi);
	ELSE IF last.bmi THEN CALL symput('max_bmi',bmi);
RUN;

/*
Alternative 2
*/

PROC SQL NOPRINT;
	SELECT min(bmi), max(bmi) 
	INTO :min_bmi, :max_bmi
	FROM irlpolwomen;
QUIT;

PROC MIXED DATA=irlpolwomen;
	CLASS country;
	MODEL vitd = bmi country / s intercept;
	ODS OUTPUT solutionf=estimates;
RUN;

DATA estimates;
	SET estimates;
	IF effect = 'Intercept' THEN CALL symput('intercept',estimate);
	ELSE IF effect = 'bmi' THEN CALL symput('slope',estimate); 
	IF country=4 THEN CALL symput('intercept_difference',estimate);
RUN;

/*
Creating data with the two endpoints of each regression line in order to plot
*/

DATA plot;
	x=&min_bmi;
	y=&intercept+&slope*&min_bmi;
	line=1;
	output;
	x=&max_bmi;
	y=&intercept+&slope*&max_bmi;
	line=1;
	output;
	x=&min_bmi;
	y=(&intercept+&intercept_difference)+&slope*&min_bmi;
	line=2;
	OUTPUT;
	x=&max_bmi;
	y=(&intercept+&intercept_difference)+&slope*&max_bmi;
	line=2;
	OUTPUT;
RUN;

/*
Defining axes and plotting the two regression lines;
*/

AXIS1 ORDER = (15 to 45 by 5)
	  VALUE = (TICK = 1  '15' 
			   TICK = 2  '20' 
			   TICK = 3  '25' 
			   TICK = 4  '30'
			   TICK = 5  '40') 
	  MINOR = NONE
	  LABEL = ('BMI');

AXIS2 ORDER = (10 to 60 by 10)
	  MINOR = NONE
	  LABEL = (a=90 'Vitamin D');

SYMBOL1 VALUE=NONE I=join LINE=20 CI=black REPEAT=1;
SYMBOL2 VALUE=NONE I=join LINE=1 CI=black REPEAT=1;

LEGEND1 LABEL=('Country:')
		ORDER=(1 2)
		VALUE=('Poland' 'Ireland');

PROC GPLOT DATA = plot;
	PLOT y*x=line / HAXIS = AXIS1
					VAXIS = AXIS2
					LEGEND = LEGEND1;
RUN;
QUIT;