The Code

Run it in the browser

The notebook implements the ANCOVA at both timepoints, the interaction models, the adjusted means, and the checks the site adds: retention by arm, assumption tests on the reported model's residuals, and cell sizes behind every subgroup estimate.

Open in Colab View on GitHub

The notebook runs on synthetic data. The trial data is restricted patient-level clinical information and cannot be redistributed, so the notebook generates a dataset from this study's published aggregates: the same four arms, the same baseline-to-outcome relationship, and comparable retention.

Numbers it produces are therefore simulated and are not the study's results. The real figures are on the analysis page and in the executed SAS output. A flag at the top switches the notebook onto the real data for anyone with an authorised copy.

The pipeline

Import and subset

Read the trial dataset and take one timepoint at a time. Month 3 and month 6 are analysed as separate cross-sections.

Describe

PROC MEANS by arm, then homogeneity of variance with Levene's test.

Model

PROC GLM with the outcome on arm plus the baseline covariate. This is the model the results come from.

Interactions

Two-way and three-way models adding gender and race, with LSMEANS for the cell estimates.

Longitudinal

PROC MIXED across months 3 and 6 with an autoregressive covariance structure and patient as the subject effect.

Two changes were made before publishing the program. The LIBNAME pointed at a Windows path on a university virtual desktop and contained an account identifier; it is now a &datapath. macro variable set at the top. A header comment and a title statement carrying an internal project label were removed, the title because it would otherwise print on every page of the output.

Two proc print statements that would display patient rows are present in the original and are commented out. They must stay that way: their output would not be publishable.

The SAS program

psychosocial_treatments_analysis.sas, 621 lines.

Source on GitHub Executed output

Show the full program (621 lines)

psychosocial_treatments_analysis.sasSAS

/* Analysis of Psychosocial Treatments for Cocaine Dependence
   Data: NIDA Collaborative Cocaine Treatment Study.

   Set this to the directory holding the trial dataset. The data is restricted and
   is not distributed with this repository. */
%let datapath = C:\path\to\trial\data;


ODS HTML CLOSE;
ODS HTML;
dm 'clear log';
dm 'clear output';
options ls=90 ps=50 formdlim=' ';

LIBNAME mylib "&datapath.";

/* Prepare data for cross-sectional analysis at month 3 and month 6*/
/*
proc print data = mylib.finalcoke(obs=30);
run;
*/

proc format;
	value tx_condfmt 1='IDC+GDC' 2='CT+GDC' 3='SE+GDC' 4='GDC';
	value genderfmt 0='Male' 1='Female';
	value racefmt 0='Caucasian' 1='Non-Caucasian';
	value crackfmt 1='CRACK' 0='Snorts/IV';
run;

data cokeM3; 
	set mylib.finalcoke;
	if MONTH=3;
run;

/*
proc print data = cokeM3(obs=30);
    var BEH_ARS ps_BEH_ARS;
run;
*/

data cokeM6; 
	set mylib.finalcoke; 
	if MONTH=6;
run;


/* Descriptive Statistics: Look at plots and averages to get an overview for outcomes. */
proc means data=cokeM3 n mean std min max;
    class TX_COND;
    var BEH_ARS;
	format TX_COND tx_condfmt.;
    title 'Month 3: Descriptive Statistics for BEH_ARS vs TX_COND';
run;


proc sgplot data=cokeM3;
    vbox BEH_ARS / category=TX_COND;
	format TX_COND tx_condfmt.;
    title 'Month 3: Boxplot of BEH_ARS vs TX_COND';
run;

proc means data=cokeM6 n mean std min max;
    class TX_COND;
    var BEH_ARS;
	format TX_COND tx_condfmt.;
    title 'Month 6: Descriptive Statistics for BEH_ARS vs TX_COND';
run;


proc sgplot data=cokeM6;
    vbox BEH_ARS / category=TX_COND;
	format TX_COND tx_condfmt.;
    title 'Month 6: Boxplot of BEH_ARS vs TX_COND';
run;

proc sgpanel data=cokeM3;
    panelby TX_COND / layout=columnlattice;
    vline RACE / response=BEH_ARS group=GENDER stat=mean markers;
    keylegend / title="Gender";
    label AAE = "Mean BEH_ARS Score";
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month 3: Plot of BEH_ARS vs RACE vs GENDER';
run;

proc sgpanel data=cokeM6;
    panelby TX_COND / layout=columnlattice;
    vline RACE / response=BEH_ARS group=GENDER stat=mean markers;
    keylegend / title="Gender";
    label AAE = "Mean BEH_ARS Score";
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month 6: Plot of BEH_ARS vs RACE vs GENDER';
run;


/* Check for Homogeneity of Variance and Normality of residuals */

proc glm data=cokeM3;
        class TX_COND;
        model BEH_ARS = TX_COND;
        means TX_COND / HOVTEST=LEVENE; 
        title 'Month 3: Homogeneity of Variance Check';
run;
quit;

proc glm data=cokeM3 NOPRINT;
	class TX_COND;
	model BEH_ARS = TX_COND;
	output OUT=resid_data R=residuals;
RUN;
QUIT;

proc univariate data=resid_data NORMAL;
        var residuals;
        Histogram residuals / NORMAL;
        QQPLOT residuals / NORMAL(MU=EST SIGMA=EST);
        title 'Month3: Normality Check of Residuals';
run;

proc glm data=cokeM6;
        class TX_COND;
        model BEH_ARS = TX_COND;
        means TX_COND / HOVTEST=LEVENE; 
        title 'Month 6: Homogeneity of Variance Check';
run;
quit;

proc glm data=cokeM6 NOPRINT;
	class TX_COND;
	model BEH_ARS = TX_COND;
	output OUT=resid_data R=residuals;
run;
quit;

proc univariate data=resid_data NORMAL;
        var residuals;
        Histogram residuals / NORMAL;
        QQPLOT residuals / NORMAL(MU=EST SIGMA=EST);
        title 'Month6: Normality Check of Residuals';
run;

/* Month 6 normality is volated; lets run boxcox to find optimal lambda */
ODS GRAPHICS ON;
proc transreg data=cokeM6 maxiter=100 plots=(boxcox) DETAILS;
    model BOXCOX(BEH_ARS / lambda=-3 to 3 by 0.1 CONVENIENT) = CLASS(TX_COND);
    title 'Month 6: Box-Cox Transformation Analysis';
run;
/* Lambda = 1; no transformation needed */

/******************/
/* Models Month 3 */
/******************/
proc glm data=cokeM3;
	class TX_COND;
	model BEH_ARS = TX_COND ps_BEH_ARS / SOLUTION;
	lsmeans TX_COND / PDIFF ADJUST=BON; 
	format TX_COND tx_condfmt.;
	title 'Month 3: Model 1: Primary ANCOVA BEH_ARS = TX_COND ps_BEH_ARS';
run;
quit;


proc glm data=cokeM3;    
	class TX_COND GENDER;
	model BEH_ARS = TX_COND|GENDER ps_BEH_ARS / SOLUTION;
	/* Slice interaction to see TX_CODE effect within Genders */
	LSMEANS TX_COND*GENDER / SLICE=GENDER PDIFF ADJUST=TUKEY;        
	/* Estimate specific contrast*/
	ESTIMATE 'IDC+GDC versus GDC' TX_COND 1 0 0 -1;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month3: Model 2: 2-way Interaction BEH_ARS = TX_COND*GENDER ps_BEH_ARS';
run;
quit;

proc glm data=cokeM3;
	class TX_COND GENDER RACE;
	model BEH_ARS = TX_COND|GENDER|RACE ps_BEH_ARS / SOLUTION;
	LSMEANS TX_COND*GENDER*RACE / PDIFF=ALL;
	lsmeans TX_COND*GENDER*RACE / slice=GENDER*RACE;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month3: Model 3: 3-way Interaction BEH_ARS = TX_COND*GENDER*RACE ps_BEH_ARS';
run;
quit;

/******************/
/* Models Month 6 */
/******************/
proc glm data=cokeM6;
	class TX_COND;
	model BEH_ARS = TX_COND ps_BEH_ARS / SOLUTION;
	lsmeans TX_COND / PDIFF ADJUST=BON; 
	format TX_COND tx_condfmt.;
	title 'Month 6: Model 1: Primary ANCOVA BEH_ARS = TX_COND ps_BEH_ARS';
run;
quit;


proc glm data=cokeM6;    
	class TX_COND GENDER;
	model BEH_ARS = TX_COND|GENDER ps_BEH_ARS / SOLUTION;
	/* Slice interaction to see TX_CODE effect within Genders */
	LSMEANS TX_COND*GENDER / SLICE=GENDER PDIFF ADJUST=TUKEY;        
	/* Estimate specific contrast*/
	ESTIMATE 'IDC+GDC versus GDC' TX_COND 1 0 0 -1;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month6: Model 2: 2-way Interaction BEH_ARS = TX_COND*GENDER ps_BEH_ARS';
run;
quit;

proc glm data=cokeM6;
	class TX_COND GENDER RACE;
	model BEH_ARS = TX_COND|GENDER|RACE ps_BEH_ARS / SOLUTION;
	LSMEANS TX_COND*GENDER*RACE / PDIFF=ALL;
	lsmeans TX_COND*GENDER*RACE / slice=GENDER*RACE;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month6: Model 3: 3-way Interaction BEH_ARS = TX_COND*GENDER*RACE ps_BEH_ARS';
run;
quit;

/***********************/
/* Mixed effects model */
/***********************/

data cokeM3_6;
 	set mylib.finalcoke;
	if MONTH=3 or MONTH=6;
run;

proc mixed data=cokeM3_6;
    CLASS TX_COND PATNO MONTH;
    MODEL BEH_ARS = TX_COND MONTH TX_COND*MONTH ps_BEH_ARS / SOLUTION;
    REPEATED MONTH / SUBJECT=PATNO TYPE=AR(1); /* Auto-regressive covariance */
    LSMEANS TX_COND*MONTH / PDIFF;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
    TITLE 'Mixed Effects Model';
run;



proc glm data=cokeM3 PLOTS=DIAGNOSTICS;
	class TX_COND GENDER RACE;
	/* MODEL: Main Effects + Covariate + Interaction */
    model BEH_ARS = TX_COND|GENDER ps_BEH_ARS / SOLUTION EFFECTSIZE;

	/* 1. LSMEANS with Plots to Visualize Interactions */
	/* This visualizes if lines cross (Interaction) */
	LSMEANS TX_COND*GENDER / PDIFF ADJUST=TUKEY PLOT=MEANPLOT(SLICEBY=GENDER JOIN);

	/* 2. SLICE statement to decompose interaction */
	/* Tests "Is there a Tx effect specifically for Men? Specifically for Women?" */
	LSMEANS TX_COND*GENDER / SLICE=GENDER;

	/* 3. DEEP DIVE: ESTIMATES & CONTRASTS (Determining the Driver) */
    /* Hypothesis A: Does IDC (Tx 1) outperform the average of all other treatments? */
    /* Coefficients: 3 for IDC, -1 for others. Divisor=3 averages the others. */
    ESTIMATE 'Driver: IDC vs All Others' TX_COND 3 -1 -1 -1 / DIVISOR=3;

    /* Hypothesis B: Do Professional Psychotherapies (CT & SE) outperform simple GDC? */
    /* Compares Tx 2 & 3 vs Tx 4 */
    ESTIMATE 'Psychotherapy vs GDC Control' TX_COND 0 1 1 -2 / DIVISOR=2;

    /* Hypothesis C: Does adding ANYTHING to GDC help? (IDC/CT/SE vs GDC) */
    CONTRAST 'Any Add-on vs GDC Alone' TX_COND 1 1 1 -3;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month 3: Additional GLM Analysis';
run;
quit;



proc glm data=cokeM6 PLOTS=DIAGNOSTICS;
	class TX_COND GENDER RACE;
	/* MODEL: Main Effects + Covariate + Interaction */
    model BEH_ARS = TX_COND|GENDER ps_BEH_ARS / SOLUTION EFFECTSIZE;

	/* 1. LSMEANS with Plots to Visualize Interactions */
	/* This visualizes if lines cross (Interaction) */
	LSMEANS TX_COND*GENDER / PDIFF ADJUST=TUKEY PLOT=MEANPLOT(SLICEBY=GENDER JOIN);

	/* 2. SLICE statement to decompose interaction */
	/* Tests "Is there a Tx effect specifically for Men? Specifically for Women?" */
	LSMEANS TX_COND*GENDER / SLICE=GENDER;

	/* 3. DEEP DIVE: ESTIMATES & CONTRASTS (Determining the Driver) */
    /* Hypothesis A: Does IDC (Tx 1) outperform the average of all other treatments? */
    /* Coefficients: 3 for IDC, -1 for others. Divisor=3 averages the others. */
    ESTIMATE 'Driver: IDC vs All Others' TX_COND 3 -1 -1 -1 / DIVISOR=3;

    /* Hypothesis B: Do Professional Psychotherapies (CT & SE) outperform simple GDC? */
    /* Compares Tx 2 & 3 vs Tx 4 */
    ESTIMATE 'Psychotherapy vs GDC Control' TX_COND 0 1 1 -2 / DIVISOR=2;

    /* Hypothesis C: Does adding ANYTHING to GDC help? (IDC/CT/SE vs GDC) */
    CONTRAST 'Any Add-on vs GDC Alone' TX_COND 1 1 1 -3;
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
	title 'Month 6: Additional GLM Analysis';
run;
quit;











/*
proc glm data=cokeM6;
    class TX_COND;
    model BEH_ARS = TX_COND;
    contrast 'TX_COND IDC vs Others' TX_COND 3 -1 -1 -1; 
    contrast 'TX_COND IDC vs TX_COND CT'   TX_COND 1 -1  0  0;
    title 'Month 6: Linear Contrasts';
run;
quit;



proc mixed data=cokeM3_6;
    class TX_COND GENDER PATNO MONTH; 
    

    model BEH_ARS = TX_COND MONTH TX_COND*MONTH GENDER / solution;
    

    random intercept / subject=PATNO; 
    
    lsmeans TX_COND*MONTH / pdiff;
    title 'Mixed Model (Month 3 & 6)';
run;

proc sgplot data=cokeM3_6;
    vline MONTH / group=tx_cond 
                   response=BEH_ARS 
                   stat=mean 
                   markers;
    xaxis label="Month";
    yaxis label="Mean BEH_ARS";
    keylegend / title="Treatment";
	format tx_cond tx_condfmt.;
	title "Treatment by Month";
run;

proc mixed data=cokeM3_6;
    class PATNO tx_cond month;
    model BEH_ARS = tx_cond month tx_cond*month / s;
    random PATNO; 
    lsmeans tx_cond / pdiff cl;
    lsmeans tx_cond*month / slice=month; 
	format tx_cond tx_condfmt.;
	title "BEH_ARS Mixed Model";
run;


proc glm data=cokeM6 plots=diagnostics; 
    class TX_COND GENDER;
    model BEH_ARS = TX_COND|GENDER;
    lsmeans TX_COND*GENDER / pdiff adjust=tukey plots=interaction(sliceby=GENDER); 
    title 'Month 6: Two-Way ANOVA with Interaction Plot';
run; quit;


*/

PROC FORMAT;
    VALUE tx_fmt
       1='IDC+GDC' 2='CT+GDC' 3='SE+GDC' 4='GDC';
    VALUE month_fmt
        0 = 'Intake'
        1 = '1'
        2 = '2'
        3 = '3'
        4 = '4'
        5 = '5'
        6 = '6'
        9 = '9'
        12 = '12';
RUN;

/* 3. PREPARE DATA */
DATA plot_data;
    SET mylib.finalcoke;
    
    /* Create Binary Indicator: 1 if used cocaine, 0 if abstinent */
    /* Assuming DAYSUSED is the days used in past 30 days */
    if DAYSUSED > 0 then AnyCocaine = 1;
    else if DAYSUSED = 0 then AnyCocaine = 0;
    
    /* Apply Formats */
    FORMAT TX_COND tx_fmt. MONTH month_fmt.;
RUN;

/* 4. CALCULATE PROPORTIONS (MEANS) */
PROC MEANS DATA=plot_data NWAY NOPRINT;
    CLASS MONTH TX_COND;
    VAR AnyCocaine;
    OUTPUT OUT=graph_stats MEAN=PropUsed;
RUN;

/* 5. GENERATE GRAPH */
TITLE 'Proportion of Patients Using Cocaine';
PROC SGPLOT DATA=graph_stats;
    /* Line Plot with Markers */
    STYLEATTRS DATASYMBOLS=(DiamondFilled SquareFilled TriangleFilled CircleFilled)
               DATACONTRASTCOLORS=(Black Red Yellow Green); /* Grayscale style */
               
    SERIES X=MONTH Y=PropUsed / GROUP=TX_COND 
                               MARKERS 
                               MARKERATTRS=(SIZE=9)
                               LINEATTRS=(THICKNESS=2);

    /* X-Axis Configuration */
    XAXIS LABEL="Months" 
          VALUES=(0 1 2 3 4 5 6) 
          OFFSETMIN=0.05 OFFSETMAX=0.05;

    /* Y-Axis Configuration */
    YAXIS LABEL="Proportion of Patients Using Cocaine" 
          MIN=0.3 MAX=1.0 
          VALUES=(0.3 to 1.0 by 0.1)
          GRID; /* Optional: Remove if exact white background needed */

    /* Legend Configuration */
    KEYLEGEND / LOCATION=INSIDE POSITION=TOPRIGHT ACROSS=1 NOBORDER;
RUN;
TITLE;

/* 4. CALCULATE MEANS FOR BEH_ARS */
PROC MEANS DATA=plot_data NWAY NOPRINT;
    CLASS MONTH TX_COND;
    VAR BEH_ARS;
    OUTPUT OUT=graph_stats MEAN=MeanBEH;
RUN;

/* 5. GENERATE GRAPH */
TITLE 'Mean BEH_ARS by Month';
PROC SGPLOT DATA=graph_stats;
    /* Line Plot with Markers */
    STYLEATTRS DATASYMBOLS=(DiamondFilled SquareFilled TriangleFilled CircleFilled)
               DATACONTRASTCOLORS=(Black Red Yellow Green); /* Grayscale style */
               
    SERIES X=MONTH Y=MeanBEH / GROUP=TX_COND 
                               MARKERS 
                               MARKERATTRS=(SIZE=9)
                               LINEATTRS=(THICKNESS=2);

    /* X-Axis Configuration */
    XAXIS LABEL="Months" 
          VALUES=(0 1 2 3 4 5 6) 
          OFFSETMIN=0.05 OFFSETMAX=0.05;

    /* Y-Axis Configuration */
    /* Removed fixed 0.3-1.0 range as BEH_ARS is on a different scale (approx 40-60) */
    YAXIS LABEL="Mean BEH_ARS Score" 
	      MIN=40 MAX=60
		  VALUES=(40 to 60 by 5)
          GRID; 

    /* Legend Configuration */
    KEYLEGEND / LOCATION=INSIDE POSITION=TOPRIGHT ACROSS=1 NOBORDER;
RUN;
TITLE;



PROC SORT DATA=cokeM3_6;
    BY TX_COND MONTH;
RUN;


/* CLUSTERED SIDE-BY-SIDE PLOTS */
ODS GRAPHICS ON / WIDTH=7in HEIGHT=6in;

TITLE 'Comparison of BEH_ARS by Intervention: Month 3 vs Month 6';
PROC SGPLOT DATA=cokeM3_6;
    /* VBOX with Grouping and Clustering */
    /* CATEGORY=TX_COND puts Interventation on the X-axis */
    /* GROUP=MONTH splits each Interventation into Month 3 and 6 */
    /* GROUPDISPLAY=CLUSTER places them side-by-side instead of stacked */
    VBOX BEH_ARS / CATEGORY=TX_COND 
                   GROUP=MONTH 
                   GROUPDISPLAY=CLUSTER
                   MEANATTRS=(SYMBOL=DiamondFilled)
                   OUTLIERATTRS=(SYMBOL=Circle); 
                   
    STYLEATTRS DATACONTRASTCOLORS=(CornflowerBlue Firebrick) 
               DATACOLORS=(AliceBlue MistyRose);

    YAXIS LABEL="BEH_ARS" GRID;
    XAXIS LABEL="Intervention";
    
    KEYLEGEND / POSITION=TOPRIGHT NOBORDER TITLE="Month";
	FORMAT TX_COND tx_fmt. MONTH month_fmt.;
RUN;
TITLE; 

ODS GRAPHICS OFF;


ODS GRAPHICS ON / WIDTH=8in HEIGHT=6in;

TITLE 'BEH_ARS: Month 3 vs Month 6 (Interventions Combined)';
PROC SGPLOT DATA=cokeM3_6;
    /* VBOX with CATEGORY=MONTH places the two timepoints on the X-axis */
    VBOX BEH_ARS / CATEGORY=MONTH 
                   MEANATTRS=(SYMBOL=DiamondFilled SIZE=10)
                   OUTLIERATTRS=(SYMBOL=Circle);
                   
    STYLEATTRS DATACONTRASTCOLORS=(DarkBlue) 
               DATACOLORS=(CornflowerBlue);

    YAXIS LABEL="BEH_ARS" GRID;
    XAXIS LABEL="Month";
RUN;
TITLE; 

ODS GRAPHICS OFF;


PROC MEANS DATA=cokeM3_6 NWAY NOPRINT;
    CLASS MONTH TX_COND GENDER;
    VAR BEH_ARS;
    OUTPUT OUT=means_data MEAN=MeanScore;
RUN;

PROC SORT DATA=means_data;
    BY TX_COND MONTH GENDER;
RUN;


/* 
 Interaction Plot: Intervention x Gender over Month

 Check if Intervention effect depends on Gender.
 Crossing lines (Interaction) vs Parallel lines (No Interaction).

*/

TITLE 'Interaction Plot: Intervention x Gender over Month';
PROC SGPANEL DATA=means_data;
    /* Panel by Treatment to see the trajectory for each condition side-by-side */
    PANELBY TX_COND / COLUMNS=4 NOVARNAME;
    
    /* X-Axis = Month (Timepoint), Group = Gender */
    /* This visualizes the Gender x Time interaction within each Intervention */
    SERIES X=MONTH Y=MeanScore / GROUP=GENDER 
                                 MARKERS 
                                 MARKERATTRS=(SIZE=10 SYMBOL=CircleFilled)
                                 LINEATTRS=(THICKNESS=3);
                                   
    ROWAXIS LABEL="Mean BEH_ARS" GRID;
    COLAXIS LABEL="Month" MIN=3 MAX=6 VALUES=(3 to 6 by 3);
    KEYLEGEND / TITLE="Gender";
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
RUN;
TITLE; 



PROC MEANS DATA=cokeM3_6 NWAY NOPRINT;
    CLASS MONTH TX_COND RACE;
    VAR BEH_ARS;
    OUTPUT OUT=means_data2 MEAN=MeanScore;
RUN;

PROC SORT DATA=means_data2;
    BY TX_COND MONTH RACE;
RUN;


/* 
 Interaction Plot: Intervention x Race over Month

 Check if Intervention effect depends on Race.
 Crossing lines (Interaction) vs Parallel lines (No Interaction).

*/

TITLE 'Interaction Plot: Intervention x Race over Month';
PROC SGPANEL DATA=means_data2;
    /* Panel by Treatment to see the trajectory for each condition side-by-side */
    PANELBY TX_COND / COLUMNS=4 NOVARNAME;
    
    /* X-Axis = Month, Group = Race */
    /* This visualizes the Race x Time interaction within each Intervention */
    SERIES X=MONTH Y=MeanScore / GROUP=RACE 
                                 MARKERS 
                                 MARKERATTRS=(SIZE=10 SYMBOL=CircleFilled)
                                 LINEATTRS=(THICKNESS=3);
                                   
    ROWAXIS LABEL="Mean BEH_ARS" GRID;
    COLAXIS LABEL="Month" MIN=3 MAX=6 VALUES=(3 to 6 by 3);
    KEYLEGEND / TITLE="Race";
	format GENDER genderfmt.;
	format RACE racefmt.;
	format TX_COND tx_condfmt.;
RUN;
TITLE;

SAS to Python

The equivalences the notebook relies on.
SAS	Python	Note
`PROC GLM ... MODEL y = A x`	`smf.ols("y ~ C(A) + x")`	ANCOVA. Both drop rows missing any model term
`CLASS A`	`C(A, Treatment(ref))`	SAS and statsmodels pick different default references; set it explicitly or contrasts will not match
`LSMEANS A / PDIFF`	`get_prediction` on a covariate-mean grid	Adjusted means
`MEANS A / HOVTEST=LEVENE`	`scipy.stats.levene(..., center="median")`	SAS defaults to the median form
`PROC UNIVARIATE NORMAL`	`scipy.stats.shapiro`	Test the residuals of the model you report
`PROC MIXED / TYPE=AR(1)`	`MixedLM`, or GEE	statsmodels has no direct AR(1) repeated-measures equivalent
`OUTPUT OUT=d R=residuals`	`fit.resid`

Rebuilding this site

Stages that need the trial data write aggregate results into tools/derived/, which is committed, so the pages rebuild on any checkout without it.

From the repository rootshell

bash tools/build.sh          # rebuild the pages from committed data
bash tools/build.sh --all    # also recompute from the trial data