libname in 'e:\singer';
options nocenter nodate nonumber nolabel;

/*Descriptive statistics of level-1 variables*/
proc means data = in.hsb12 mean std;
  var ses mathach;
run;

/*Descriptive statistics of level-2 variables*/
proc sort data = in.hsb12;
  by school;
run;
data level2;
  set in.hsb12;
  by school;
  if last.school;
  keep school meanses sector;
run;
proc means data = level2 mean std;
  var meanses sector;
run;

/***********Model 1*******************************/

title "Unconditional Means Model";
proc mixed data = in.hsb12 covtest noclprint;
   class school;
   model mathach = / solution;
   random intercept / subject = school;
run;

/*************Model 2****************************/

title "Model 2: Including Effects of School Level (level 2) Predictors";
title2 "-- predicting mathach from meanses-- page 331";

proc mixed data = in.hsb12 covtest noclprint;
   class school;
   model mathach = meanses / solution ddfm = bw;
   random intercept / subject = school;
run;
/*comparing the degrees of freedom with and w/o option ddfm = bw*/ 
proc mixed data = in.hsb12 covtest noclprint;
   class school;
   model mathach = meanses / solution ;
   random intercept / subject = school;
run;

/***********Model 3****************************/

title "Model 3: Including Effects of Student-Level Predictors";
title2 "--predicting mathach from cses -- page 334.";

data hsbc;
  set in.hsb12;
    cses = ses - meanses;
run;

proc mixed data = hsbc noclprint covtest noitprint;
  class school;
  model mathach = cses / solution ddfm = bw notest;
  random intercept cses / subject = school type = un gcorr;
run;

/*************Model 4**************************/

title "Model 4: Including Both Level-1 and Level-2 Predictors";
title2 "--predicting mathach from meanses, sector, cses and ";
title3 "cross level interaction of  meanses and sector with cses-- page 337.";

proc mixed data = hsbc noclprint covtest noitprint;
  class school;
  model mathach = meanses sector cses meanses*cses sector*cses 
                  / solution ddfm = bw notest;
  random intercept cses / subject = school type = un;
run;

/**********Comment on Model 4, creating plot**************/

proc univariate data = hsbc;
  var meanses;
run;
data toplot;
  set hsbc;
  if meanses <= -0.317 then do;
		ms = -0.317;
 		strata = "Low";   end;
  else if meanses >= 0.333 then do;
		ms = 0.333;
		strata = "Hig";   end;
  else do; ms = 0.038; strata = "Med" ; end;
  predicted = 12.1136 + 5.3391*ms * 1.2167*sector + 2.9388*cses +
              1.0389*ms*cses - 1.6426*sector*cses;
run;
proc sort data = toplot;
   by strata;
run;
goptions reset = all;
symbol1 v = none i = join c = red ;
symbol2 v = none i = join c = blue  ;
axis1 order = (-4 to 3 by 1) minor = none label=("Group Centered SES");
axis2 order = (0 to 22 by 2) minor = none label=(a = 90 "Math Achievement Score");
proc gplot data = toplot;
   by strata;
   plot predicted*cses = sector / vaxis = axis2 haxis = axis1; 
run;
quit; 

/**********Comment on Model , possible interactions*******************************/

proc mixed data = hsbc noclprint covtest noitprint;
  class school;
  model mathach = meanses sector cses meanses*sector 
                  meanses*cses sector*cses meanses*sector*cses 
                  / solution ddfm = bw notest;
  random intercept cses / subject = school type = un;
run;
 
/********Comment on Model , using simpler model********************/

proc mixed data = hsbc noclprint covtest noitprint;
  class school;
  model mathach = meanses sector cses meanses*cses sector*cses / solution ddfm = bw notest;
  random intercept / subject = school;
run;

/*******************Comment on Comparing models*******************/

data pvalue;
  df = 2; chisq = 46504.8 - 46503.7;
  pvalue = 1 - probchi(chisq, df);
run;
proc print data = pvalue noobs;
run;

/*******Linear Growth Model**************************/

title "Model 1: Unconditional Linear Growth Model -- page 340";

proc mixed data = in.willett noclprint covtest;
  class id;
  model y = time /solution ddfm = bw notest;
  random intercept time / subject = id type = un;
run;
/*******Comment on Model 1, creating plot****************/
goptions reset = all;
symbol i = join r = 20;
proc gplot data = in.willett; 
  plot y*time = id;
  where id <=20;
run;
quit;
/************Model 2**************************/          

title "Model 2: A Linear Growth Model with a Person-Level Covariance";
title2 "-- predicting y by time and centered covar -- page 344";

data willett;
  set in.willett;
  wave = time;
  ccovar = covar -  113.4571429;
run;
proc mixed data = willett noclprint covtest;
  class id;
  model y = time ccovar time*ccovar /solution ddfm = bw notest;
  random intercept time / subject = id type = un gcorr;
run;

/************Model 3: covarinace structure within persons***************/

title "A. Compound Symmetry";
proc mixed data = willett covtest noitprint;
  class id wave;
  model y = time / s notest;
  repeated wave /type = cs subject = id r;
run;

title "B.Unstructured";
proc mixed data = willett covtest noitprint;
  class id wave;
  model y = time / s notest;
  repeated wave /type = un subject = id r;
run;

title "C. AR(1)";
proc mixed data = willett covtest noitprint;
  class id wave;
  model y = time / s notest;
  repeated wave /type = ar(1) subject = id r;
run;

title "D: AR(1) with random effect covariates";
proc mixed data = willett noclprint covtest;
  class id wave;
  model y = time ccovar time*ccovar /solution ddfm = bw notest;
  random intercept time / subject = id type = un gcorr;
  repeated wave /type = ar(1) subject = id r;
run;