/**************************************
 ANCOVA Using Proc Reg and Proc GLM
 FILENAME: ancova_example.sas
****************************************/
options pageno=1 nodate;
title;

/*Import data from SPSS portable file*/
filename file1 "E:\510\2006\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=red interpol=rl value=dot line=1;
symbol2 color=blue 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 using Proc Reg*/
title "ANCOVA for Males and Females";
title2 "Relationship of Height to Age";
proc reg data=htwt2;
  where age <=19;
  model height = age female fem_age;
  plot rstudent. * predicted.;
  output out=outreg1 p=predict1 r=resid1 rstudent=rstud1;
run;
proc univariate data=outreg1 normal;
  var resid1;
  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 = centage female 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;

/**************************************
 Another ANCOVA Model for the Cars
 Data Set
****************************************/
/*Import data from SPSS portable file*/
filename file2 "C:\Documents and Settings\Kwelch\Desktop\b510\cars.por";
proc convert spss=file2 out=cars;
run;
title "Contents of Cars Data Set";
proc contents data=cars;
run;

title "Descriptive Statistics for Cars Data Set";
proc means data=cars;
run;


title "Oneway Frequency Tabulation of Cylinder for Cars Data Set";
proc freq data=cars;
  tables cylinder;
run;

goptions reset = all;
goptions device=win target=winprtm;
symbol1 color=black interpol=rl value=dot;
symbol2 color=black interpol=rl value=star;
symbol3 color=black interpol=rl value=square;
title "Scatterplot of Horsepower by Weight";
title2 "For 4, 6, and 8 Cylinder Vehicles";
proc gplot data=cars;
  where cylinder in(4,6,8);
  plot horse * weight  = cylinder;
run; quit;

/*Create New data Set with Dummy Variables*/

data cars2;
  set cars;
  
  if cylinder not=. then do;
    fourcyl = 0;
	if cylinder = 4 then fourcyl=1;
    sixcyl  = 0;
    if cylinder = 6 then sixcyl=1;
    eightcyl = 0;
    if cylinder = 8 then eightcyl=1;
  end;
  
   fourcyl_wt  = fourcyl*weight;
   sixcyl_wt   = sixcyl *weight;
   eightcyl_wt = eightcyl*weight;

   weight1000 = weight/1000;
run; 

title "ANCOVA Model for Cars Data";
title2 "Using Proc Reg";
proc reg data=cars2;
  where cylinder in (4,6,8);
  model horse = sixcyl eightcyl WEIGHT
                sixcyl_wt eightcyl_wt;
run;quit;
 
title "ANCOVA Model for Cars Data";
title2 "Using Proc GLM";
proc glm data=cars2;
  where cylinder in (4,6,8);
  class cylinder;
  model horse = cylinder weight cylinder*weight;
  estimate "Slope for 4 Cylinders" weight 1 cylinder*weight 1 0 0;
  estimate "Slope for 6 Cylinders" weight 1 cylinder*weight 0 1 0;
  estimate "Slope for 8 Cylinders" weight 1 cylinder*weight 0 0 1;
run; quit;

proc sort data=cars2;
  by cylinder;
run;

title "Separate Regression for Each Number of Cylinders";
proc reg data=cars2;
 where cylinder in (4,6,8);
 by cylinder;
 model horse = weight;
run; quit;

title "Separate Regression for Each Number of Cylinders";
title2 "Using Weight/1000";
proc reg data=cars2;
 where cylinder in (4,6,8);
 by cylinder;
 model horse = weight1000;
run; quit;