/**************************************
 ANCOVA Using Proc Reg and Proc GLM
 FILENAME: ancova_example.sas
****************************************/
OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";
options pageno=1 nodate;
title;

/*Import data from SPSS portable file*/
filename file1 "d:\510\2007\data\htwt.por";
proc convert spss=file1 out=htwt;
run;
title "Contents of HTWT Data Set";
proc contents data=htwt;
run;
title "Descriptive Statistics for HTWT Data Set";
proc means data=htwt;
run;
title "Oneway Frequency Tabulation for Sex for HTWT Data Set";
proc freq data=htwt;
tables sex;
run;


/*Create a new data set with new variables*/
data htwt2;
  set htwt;

  /*Create dummy variables for female*/
  if sex="f" then female=1;
  if sex="m" then female=0;
  
  /*Center age at 16 years*/
  centage = age - 16.5;

  /*Create interactions*/
  fem_age = female * age;
  fem_centage = female * centage;

run;

goptions reset = all;
goptions device=win target=winprtc;
symbol1 color=black interpol=rl value=dot line=1;
symbol2 color=black interpol=rl value=star line=2;
title "Scatterplot of Height by Age";
title2 "For Males and Females";
proc gplot data=htwt2;
  where age <= 19;
  plot height * age = sex;
run; quit;

/*ANCOVA model using Proc Reg with Dummy Variables and Interactions*/title "ANCOVA using Proc Reg for Males and Females";
title2 "Relationship of Height to Age";
proc reg data=htwt2;
  where age <=19;
  model height =female  age fem_age;
  plot rstudent. * predicted.;
  output out=outreg1 p=predict1 r=resid1 rstudent=rstud1;
run; quit;
proc univariate data=outreg1 normal;
  var rstud1;
  histogram / normal;
  probplot / normal (mu=est sigma=est);
run;

goptions reset = all;
goptions device=win target=winprtm;
symbol1 color=black interpol=rl value=dot;
symbol2 color=black interpol=rl value=star;
axis1 order = 0 to 20 by 2;
axis2 order = 0 to 80 by 10;
title "Scatterplot Showing Origin";
proc gplot data=htwt2;
  where age <= 19;
  plot height * age = sex / haxis = axis1
                            vaxis = axis2;
run; quit;


title "ANCOVA for Males and Females";
title2 "Relationship of Height to Centered Age";
proc reg data=htwt2;
  where age <=19;
  model height = female centage  fem_centage;
  plot rstudent. * predicted.;
  output out=outreg2 p=predict2 r=resid2 rstudent=rstud2;
run;

/*GLM model with interactions*/
title "ANCOVA model using GLM";
title "Relationship of Height to Centered AGE";
proc glm data=htwt2;
  where age <=19;
  class sex;
  model height = sex centage sex*centage / solution;
  estimate "intercept males" intercept 1 sex 0 1;
  estimate "slope males" centage 1 sex*centage 0 1;
  estimate "intercept females" intercept 1 sex 1 0;
  estimate "slope females" centage 1 sex*centage 1 0;
run; quit;

/*Separate Regression Models for Females and Males*/
proc sort data = htwt2;
  by sex;
run;
title "Separate Regressions for Females and Males";
proc reg data = htwt2;
  by sex;
  model height = centage;
run; quit;