* Output options: all can be turned on/off by adding or removing the NO;
* page number, date, centering, or page breaks, page length and width;
OPTIONS NOnumber NOdate NOcenter FormDlim=' ' PageSize=MAX LineSize=MAX;
* Eliminate SAS default titles and names of tables in output (TRACE ON to show);
TITLE; ODS TRACE OFF; ODS LISTING CLOSE;

***********************************************************************************;
*******   MACRO PROGRAMS TO AUTOMATE CALCULATIONS FOR CHAPTER 6 EXAMPLE     *******;
*******               NOTHING IN HERE NEEDS TO BE CHANGED                   *******;
***********************************************************************************;

* To use FitTest macro;
* FitFewer = Name of ODS InfoCrit table for nested model;
* FitMore  = Name of ODS InfoCrit table for comparison model;
%MACRO FitTest(FitFewer=,FitMore=);
DATA &FitFewer.; LENGTH Name $30.; SET &FitFewer.; Name="&FitFewer."; RUN;
DATA &FitMore.;  LENGTH Name $30.; SET &FitMore.;  Name="&FitMore.";  RUN;
DATA FitCompare; LENGTH Name $30.; SET &FitFewer. &FitMore.; RUN;
DATA FitCompare; SET FitCompare; DevDiff=Lag1(Neg2LogLike)-Neg2LogLike;
     DFdiff=Parms-LAG1(Parms); Pvalue=1-PROBCHI(DevDiff,DFdiff);
     DROP AICC HQIC CAIC; RUN;
TITLE9 "Likelihood Ratio Test for &FitFewer. vs. &FitMore.";
PROC PRINT NOOBS DATA=FitCompare; RUN; TITLE9;
%MEND FitTest;

* To use PseudoR2 macro;
* Ncov =     TOTAL # entries in covariance parameter estimates table;
* CovFewer = Name of ODS CovParms table for nested model;
* CovMore =  Name of ODS CovParms table for comparison model;
%MACRO PseudoR2(NCov=,CovFewer=,CovMore=);
DATA &CovFewer.; LENGTH Name $30.; SET &CovFewer.; Name="&CovFewer."; RUN;
DATA &CovMore.;  LENGTH Name $30.; SET &CovMore.;  Name="&CovMore.";  RUN;
DATA CovCompare; LENGTH Name $30.; SET &CovFewer. &CovMore.; RUN;
DATA CovCompare; SET CovCompare;
     PseudoR2=(LAG&Ncov.(Estimate)-Estimate)/LAG&Ncov.(Estimate); RUN;
DATA CovCompare; SET CovCompare;
     IF CovParm IN("UN(2,1)","UN(3,1)","UN(4,1)","UN(3,2)","UN(4,2)","UN(4,3)")
     THEN DELETE; RUN;
TITLE9 "PsuedoR2 (% Reduction) for &CovFewer. vs. &CovMore.";
PROC PRINT NOOBS DATA=CovCompare; RUN; TITLE9;
%MEND PseudoR2;

***********************************************************************************;
*******             BEGIN DATA MANIPULATION OF CHAPTER 6 EXAMPLE            *******;
*******               CHANGE "filesave" to your directory                   *******;
***********************************************************************************;

* Defining global variable for file location to be replaced in code below;
%LET filesave= C:\Dropbox\PilesOfVariance\Chapter6\SAS;
* Location for SAS files for these models (uses macro variable filesave);
LIBNAME filesave "&filesave.";

* Import chapter 6 stacked data and create time predictor variables;
DATA work.Chapter6; SET filesave.SAS_Chapter6;
* Predictors for polynomial time models;
time1 = session - 1;
time6 = session - 6;
LABEL
time1 = "time1: Session (0=1)"
time6 = "time6: Session (0=6)";
* Predictors for piecewise time models;
     IF session = 1 THEN DO; slope12 = 0; slope26 = 0; END;
ELSE IF session = 2 THEN DO; slope12 = 1; slope26 = 0; END;
ELSE IF session = 3 THEN DO; slope12 = 1; slope26 = 1; END;
ELSE IF session = 4 THEN DO; slope12 = 1; slope26 = 2; END;
ELSE IF session = 5 THEN DO; slope12 = 1; slope26 = 3; END;
ELSE IF session = 6 THEN DO; slope12 = 1; slope26 = 4; END;
LABEL
slope12 = "slope12: Early Practice Slope (Session 1-2)"
slope26 = "slope26: Later Practice Slope (Session 2-6)";
* Piecewise predictors for UNBALANCED DATA (uses per-subject time variables instead);
*      IF session LE 1 THEN DO; * slope12=Time1-Time1; * slope26=Time2-Time2; * END;
* ELSE IF session LE 2 THEN DO; * slope12=Time2-Time1; * slope26=Time2-Time2; * END;
* ELSE IF session LE 3 THEN DO; * slope12=Time2-Time1; * slope26=Time3-Time2; * END;
* ELSE IF session LE 4 THEN DO; * slope12=Time2-Time1; * slope26=Time4-Time2; * END;
* ELSE IF session LE 5 THEN DO; * slope12=Time2-Time1; * slope26=Time5-Time2; * END;
* ELSE IF session LE 6 THEN DO; * slope12=Time2-Time1; * slope26=Time6-Time2; * END;
* LABEL 
* slope12 = "slope12: Early Practice Slope (Session 1-2)"
* slope26 = "slope26: Later Practice Slope (Session 2-6)"; 
RUN;

***********************************************************************************;
*******                        BEGIN CHAPTER 6 MODELS                       *******;
***********************************************************************************;

* Open output directory to save results to;
ODS HTML FILE="&filesave.\SAS_Chapter6_Output.html"
         (URL="SAS_Chapter6_Output.html") STYLE=HTMLBlue;

TITLE1 "Chapter 6 Example: Means by session for RT outcome";
* CLASS= means per session, WAYS= means overall=0 and per session=1;
PROC MEANS MEAN STDERR MIN MAX DATA=work.Chapter6;
     CLASS session;
     WAYS 0 1;
     VAR rt;
RUN; TITLE1;

TITLE1 'Ch 6: 0: Saturated Means, Unstructured Variance Model';
TITLE2 'TOTAL ANSWER KEY';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = session / SOLUTION CL CHISQ DDFM=Satterthwaite;
     REPEATED session / R RCORR TYPE=UN SUBJECT=PersonID;
     LSMEANS session / DIFF=ALL CL;
RUN; TITLE1; TITLE2;

TITLE1 'Ch 6: 1a: Empty Means, E-Only Variance Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt =  / SOLUTION CL CHISQ DDFM=Satterthwaite;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
RUN; TITLE1;

TITLE1 'Ch 6: 1b: Empty Means, Random Intercept Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt =  / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovEmpty;
RUN; TITLE1;

TITLE1 'Ch 6: 2a: Fixed Linear Time, Random Intercept Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovFixLin InfoCrit=FitFixLin;
     ESTIMATE 'Intercept at Session=1 Time=0'    intercept 1 time1 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'    intercept 1 time1 1 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'    intercept 1 time1 2 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'    intercept 1 time1 3 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'    intercept 1 time1 4 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'    intercept 1 time1 5 / CL;
RUN; TITLE1;
     * Call macro to calculate pseudo R2;
     %PseudoR2(NCov=2, CovFewer=CovEmpty, CovMore=CovFixLin);

TITLE1 'Ch 6: 2b: Random Linear Time Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovRandLin InfoCrit=FitRandLin;
     ESTIMATE 'Intercept at Session=1 Time=0'    intercept 1 time1 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'    intercept 1 time1 1 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'    intercept 1 time1 2 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'    intercept 1 time1 3 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'    intercept 1 time1 4 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'    intercept 1 time1 5 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitFixLin, FitMore=FitRandLin);

TITLE1 'Ch 6: 3a: Fixed Quadratic, Random Linear Time Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 time1*time1 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovFixQuad InfoCrit=FitFixQuad;
     ESTIMATE 'Intercept at Session=1 Time=0'         intercept 1 time1 0 time1*time1 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'         intercept 1 time1 1 time1*time1 1 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'         intercept 1 time1 2 time1*time1 4 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'         intercept 1 time1 3 time1*time1 9 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'         intercept 1 time1 4 time1*time1 16 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'         intercept 1 time1 5 time1*time1 25 / CL;
     ESTIMATE 'Linear Slope at Session=1 Time=0'      time1 1 time1*time1 0 / CL;
     ESTIMATE 'Linear Slope at Session=2 Time=1'      time1 1 time1*time1 2 / CL;
     ESTIMATE 'Linear Slope at Session=3 Time=2'      time1 1 time1*time1 4 / CL;
     ESTIMATE 'Linear Slope at Session=4 Time=3'      time1 1 time1*time1 6 / CL;
     ESTIMATE 'Linear Slope at Session=5 Time=4'      time1 1 time1*time1 8 / CL;
     ESTIMATE 'Linear Slope at Session=6 Time=5'      time1 1 time1*time1 10 / CL;
RUN; TITLE1;
     * Call macro to calculate pseudo R2;
     %PseudoR2(NCov=4, CovFewer=CovRandLin, CovMore=CovFixQuad);

TITLE1 'Ch 6: 3b: Random Quadratic Time Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 time1*time1 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 time1*time1 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitRandQuad;
     ESTIMATE 'Intercept at Session=1 Time=0'         intercept 1 time1 0 time1*time1 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'         intercept 1 time1 1 time1*time1 1 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'         intercept 1 time1 2 time1*time1 4 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'         intercept 1 time1 3 time1*time1 9 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'         intercept 1 time1 4 time1*time1 16 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'         intercept 1 time1 5 time1*time1 25 / CL;
     ESTIMATE 'Linear Slope at Session=1 Time=0'      time1 1 time1*time1 0 / CL;
     ESTIMATE 'Linear Slope at Session=2 Time=1'      time1 1 time1*time1 2 / CL;
     ESTIMATE 'Linear Slope at Session=3 Time=2'      time1 1 time1*time1 4 / CL;
     ESTIMATE 'Linear Slope at Session=4 Time=3'      time1 1 time1*time1 6 / CL;
     ESTIMATE 'Linear Slope at Session=5 Time=4'      time1 1 time1*time1 8 / CL;
     ESTIMATE 'Linear Slope at Session=6 Time=5'      time1 1 time1*time1 10 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitFixQuad, FitMore=FitRandQuad);

TITLE1 'Ch 6: 3b: Random Quadratic Time Model (0=Session 6)';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time6 time6*time6 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time6 time6*time6 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ESTIMATE 'Intercept at Session=1 Time=-5'        intercept 1 time6 -5 time6*time6 25 / CL;
     ESTIMATE 'Intercept at Session=2 Time=-4'        intercept 1 time6 -4 time6*time6 16 / CL;
     ESTIMATE 'Intercept at Session=3 Time=-3'        intercept 1 time6 -3 time6*time6 9 / CL;
     ESTIMATE 'Intercept at Session=4 Time=-2'        intercept 1 time6 -2 time6*time6 4 / CL;
     ESTIMATE 'Intercept at Session=5 Time=-1'        intercept 1 time6 -1 time6*time6 1 / CL;
     ESTIMATE 'Intercept at Session=6 Time= 0'        intercept 1 time6 0  time6*time6 0 / CL;
     ESTIMATE 'Linear Slope at Session=1 Time=-5'     time6 1 time6*time6 -10 / CL;
     ESTIMATE 'Linear Slope at Session=2 Time=-4'     time6 1 time6*time6 -8 / CL;
     ESTIMATE 'Linear Slope at Session=3 Time=-3'     time6 1 time6*time6 -6 / CL;
     ESTIMATE 'Linear Slope at Session=4 Time=-2'     time6 1 time6*time6 -4 / CL;
     ESTIMATE 'Linear Slope at Session=5 Time=-1'     time6 1 time6*time6 -2 / CL;
     ESTIMATE 'Linear Slope at Session=6 Time= 0'     time6 1 time6*time6 0 / CL;
RUN; TITLE1;

TITLE1 'Ch 6: 4a: Fixed Slope12, Fixed Slope26, Random Intercept Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = slope12 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovFix12Fix26 InfoCrit=FitFix12Fix26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 slope12 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 slope12 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 slope12 1 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 slope12 1 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 slope12 1 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 slope12 1 slope26 4 / CL;
     ESTIMATE 'Difference between slope12 and slope26'     slope12 -1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate pseudo R2;
     %PseudoR2(NCov=2, CovFewer=CovEmpty, CovMore=CovFix12Fix26);

TITLE1 'Ch 6: 4b: Random Slope12, Fixed Slope26 Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = slope12 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT slope12 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitRand12Fix26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 slope12 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 slope12 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 slope12 1 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 slope12 1 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 slope12 1 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 slope12 1 slope26 4 / CL;
     ESTIMATE 'Difference between slope12 and slope26'     slope12 -1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitFix12Fix26, FitMore=FitRand12Fix26);

TITLE1 'Ch 6: 4c: Random Slope12, Random Slope26 Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = slope12 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT slope12 slope26 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitRand12Rand26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 slope12 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 slope12 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 slope12 1 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 slope12 1 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 slope12 1 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 slope12 1 slope26 4 / CL;
     ESTIMATE 'Difference between slope12 and slope26'     slope12 -1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitRand12Fix26, FitMore=FitRand12Rand26);

TITLE1 'Ch 6: Random Slope12, Random Slope26 Model + Fixed Quadratic Slope26';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = slope12 slope26 slope26*slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT slope12 slope26 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
RUN; TITLE1;

TITLE1 'Ch 6: 5a: Fixed Time, Fixed Slope26, Random Intercept Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT CovParms=CovFix16Fix26 InfoCrit=FitFix16Fix26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 time1 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 time1 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 time1 2 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 time1 3 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 time1 4 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 time1 5 slope26 4 / CL;
     ESTIMATE 'Rate of change from session 2 to 6'         time1 1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate pseudo R2;
     %PseudoR2(NCov=2, CovFewer=CovEmpty, CovMore=CovFix16Fix26);

TITLE1 'Ch 6: 5b: Random Time, Fixed Slope26 Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitRand16Fix26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 time1 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 time1 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 time1 2 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 time1 3 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 time1 4 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 time1 5 slope26 4 / CL;
     ESTIMATE 'Rate of change from session 2 to 6'         time1 1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitFix16Fix26, FitMore=FitRand16Fix26);

TITLE1 'Ch 6: 5c: Random Time, Random Slope26 Model';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=REML;
     CLASS PersonID session;
     MODEL rt = time1 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 slope26 / G GCORR V VCORR TYPE=UN SUBJECT=PersonID;
     REPEATED session / R RCORR TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitRand16Rand26;
     ESTIMATE 'Intercept at Session=1 Time=0'              intercept 1 time1 0 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=2 Time=1'              intercept 1 time1 1 slope26 0 / CL;
     ESTIMATE 'Intercept at Session=3 Time=2'              intercept 1 time1 2 slope26 1 / CL;
     ESTIMATE 'Intercept at Session=4 Time=3'              intercept 1 time1 3 slope26 2 / CL;
     ESTIMATE 'Intercept at Session=5 Time=4'              intercept 1 time1 4 slope26 3 / CL;
     ESTIMATE 'Intercept at Session=6 Time=5'              intercept 1 time1 5 slope26 4 / CL;
     ESTIMATE 'Rate of change from session 2 to 6'         time1 1 slope26 1 / CL;
RUN; TITLE1;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitRand16Fix26, FitMore=FitRand16Rand26);

TITLE1 'Ch 6: 6a Fixed Asymptote Fixed Amount Fixed Rate Model';
PROC NLMIXED DATA=work.Chapter6 METHOD=GAUSS TECH=NEWRAP GCONV=1e-12;
* Must define all parameters to be estimated and provide start values;
* First line is fixed effects, second line is variances (in SD metric);
  PARMS fasymp= 1600 famount=300 frate=-1
        SDE=600;
* Setting up level-2 equations;
  b0i = fasymp;
  b1i = famount;
  b2i = frate;
* Setting up level-1 equation WITHOUT level-1 residual;
  PredY = (b0i) + (b1i*EXP(b2i*time1));
* Telling it which DV, defining level-1 residual;
* RTs is normally distributed with a mean of 'PredY' and a variance of 'VarE';
  MODEL rt ~ normal(PredY, sdE*sdE);
* Labeling estimated parameters;
  ESTIMATE 'Fixed Asymptote'                    fasymp;
  ESTIMATE 'Fixed Amount'                       famount;
  ESTIMATE 'Fixed Rate'                         frate;
  ESTIMATE 'Level-1 Residual E Variance'        sdE*sdE;
* Creating extra parameters and predicted means;
  ESTIMATE 'Fixed Intercept'                    fasymp+famount;
  ESTIMATE 'Session=1 Time=0 Predicted Mean'    fasymp+(famount*EXP(frate*0));
  ESTIMATE 'Session=2 Time=1 Predicted Mean'    fasymp+(famount*EXP(frate*1));
  ESTIMATE 'Session=3 Time=2 Predicted Mean'    fasymp+(famount*EXP(frate*2));
  ESTIMATE 'Session=4 Time=3 Predicted Mean'    fasymp+(famount*EXP(frate*3));
  ESTIMATE 'Session=5 Time=4 Predicted Mean'    fasymp+(famount*EXP(frate*4));
  ESTIMATE 'Session=6 Time=5 Predicted Mean'    fasymp+(famount*EXP(frate*5));
RUN; TITLE1;

TITLE1 'Ch 6: 6b Random Asymptote Fixed Amount Fixed Rate Model';
PROC NLMIXED DATA=work.Chapter6 METHOD=GAUSS TECH=NEWRAP GCONV=1e-12;
* Must define all parameters to be estimated and provide start values;
* First line is fixed effects, second line is variances (in SD metric);
  PARMS fasymp= 1675 famount=284 frate=-.7
        sdE=474 sdU0=10;
* Setting up level-2 equations;
  b0i = fasymp + U0i;
  b1i = famount;
  b2i = frate;
* Setting up level-1 equation WITHOUT level-1 residual;
  PredY = (b0i) + (b1i*EXP(b2i*time1));
* Telling it which DV, defining level-1 residual;
* RTs is normally distributed with a mean of 'PredY' and a variance of 'VarE';
  MODEL rt ~ normal(PredY, sdE*sdE);
* Defining random effects: normally distributed with means and variances;
  RANDOM U0i ~ normal([0],[sdU0*sdU0]) SUBJECT=PersonID;
* Labeling estimated parameters;
  ESTIMATE 'Fixed Asymptote'                        fasymp;
  ESTIMATE 'Fixed Intercept'                        fasymp+famount;
  ESTIMATE 'Fixed Amount'                           famount;
  ESTIMATE 'Fixed Rate'                             frate;
  ESTIMATE 'Level-1 Residual E Variance'            sdE*sdE;
  ESTIMATE 'Level-2 Random Asymptote U0 Variance'   sdU0*sdU0;
* Creating extra parameters and predicted means;
  ESTIMATE 'Fixed Intercept'                        fasymp+famount;
  ESTIMATE 'Session=1 Time=0 Predicted Mean'        fasymp+(famount*EXP(frate*0));
  ESTIMATE 'Session=2 Time=1 Predicted Mean'        fasymp+(famount*EXP(frate*1));
  ESTIMATE 'Session=3 Time=2 Predicted Mean'        fasymp+(famount*EXP(frate*2));
  ESTIMATE 'Session=4 Time=3 Predicted Mean'        fasymp+(famount*EXP(frate*3));
  ESTIMATE 'Session=5 Time=4 Predicted Mean'        fasymp+(famount*EXP(frate*4));
  ESTIMATE 'Session=6 Time=5 Predicted Mean'        fasymp+(famount*EXP(frate*5));
RUN; TITLE1;

TITLE1 'Ch 6: 6c Random Asymptote Random Amount Fixed Rate Model';
PROC NLMIXED DATA=work.Chapter6 METHOD=GAUSS TECH=NEWRAP GCONV=1e-12;
* Must define all parameters to be estimated and provide start values;
* First line is fixed effects, second line is variances (in SD metric);
  PARMS fasymp=1675 famount=284 frate=-.7
        sdE=184 sdU0=447 sdU01=1 sdU1=10;
* Setting up level-2 equations;
  b0i = fasymp  + U0i;
  b1i = famount + U1i;
  b2i = frate;
* Setting up level-1 equation WITHOUT level-1 residual;
  PredY = (b0i) + (b1i*EXP(b2i*time1));
* Telling it which DV, defining level-1 residual;
* RTs is normally distributed with a mean of 'y' and a variance of 'VarE';
  MODEL rt ~ normal(PredY, sdE*sdE);
* Defining random effects: normally distributed with means and variances;
  RANDOM U0i U1i ~ normal([0,0],[sdU0*sdU0,sdU01*sdU01,sdU1*sdU1]) SUBJECT=PersonID;
* Labeling estimated parameters;
  ESTIMATE 'Fixed Asymptote'                            fasymp;
  ESTIMATE 'Fixed Amount'                               famount;
  ESTIMATE 'Fixed Rate'                                 frate;
  ESTIMATE 'Level-1 Residual E Variance'                sdE*sdE;
  ESTIMATE 'Level-2 Random Asymptote U0 Variance'       sdU0*sdU0;
  ESTIMATE 'Level-2 Asymptote-Amount U01 Covariance'    sdU01*sdU01;
  ESTIMATE 'Level-2 Random Amount U1 Variance'          sdU1*sdU1;
  ESTIMATE 'Level-2 Asymptote-Amount Correlation'       (sdU01*sdU01)/(sdU0*sdU1);
* Creating extra parameters and predicted means;
  ESTIMATE 'Fixed Intercept'                            fasymp+famount;
  ESTIMATE 'Session=1 Time=0 Predicted Mean'            fasymp+(famount*EXP(frate*0));
  ESTIMATE 'Session=2 Time=1 Predicted Mean'            fasymp+(famount*EXP(frate*1));
  ESTIMATE 'Session=3 Time=2 Predicted Mean'            fasymp+(famount*EXP(frate*2));
  ESTIMATE 'Session=4 Time=3 Predicted Mean'            fasymp+(famount*EXP(frate*3));
  ESTIMATE 'Session=5 Time=4 Predicted Mean'            fasymp+(famount*EXP(frate*4));
  ESTIMATE 'Session=6 Time=5 Predicted Mean'            fasymp+(famount*EXP(frate*5));
RUN;

* No convergence;
TITLE1 'Ch 6: 6d Random Asymptote Random Amount Random Rate Model';
PROC NLMIXED DATA=work.Chapter6 METHOD=GAUSS TECH=NEWRAP GCONV=1e-12;
* Must define all parameters to be estimated and provide start values;
* First line is fixed effects, second line is variances (in SD metric);
  PARMS fasymp=1678.25 famount=-282.72 frate=-.7323
        sdE=130.88 sdU0=426.79 sdU01=87.68 sdU1=290.57
        sdU02=5.1343 sdU12=.08921 sdU2=1.3483;
* Setting up level-2 equations;
  b0i = fasymp  + U0i;
  b1i = famount + U1i;
  b2i = frate   + U2i;
* Setting up level-1 equation WITHOUT level-1 residual;
  PredY = (b0i) + (b1i*EXP(b2i*time1));
* Telling it which DV, defining level-1 residual;
* RTs is normally distributed with a mean of 'PredY' and a variance of 'VarE';
MODEL rt ~ normal(PredY , sdE*sdE);
* Defining random effects: normally distributed with means and variances;
RANDOM  U0i U1i U2i ~ normal([0,0,0],[sdU0*sdU0,sdU01*sdU01,sdU1*sdU1,
                                     sdU02*sdU02,sdU12*sdU12,sdU2*sdU2]) SUBJECT=PersonID;
* Labeling estimated parameters;
  ESTIMATE 'Fixed Asympote'                             fasymp;
  ESTIMATE 'Fixed Amount'                               famount;
  ESTIMATE 'Fixed Rate'                                 frate;
  ESTIMATE 'Level-1 Residual E Variance'                sdE*sdE;
  ESTIMATE 'Level-2 Random Asymptote U0 Variance'       sdU0*sdU0;
  ESTIMATE 'Level-2 Asymptote-Amount U01 Covariance'    sdU01*sdU01;
  ESTIMATE 'Level-2 Random Amount U1 Variance'          sdU1*sdU1;
  ESTIMATE 'Level-2 Asymptote-Rate U02 Variance'        sdU02*sdU02;
  ESTIMATE 'Level-2 Amount-Rate U12 Covariance'         sdU12*sdU12;
  ESTIMATE 'Level-2 Random Rate U2 Variance'            sdU2*sdU2;
* Creating extra parameters and predicted means;
  ESTIMATE 'Fixed Intercept'                            fasymp+famount;
  ESTIMATE 'Session=1 Time=0 Predicted Mean'            fasymp+(famount*EXP(frate*0));
  ESTIMATE 'Session=2 Time=1 Predicted Mean'            fasymp+(famount*EXP(frate*1));
  ESTIMATE 'Session=3 Time=2 Predicted Mean'            fasymp+(famount*EXP(frate*2));
  ESTIMATE 'Session=4 Time=3 Predicted Mean'            fasymp+(famount*EXP(frate*3));
  ESTIMATE 'Session=5 Time=4 Predicted Mean'            fasymp+(famount*EXP(frate*4));
  ESTIMATE 'Session=6 Time=5 Predicted Mean'            fasymp+(famount*EXP(frate*5));
RUN; TITLE1;

TITLE1 'Ch 6: 0: Saturated Means, Unstructured Variance Model';
TITLE2 'Using ML Instead of REML';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=ML;
     CLASS PersonID session;
     MODEL rt = session / SOLUTION CL CHISQ DDFM=Satterthwaite;
     REPEATED session / TYPE=UN SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitMLSatUN;
RUN; TITLE1; TITLE2;

TITLE1 'Ch 6: 1b: Empty Means, Random Intercept Model';
TITLE2 'Using ML Instead of REML';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=ML;
     CLASS PersonID session;
     MODEL rt =  / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT / TYPE=UN SUBJECT=PersonID;
     REPEATED session / TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitMLEmptyRI;
RUN; TITLE1; TITLE2;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitMLEmptyRI, FitMore=FitMLSatUN);

TITLE1 'Ch 6: 3b: Random Quadratic Time Model';
TITLE2 'Using ML instead of REML';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=ML;
     CLASS PersonID session;
     MODEL rt = time1 time1*time1 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT time1 time1*time1 / TYPE=UN SUBJECT=PersonID;
     REPEATED session / TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitMLRandQuad;
RUN; TITLE1; TITLE2;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitMLEmptyRI, FitMore=FitMLRandQuad);
     %FitTest(FitFewer=FitMLRandQuad, FitMore=FitMLSatUN);

TITLE1 'Ch 6: 4c: Random Slope12, Random Slope26 Model';
TITLE2 'Using ML instead of REML';
PROC MIXED DATA=work.Chapter6 COVTEST NOCLPRINT NAMELEN=100 IC METHOD=ML;
     CLASS PersonID session;
     MODEL rt = slope12 slope26 / SOLUTION CL CHISQ DDFM=Satterthwaite;
     RANDOM INTERCEPT slope12 slope26 / TYPE=UN SUBJECT=PersonID;
     REPEATED session / TYPE=VC SUBJECT=PersonID;
     ODS OUTPUT InfoCrit=FitMLRand12Rand26;
RUN; TITLE1; TITLE2;
     * Call macro to calculate LRT for nested models;
     %FitTest(FitFewer=FitMLEmptyRI, FitMore=FitMLRand12Rand26);
     %FitTest(FitFewer=FitMLRand12Rand26, FitMore=FitMLSatUN);

****** END CHAPTER 6 MODELS ******;

* Close output directory;
ODS HTML CLOSE;