Thanks Kurt, I am copying the macro in SAS code window below. I used this syntax to CALL it; %macro itsa(model=, dataset=, outcome=, time=, interrupt=, lag=);
%if (&model=sitsa) %then %do;
* Create SITSA dummy variables;
data sitsa_vars;
set &dataset;
t=_n_;
if &time<&interrupt then x=0;
else x=1;
if &time<&interrupt then tx=0;
else tx+1;
keep &time &outcome t x tx;
run;
* Compute lag+1 for proc model;
%let lagl=%eval(&lag + 1);
* If lag is set to 0 or missing then conduct SITSA with OLS regression;
* If lag>0 then conduct SITSA with OLS regression + Newey-West standard errors;
%if (&lag<1) %then %do;
ods output ParameterEstimates=sitsa_stat TestResults=sitsa_stat2 CovB=sitsa_covar;
proc model data=sitsa_vars;
parms b0 b1 b2 b3;
&outcome = b0+(b1*t)+(b2*x)+(b3*tx);
fit &outcome / covb;
test b1+b3;
run; quit;
%end;
%else %do;
ods output ParameterEstimates=sitsa_stat TestResults=sitsa_stat2 CovB=sitsa_covar;
proc model data=sitsa_vars;
parms b0 b1 b2 b3;
&outcome = b0+(b1*t)+(b2*x)+(b3*tx);
fit &outcome / covb gmm kernel=(bart,&lagl,0) vardef=n;
test b1+b3;
run; quit;
%end;
* Produce tables for model coefficients;
data sitsa_stat;
set sitsa_stat;
keep parameter estimate stderr probt;
rename probt=prob;
run;
data sitsa_stat2;
set sitsa_stat2;
keep label probchisq;
rename label=parameter probchisq=prob;
run;
* Append stats from outputs;
data sitsa_stat;
set sitsa_stat sitsa_stat2;
run;
* Calculate and save sum of b1+b3 coefficients;
proc means data=sitsa_stat(where=(parameter='b1' or parameter='b3')) sum noprint;
output out=b1b3_est sum= / autoname;
var estimate;
run;
* Save sum to macro variable;
data _null_;
set b1b3_est;
call symputx("Sum", Estimate_Sum);
run;
* Add sum as b1+b3 parameter estimate;
data sitsa_stat;
set sitsa_stat;
if parameter='b1+b3' then estimate=∑
run;
* Save covariance to macro variable;
data _null_;
set sitsa_covar;
if parameter='b3' then call symputx("Covar", b1);
run;
* Save b1 and b3 standard errors as macro variables;
data _null_;
set sitsa_stat;
if parameter='b1' then call symputx("SE_1", StdErr);
if parameter='b3' then call symputx("SE_3", StdErr);
run;
* Calculate b1+b3 standard error;
%let sum_se=%sysfunc(sqrt(&SE_1.**2+&SE_3.**2+2*&Covar.));
* Add b1+b3 standard error;
* Calculate 95% confidence intervals for all results;
data sitsa_stat;
set sitsa_stat;
if parameter='b1+b3' then StdErr=&sum_se;
CI_LL=(estimate-1.96*StdErr);
CI_UL=(estimate+1.96*StdErr);
run;
* Calculate residual/predicted values;
data _null_;
set sitsa_stat;
if parameter='b0' then call symputx("b0", Estimate);
if parameter='b1' then call symputx("b1", Estimate);
if parameter='b2' then call symputx("b2", Estimate);
if parameter='b3' then call symputx("b3", Estimate);
run;
data sitsa_vars;
set sitsa_vars;
pred=&b0+(&b1*t)+(&b2*x)+(&b3*tx);
res=(&outcome-pred);
run;
* Print coefficient estimates;
proc print data=sitsa_stat;
id parameter;
run;
* Define plot characteristics;
symbol1 color=black line=1 interpol=join value=star;
symbol2 color=black line=1 interpol=join value=dot;
axis1 label=('Time') offset=(1,1);
axis2 label=('Outcome') offset=(0,0);
legend1 label=('Phase:') value=('Pre-Interruption' 'Post-Interruption');
pattern1 color=gray;
pattern2 color=black;
* Produce 2x plots;
proc gplot data=sitsa_vars;
plot &outcome*t = x / haxis=axis1 vaxis=axis2 cvref=black legend=legend1;
run; quit;
proc gplot data=sitsa_vars;
plot &outcome*t = x / areas=2 haxis=axis1 vaxis=axis2 cvref=black legend=legend1;
run; quit;
* Delete additional datasets;
proc datasets library=work noprint;
delete sitsa_stat sitsa_stat2 b1b3_est sitsa_covar;
run; quit;
%end;
%mend; %INCLUDE itsa; How do I execute it. You are right, I am not familiar with macro language but using this macro is the best way of analysing the data I have, so I am trying to utilise given resources as best as I can. Thanks for all your help. S
... View more