*--------------------------------------------------------------*
*                                                              *
*  SMSUB 3.1                                                   *
*                                                              *
*--------------------------------------------------------------*;
%macro smsub(
data=,
statistics=,
strata=-99,
cluster=-99,
popsize=-99,
alpha=0.05,
weight=,
tables=,
var=99,
ratio=,
contrast=,
subpop=-99,
format=,
output=table=_tables_ summary=_summary_ contrast=_contrast_,
colwidth=26,
options=,
title=The SAS System);

options nocenter nonumber nodate nonotes nostimer formdlim=' '
        formchar="|----|+|---+=|-/\<>*" ps=max ls=94;
options &options;

 *---check for iml license---*;
data _null_;
imlproduct=sysprod('IML');
call symput('ERR1',imlproduct);
run;
%if (&ERR1^=1) %then %do;
 %put; %put; %put;
 %put ERROR:  The SMSUB macro requires SAS/IML to be licensed.;
 %put;
 %goto exit;
%end;

%let usertitle=%sysfunc(compress(&title,"'"));
title;
 *---initialize message indicators---*;
%let ERR101=0; %let ERR102=0; %let ERR103=0; %let ERR104=0; %let ERR105=0;
%let ERR106=0; %let ERR107=0; %let ERR108=0; %let ERR109=0;
%let ERR201=0; %let ERR202=0; %let ERR203=0; %let ERR204=0;
%let ERR301=0;

 *---err check (1)---*;
%let ERR103= %index(&var,-);
%if (&ERR103) %then %do;
 %put; %put; %put;
 %put ERROR:  Numbered range lists are not allowed in the VAR= parameter.;
 %put;
 %goto exit;
%end;

 *---parameter indicators---*;
%if (&var ne 99) %then %let ivar=1; %else %let ivar=0;
%if (%length(&ratio)) %then %let iratio=1; %else %let iratio=0;
%if (&strata ne -99) %then %let istrata=1; %else %let istrata=0;
%if (&cluster ne -99) %then %let icluster=1; %else %let icluster=0;
%if (&popsize ne -99) %then %let ipopsize=1; %else %let ipopsize=0;
%if (&icluster and &ipopsize) %then %let iclupop=1; %else %let iclupop=0;
%if (%length(&tables)) %then %let itables=1; %else %let itables=0;
%if (%length(&statistics)) %then %let istatistics=1; %else %let istatistics=0;
%if (%length(&contrast)) %then %let icontrast=1; %else %let icontrast=0;
%if ("&subpop" ne "-99") %then %let isubpop=1; %else %let isubpop=0;

 *---syntax check---*;
%if (%length(&data)=0) or (%length(&weight)=0) or (&ivar=0 and &iratio=0)
  %then %let ERR101=1;
%if (&ivar and &iratio) %then %let ERR102=1;
%let ERR104=%length(%scan(&strata,2));
%let ERR105=%length(%scan(&cluster,2));

 *---process statistics---*;
%if (&istatistics=0) %then %do;
 %if (&ivar=1 and &icontrast=0) %then %let statistics=nobs mean semean clmean;
 %else %if (&ivar=1 and &icontrast=1) %then
                                 %let statistics=nobs mean semean tmean pmean;
 %else %if (&iratio=1 and &icontrast=0) %then
                                   %let statistics=nobs ratio seratio clratio;
 %else %if (&iratio=1 and &icontrast=1) %then
                             %let statistics=nobs ratio seratio tratio pratio;
 %end;
%statistics;

 *---process contrast---*;
%let coef1=-99; %let coef2=-99; %let coef3=-99; %let coef4=-99;
%let coef5=-99; %let coef6=-99; %let coef7=-99; %let coef8=-99;
%let coef9=-99; %let coef10=-99;
%if (&icontrast=0) %then %do;
 %let cvar=-99;
 %let outcontrastlabel=;
 %let outcontrastindex=;
%end;
%else %do;
 %contrast;
 %let outcontrastlabel=CONTRASTLABEL;
 %let outcontrastindex=CONTRASTINDEX;
%end;

 *---process subpop---*;
%let inoperator=0;
%let inelements=-99;
%if (&isubpop=0) %then %do;
 %let subpopvar=one;
 %let subpop2=one=1;
 %let subpoplabel=;
 %let subpopwsumlabel=;
 %let subpopobslabel=;
 %let outsubpop=;
 %let outsubpopobs=;
 %let outsubpopwsum=;
%end;
%else %do;
 %subpop;
 %let outsubpop=subpop;
 %let outsubpopobs=subpopobs;
 %let outsubpopwsum=subpopwsum;
 %let subpopwsumlabel=subpopwsum='Subpopulation Weight Sum';
 %let subpopobslabel=subpopobs='Subpopulation Observations';
%end;

 *---err checks (2)---*;
%if (&ERR101 or &ERR102 or &ERR104 or &ERR105) %then %do;
 %put; %put; %put;
 %if (&ERR101) %then %do;
  %put ERROR:  Missing one or more required parameters (DATA=, WEIGHT=, and;
  %put %str(        )either VAR= or RATIO= ).;
  %end;
 %if (&ERR102) %then %put
  ERROR:  The VAR= and RATIO= parameters cannot be used simultaneously.;
 %if (&ERR104) %then %put
  ERROR:  You may only specify one variable in the STRATA= parameter.;
 %if (&ERR105) %then %put
  ERROR:  You may only specify one variable in the CLUSTER= parameter.;
 %put;
 %goto exit;
%end;
%if (&ERR201 or &ERR202 or &ERR203 or &ERR204 OR &ERR301) %then %do;
 %put; %put; %put;
 %if (&ERR201) %then %do;
  %put ERROR:  The following statistics are not available with the VAR=
parameter:;
  %put %str(          )NUMTOTAL DENTOTAL RATIO SERATIO CLRATIO CVRATIO.;
  %end;
 %if (&ERR202) %then %do;
  %put ERROR:  The following statistics are not available with the RATIO=
parameter:;
%put %str(          )MEAN SEMEAN CLMEAN CVMEAN TOTAL SETOTAL CLTOTAL CVTOTAL.;
  %end;
 %if (&ERR203) %then %do;
%put ERROR:  The following specified statistic-keyword(s) are not recognized:;
%put %str(          )&unknown..;
  %end;
 %if (&ERR204) %then %do;
%put ERROR:  The following statistics require the CONTRAST= parameter:;
%put %str(          )TMEAN PMEAN TRATIO PRATIO;
  %end;
 %if (&ERR301) %then %put
  ERROR:  Syntax error in the SUBPOP= parameter.;
 %put;
 %goto exit;
%end;

 *---process var,ratio ---*;
%if (&ivar) %then %do;
 %let indextype=varindex;
 %let statvar=variable;
 %let nratio=-99;
 %let topvar=-99; %let topequal=-99; %let topvalue=-99;
 %let botvar=-99; %let botequal=-99; %let botvalue=-99;
%end;
%else %do;
 %let indextype=ratioindex;
 %let statvar=rationame;
 %ratio;
%end;

 *---process output---*;
%let outtable= _tables_;
%let outsummary=_summary_;
%let outcontrast=_contrast_;
%let output=%upcase(&output);
 *---first dataset name---*;
%let equal= %index(&output,=);
%let leftside=%upcase(%substr(&output,1,&equal-1));
%let rest=%substr(&output,&equal+1);
%let rightside=%scan(&rest,1);
%if (&leftside=TABLES) %then         %let outtable=&rightside;
%else %if (&leftside=SUMMARY) %then  %let outsummary=&rightside;
%else %if (&leftside=CONTRAST) %then %let outcontrast=&rightside;
 *---second dataset name---*;
%let rest=&rest z_z_z;
%let nextword=%scan(&rest,2);
%if (&nextword ne z_z_z) %then %do;
 %let pnextword=%index(&rest,&nextword);
 %let truncate=%substr(&rest,&pnextword);
 %let equal= %index(&truncate,=);
 %let leftside=%upcase(%substr(&truncate,1,&equal-1));
 %let rest=%substr(&truncate,&equal+1);
 %let rightside=%scan(&rest,1);
 %if (&leftside=TABLES) %then         %let outtable=&rightside;
 %else %if (&leftside=SUMMARY) %then  %let outsummary=&rightside;
 %else %if (&leftside=CONTRAST) %then %let outcontrast=&rightside;
%end;
 *---third dataset name---*;
%let nextword=%scan(&rest,2);
%if (&nextword ne z_z_z) %then %do;
 %let pnextword=%index(&rest,&nextword);
 %let truncate=%substr(&rest,&pnextword);
 %let equal= %index(&truncate,=);
 %let leftside=%upcase(%substr(&truncate,1,&equal-1));
 %let rest=%substr(&truncate,&equal+1);
 %let pz_z_z=%index(&rest,z_z_z);
 %let rightside=%substr(&rest,1,&pz_z_z-1);
 %if (&leftside=TABLES) %then         %let outtable=&rightside;
 %else %if (&leftside=SUMMARY) %then  %let outsummary=&rightside;
 %else %if (&leftside=CONTRAST) %then %let outcontrast=&rightside;
%end;

 *---process table options---*;
%let style=BOX;
%let tables=%sysfunc(tranwrd(&tables,*,#));
%let tables=%sysfunc(tranwrd(&tables,*,#));
%let tables=%sysfunc(tranwrd(&tables,%str(   #),#));
%let tables=%sysfunc(tranwrd(&tables,%str(  #),#));
%let tables=%sysfunc(tranwrd(&tables,%str( #),#));
%let tables=%sysfunc(tranwrd(&tables,%str(#   ),#));
%let tables=%sysfunc(tranwrd(&tables,%str(#  ),#));
%let tables=%sysfunc(tranwrd(&tables,%str(# ),#));
%let slash=%index(&tables,/);
%if (&slash) %then %do;
 %if (&slash=1)  %then %let tablerequests=;
 %else %let tablerequests=%substr(&tables,1,&slash-1);
 %let tableoptions=%upcase(%substr(&tables,&slash+1));
 %let levelmiss=%index(&tableoptions,LEVELMISS);
 %let nomarg=%index(&tableoptions,NOMARG);
 %let exclmissmarg=%index(&tableoptions,EXCLMISSMARG);
 %if (&exclmissmarg and (&levelmiss or &nomarg)) %then %let ERR106=1;
 %if (%index(&tableoptions,STYLE)) %then %do;
  %let equal=%index(&tableoptions,=);
  %let rest=%substr(&tableoptions,&equal+1);
  %let style=%scan(&rest,1);
 %end;
%end;
%else %do;
 %let tablerequests=&tables;
 %let levelmiss=0; %let nomarg=0; %let exclmissmarg=0;
%end;
 *--err check (3)---*;
%if (&ERR106) %then %do;
 %put; %put; %put;
  %put ERROR:  The EXCLMISSMARG option cannot be used with the LEVELMISS or;
  %put %str(        )NOMARG options.;
 %put;
 %goto exit;
%end;
 *---err109---*;
%if (&istrata)   %then %let strata_chk=&strata; %else %let strata_chk=;
%if (&icluster)  %then %let cluster_chk=&cluster; %else %let cluster_chk=;
%if (&ipopsize) %then %let popsize_chk=&popsize; %else %let popsize_chk=;
%if (&ivar) %then %let var_chk=&var; %else %let var_chk=;
%if (&iratio) %then %do;
 %let topvar_chk=&topvar;
 %let botvar_chk=&botvar;
 %end;
%else %do;
 %let topvar_chk=;
 %let botvar_chk=;
%end;
%let tablevars=%upcase(%sysfunc(tranwrd(&tablerequests,#,%str( ))));
%let tablevars=%sysfunc(tranwrd(&tablevars,_OVERALL_,%str( )));
%if (&isubpop) %then %let subpopvar_chk=&subpopvar; %else %let subpopvar_chk=;
%if (&icontrast) %then %let cvar_chk=&cvar; %else %let cvar_chk=;
proc contents data=&data(keep= &strata_chk &cluster_chk &popsize_chk &weight
  &var_chk &topvar_chk &botvar_chk &tablevars &subpopvar_chk &cvar_chk)
  out=_vars1(keep=type) noprint; run;
%if (&syserr) %then %do;
 %put; %put; %put;
 %put ERROR:  Input Error.;
 %put;
 %goto exit;
%end;
data _vars2;
set _vars1;
if type=2 then call symput('ERR109','1');
data _null_;
 *---err check (4)---*;
%if (&ERR109) %then %do;
 %put; %put; %put;
 %put ERROR:  Character variables are not allowed in the parameters.;
 %put;
 %goto exit;
%end;

 *---process tablerequests---*;
%if %length(&tablerequests)=0 %then %do;
 %let tabletype=0; %let tableindex=1; %let _overall_=1;
 %let numsubvars=1;
 %let subgrp1=_overall_; %let subgrp2=-99; %let subgrp3=-99;
 %let grouplist=_OVERALL_; %let tablerequest=Overall;
 %let printtablerequest=_overall_; %let outsubgrp1=&subgrp1;
 %let outsubgrp2=ignore; %let outsubgrp3=ignore;
 %let s1fmt=_commona; %let s2fmt=_commonb; %let s3fmt=_commonc;
 %let s1fmtcodes=; %let s2fmtcodes=; %let s3fmtcodes=;
 %COMPUTE;
 %if (&syserr) %then %do;
  %put; %put; %put;
  %put ERROR:  Input error.;
  %put;
  %goto exit;
 %end;
 data _null_;
 set _err;
 call symput('ERR401',ERR401);
 call symput('ERR402',ERR402);
 call symput('ERR403',ERR403);
 call symput('ERR404',ERR404);
 call symput('ERRstr',trim(left(ERRstr)));
 call symput('ERR405',ERR405);
 call symput('ERR406',ERR406);
 call symput('ERR407',ERR407);
  *---err checks (7)---*;
 data _null_;
 %if (&ERR401 or &ERR402 or &ERR403 or &ERR404 or &ERR405 or &ERR406 or
      &ERR407) %then %do;
  %if (&icluster) %then %let errpsu=cluster (&cluster);
  %else                 %let errpsu=observation;
  %put; %put; %put;
  %if (&ERR401) %then %put
   ERROR:  The POPSIZE= variable can only have one value per stratum.;
  %if (&ERR402) %then %put
   ERROR:  Subgroup variables must be greater than zero.;
  %if (&ERR403) %then %put
   ERROR:  Contrast variables must be greater than zero.;
  %if (&ERR404) %then %do;
   %put ERROR:  Cannot compute variance.  Stratum &errstr has only one;
   %put %str(        )noncertainty &errpsu..;
   %end;
  %if (&ERR405) %then %do;
   %put ERROR:  The STRATA=, CLUSTER=, POPSIZE=, and WEIGHT= variables;
   %put %str(        )cannot have missing values.;
   %end;
  %if (&ERR406) %then %put
   ERROR:  The POPSIZE= variable cannot have values less than -1.;
  %if (&ERR407) %then %put
   ERROR:  The alpha value must be between 0.0001 and 0.9999.;
  %put;
  %goto exit;
 %end;
 %OUTPUT;
 %if (&syserr) %then %do;
  %put; %put; %put;
  %put ERROR:  Input error.;
  %put;
  %goto exit;
 %end;
 %if (&noprint=0) %then %PRINT;
%end;
%else %do;
 %let tableindex=1;
 %let tablerequest=%scan(&tablerequests,&tableindex);
 %do %while (%length(&tablerequest));
   *---err checks (5)---*;
  %let varcnt=1;
  %let tablevar=%scan(&tablerequest,1,#);
  %do %while(%length(&tablevar)>0);
   %if (&icontrast) %then %do;
    %let i=1;
    %let contvar=%scan(&cvar,1);
    %do %while(%length(&contvar));
     %if (&contvar=&tablevar) %then %let ERR108=1;
     %let i=%eval(&i+1);
     %let contvar=%scan(&cvar,&i);
    %end;
   %end;
   %let varcnt=%eval(&varcnt+1);
   %let tablevar=%scan(&tablerequest,%eval(&varcnt),#);
  %end;
  %if (&varcnt>4) %then %let ERR107=1;
   *--err checks (6)---*;
  %if (&ERR107 or &ERR108) %then %do;
   %put; %put; %put;
   %if (&ERR107) %then %put
    ERROR:  Table requests cannot cross more than 3 variables.;
   %if (&ERR108) %then %put
    ERROR:  The CONTRAST= variable cannot be in the TABLES= parameter.;
   %put;
   %goto exit;
  %end;
   *---run table---*;
  %TABLE;
  %COMPUTE;
  %if (&syserr) %then %do;
   %put; %put; %put;
   %put ERROR:  Input error.;
   %put;
   %goto exit;
  %end;
  data _null_;
  set _err;
  call symput('ERR401',ERR401);
  call symput('ERR402',ERR402);
  call symput('ERR403',ERR403);
  call symput('ERR404',ERR404);
  call symput('ERRstr',trim(left(ERRstr)));
  call symput('ERR405',ERR405);
  call symput('ERR406',ERR406);
  call symput('ERR407',ERR407);
   *---err checks (7)---*;
  data _null_;
  %if (&ERR401 or &ERR402 or &ERR403 or &ERR404 or &ERR405 or &ERR406 or
      &ERR407) %then %do;
   %if (&icluster) %then %let errpsu=cluster (&cluster);
   %else                 %let errpsu=observation;
   %put; %put; %put;
   %if (&ERR401) %then %put
    ERROR:  The POPSIZE= variable can only have one value per stratum.;
   %if (&ERR402) %then %put
    ERROR:  Subgroup variables must be greater than zero.;
   %if (&ERR403) %then %put
    ERROR:  Contrast variables must be greater than zero.;
   %if (&ERR404) %then %do;
    %put ERROR:  Cannot compute variance.  Stratum &errstr has only one;
    %put %str(        )noncertainty &errpsu..;
    %end;
   %if (&ERR405) %then %do;
    %put ERROR:  The STRATA=, CLUSTER=, POPSIZE=, and WEIGHT= variables;
    %put %str(        )cannot have missing values.;
    %end;
   %if (&ERR406) %then %put
    ERROR:  The POPSIZE= variable cannot have values less than -1.;
   %if (&ERR407) %then %put
    ERROR:  The alpha value must be between 0.0001 and 0.9999.;
   %put;
   %goto exit;
  %end;
  %OUTPUT;
  %if (&syserr) %then %do;
   %put; %put; %put;
   %put ERROR:  Input error.;
   %put;
   %goto exit;
  %end;
  %if (&noprint=0) %then %PRINT;
  %let tableindex=%eval(&tableindex+1);
  %let tablerequest=%scan(&tablerequests,&tableindex);
 %end;
%end;

%exit:
options notes;
data _null_;
file log; put;
put 'NOTE:  The SMSUB MACRO has completed.';
run;
title'The SAS System'; title2; title3; title4;
%mend smsub;

*--------------------------------------------------------*
*                                                        *
*    COMPUTE                                             *
*                                                        *
*--------------------------------------------------------*;
%macro compute;
proc iml;
start compute;
 *---initialize error codes---*;
ERR401=0; ERR402=0; ERR403=0; ERR404=0; ERR405=0; ERR406=0; ERR407=0;
 *---initialize warning codes---*;
warn101=0; warn102=0; warn103=0;
 *---initialize variance---*;
vtotal=0; vmean=0; vratio=0; npsu=0;

 *---err407---*;
alpha=&alpha;
if (alpha <= .0001) | (alpha > 0.9999) then ERR407=1;
if (ERR407) then goto skip;

 *---tableindex---*;
_tableindex=&tableindex;

*---------------------------*
* READ VARIABLES            *
*---------------------------*;
use &data;
if (&ivar) then do;
 read all var{&var} into var;
 topvar=-99;
 botvar=-99;
 nvar=ncol(var);
end;
else do;
 read all var{&topvar} into topvar;
 read all var{&botvar} into botvar;
 var=-99;
 nratio=ncol(topvar);
 topequal={&topequal};
 botequal={&botequal};
 topvalue={&topvalue};
 botvalue={&botvalue};
end;
read all var{&weight} into w;
nobsread=nrow(w);
if (ncol(loc(w<=0))) then izerowt=1; else izerowt=0;
if (&istrata)  then
 read all var{&strata} into strata;
else strata=j(nobsread,1,1);
if (&icluster) then
 read all var{&cluster} into cluster;
else cluster=t(1:nobsread);
if (&ipopsize) then
 read all var{&popsize} into popsize;
else popsize=j(nobsread,1,1e100);
if (&isubpop) then read all var{&subpopvar} into &subpopvar;
else  &subpopvar=j(nobsread,1,1);
inelements={&inelements};
ninelements=ncol(inelements);
if (&icontrast) then read all var{&cvar} into cvar;
else cvar=-99;

 *---check design variables---*;
if ncol(loc(strata=.)) then ERR405=1;
if ncol(loc(cluster=.)) then ERR405=1;
if ncol(loc(popsize=.)) then ERR405=1;
if ncol(loc(w=.)) then ERR405=1;
if (min(popsize)<-1) then ERR406=1;
if (ERR405 | ERR406) then goto skip;

 *---number of strata---*;
stratadesign=design(strata);
nstrata=ncol(stratadesign);

 *---unique strata---*;
uniqstr=unique(strata);

 *---err401---*;
do h=1 to nstrata;
 popchk= popsize[loc(stratadesign[,h]=1)];
 if (ncol(unique(popchk)) > 1) then ERR401=1;
end;
if (ERR401) then goto skip;

store strata cluster popsize var topvar botvar cvar;
free strata cluster popsize var topvar botvar cvar;

 *---subgrps---*;
if (&tabletype=0) then do;
 outs1levels=0;
 outs2levels=-2;
 outs3levels=-2;
end;
else if (&tabletype=1) then do;
 if (&_overall_=1) then subgrp1=j(nobsread,1,1);
 else read all var{&subgrp1} into subgrp1;
 if ncol(loc(subgrp1=.)) then mvsubgrp1=1; else mvsubgrp1=0;
 if (mvsubgrp1) then subgrp1[loc(subgrp1=.)]=99e15;
 s1levels=unique(subgrp1);
 if (min(s1levels) <=0) then ERR402=1;
 if (ERR402) then goto skip;
 ns1col=ncol(s1levels);
 if (mvsubgrp1 & &levelmiss=0) then do;
  ns1col=ns1col-1;
  s1levels=s1levels[,1:ns1col];
 end;
 if (&nomarg) then outs1levels=s1levels;
 else              outs1levels= 0 || s1levels;
 outs2levels=-2;
 outs3levels=-2;
end;
else if (&tabletype=2) then do;
 read all var{&subgrp1} into subgrp1;
 read all var{&subgrp2} into subgrp2;
 if ncol(loc(subgrp1=.)) then mvsubgrp1=1; else mvsubgrp1=0;
 if ncol(loc(subgrp2=.)) then mvsubgrp2=1; else mvsubgrp2=0;
 if (mvsubgrp1) then subgrp1[loc(subgrp1=.)]=99e15;
 if (mvsubgrp2) then subgrp2[loc(subgrp2=.)]=99e15;
 s1levels=unique(subgrp1);
 s2levels=unique(subgrp2);
 if (min(s1levels)<=0 | min(s2levels)<=0) then ERR402=1;
 if (ERR402) then goto skip;
 ns1col=ncol(s1levels);
 ns2col=ncol(s2levels);
 if (mvsubgrp1 & &levelmiss=0) then do;
  ns1col=ns1col-1;
  s1levels=s1levels[,1:ns1col];
 end;
 if (mvsubgrp2 & &levelmiss=0) then do;
  ns2col=ns2col-1;
  s2levels=s2levels[,1:ns2col];
 end;
 if (&nomarg=0) then outs1levels=repeat(0,1,ns2col+1);
 do s1=1 to ns1col;
  if (&nomarg=0) then outs1levels= outs1levels ||
repeat(s1levels[s1],1,ns2col+1);
  else                outs1levels= outs1levels ||
repeat(s1levels[s1],1,ns2col);
 end;
 if (&nomarg=0) then outs2levels= repeat(0 || s2levels,1,ns1col+1);
 else                outs2levels= repeat(s2levels,1,ns1col);
 outs3levels= -2;
end;
else if (&tabletype=3) then do;
 read all var{&subgrp1} into subgrp1;
 read all var{&subgrp2} into subgrp2;
 read all var{&subgrp3} into subgrp3;
 if ncol(loc(subgrp1=.)) then mvsubgrp1=1; else mvsubgrp1=0;
 if ncol(loc(subgrp2=.)) then mvsubgrp2=1; else mvsubgrp2=0;
 if ncol(loc(subgrp3=.)) then mvsubgrp3=1; else mvsubgrp3=0;
 if (mvsubgrp1) then subgrp1[loc(subgrp1=.)]=99e15;
 if (mvsubgrp2) then subgrp2[loc(subgrp2=.)]=99e15;
 if (mvsubgrp3) then subgrp3[loc(subgrp3=.)]=99e15;
 s1levels=unique(subgrp1);
 s2levels=unique(subgrp2);
 s3levels=unique(subgrp3);
 if (min(s1levels)<=0 | min(s2levels)<=0 | min(s3levels)<= 0) then ERR402=1;
 if (ERR402) then goto skip;
 ns1col=ncol(s1levels);
 ns2col=ncol(s2levels);
 ns3col=ncol(s3levels);
 if (mvsubgrp1 & &levelmiss=0) then do;
  ns1col=ns1col-1;
  s1levels=s1levels[,1:ns1col];
 end;
 if (mvsubgrp2 & &levelmiss=0) then do;
  ns2col=ns2col-1;
  s2levels=s2levels[,1:ns2col];
 end;
 if (mvsubgrp3 & &levelmiss=0) then do;
  ns3col=ns3col-1;
  s3levels=s3levels[,1:ns3col];
 end;
 if (&nomarg=0) then outs3levels=repeat(0,1,(ns1col+1)#(ns2col+1));
 do s3=1 to ns3col;
  if (&nomarg=0) then outs3levels= outs3levels ||
repeat(s3levels[s3],1,(ns1col+1)#(ns2col+1));
  else                outs3levels= outs3levels ||
repeat(s3levels[s3],1,ns1col#ns2col);
 end;
 if (&nomarg=0) then outs1levels= repeat(0,1,ns2col+1);
 do s1=1 to ns1col;
  if (&nomarg=0) then outs1levels= outs1levels ||
repeat(s1levels[s1],1,ns2col+1);
  else                outs1levels= outs1levels ||
repeat(s1levels[s1],1,ns2col);
 end;
 if (&nomarg=0) then  outs1levels=  repeat(outs1levels,1,ns3col+1);
 else                 outs1levels=  repeat(outs1levels,1,ns3col);
 if (&nomarg=0) then outs2levels= repeat(0 ||
s2levels,1,(ns1col+1)#(ns3col+1));
 else                outs2levels= repeat(s2levels,1,ns1col#ns3col);
end;

*---------------------------*
*  DMATRIX                  *
*---------------------------*;
if (&tabletype=0) then dmatrix= j(nobsread,1,1);
else if (&tabletype=1) then do;
 s1matrix= design(subgrp1);
 if (mvsubgrp1 & &levelmiss=0) then s1matrix= s1matrix[,1:ns1col];
 if (&nomarg=0) then dmatrix= j(nobsread,1,1) || s1matrix;
 else                dmatrix= s1matrix;
 if (mvsubgrp1 & &exclmissmarg) then dmatrix[loc(subgrp1=99e15),1]=0;
end;
else if (&tabletype=2) then do;
 s2matrix= design(subgrp2);
 if (mvsubgrp2 & &levelmiss=0) then s2matrix= s2matrix[,1:ns2col];
 if (&nomarg=0) then workmatrix= j(nobsread,1,1) || s2matrix;
 else                workmatrix= s2matrix;
 do s1=1 to ns1col;
  blockmatrix= workmatrix#(subgrp1=s1levels[s1]);
  if (mvsubgrp2 & &exclmissmarg) then blockmatrix[loc(subgrp2=99e15),1]=0;
  appendmatrix= appendmatrix || blockmatrix;
 end;
 if (&nomarg=0) then do;
  do c=2 to ns2col+1;
   if (mvsubgrp1 & &exclmissmarg) then workmatrix[loc(subgrp1=99e15),c]=0;
  end;
  if (mvsubgrp1 & &exclmissmarg) then workmatrix[loc(subgrp1=99e15),1]=0;
  if (mvsubgrp2 & &exclmissmarg) then workmatrix[loc(subgrp2=99e15),1]=0;
  dmatrix= workmatrix || appendmatrix;
 end;
 else dmatrix=appendmatrix;
end;
else if (&tabletype=3) then do;
 s2matrix= design(subgrp2);
 if (mvsubgrp2 & &levelmiss=0) then s2matrix= s2matrix[,1:ns2col];
 if (&nomarg=0) then workmatrix= j(nobsread,1,1) || s2matrix;
 else                workmatrix= s2matrix;
 if (&nomarg=0) then do s1=1 to ns1col;
  blockmatrix= workmatrix#(subgrp1=s1levels[s1]);
  do c=1 to ns2col+1;
   if (mvsubgrp1 & &exclmissmarg) then blockmatrix[loc(subgrp1=99e15),c]=0;
   if (mvsubgrp2 & &exclmissmarg) then blockmatrix[loc(subgrp2=99e15),c]=0;
   if (mvsubgrp3 & &exclmissmarg) then blockmatrix[loc(subgrp3=99e15),c]=0;
  end;
  appendmatrix= appendmatrix || blockmatrix;
 end;
 check=0;
 do s3=1 to ns3col;
  do s1=1 to ns1col;
   if (s3^=check & &nomarg=0) then do;
    blockmatrix=workmatrix#(subgrp3=s3levels[s3]);
    if (mvsubgrp1 & &exclmissmarg) then blockmatrix[loc(subgrp1=99e15),1]=0;
    if (mvsubgrp2 & &exclmissmarg) then blockmatrix[loc(subgrp2=99e15),1]=0;
    do c=2 to ns2col+1;
     if (mvsubgrp1 & &exclmissmarg) then blockmatrix[loc(subgrp1=99e15),c]=0;
    end;
    appendmatrix= appendmatrix || blockmatrix;
    check=s3;
    end;
   blockmatrix= workmatrix#(subgrp3=s3levels[s3])#(subgrp1=s1levels[s1]);
   if (mvsubgrp2 & &exclmissmarg) then blockmatrix[loc(subgrp2=99e15),1]=0;
   appendmatrix= appendmatrix || blockmatrix;
  end;
 end;
 if (&nomarg=0) then do;
  do c=2 to ns2col+1;
   if (mvsubgrp1 & &exclmissmarg) then workmatrix[loc(subgrp1=99e15),c]=0;
   if (mvsubgrp2 & &exclmissmarg) then workmatrix[loc(subgrp2=99e15),c]=0;
   if (mvsubgrp3 & &exclmissmarg) then workmatrix[loc(subgrp3=99e15),c]=0;
  end;
  if (mvsubgrp1 & &exclmissmarg) then workmatrix[loc(subgrp1=99e15),1]=0;
  if (mvsubgrp2 & &exclmissmarg) then workmatrix[loc(subgrp2=99e15),1]=0;
  if (mvsubgrp3 & &exclmissmarg) then workmatrix[loc(subgrp3=99e15),1]=0;
  dmatrix= workmatrix || appendmatrix;
 end;
 else dmatrix=appendmatrix;
end;
if (&tabletype=1) then free s1matrix subgrp1 mvsubgrp1;
else if (&tabletype=2) then
 free appendmatrix workmatrix blockmatrix s1matrix s2matrix
      subgrp1 subgrp2 mvsubgrp1 mvsubgrp2;
else if (&tabletype=3) then
 free appendmatrix workmatrix blockmatrix s1matrix s2matrix
      subgrp1 subgrp2 subgpr3 mvsubgrp1 mvsubgrp2 mvsubgrp3;

 *---if subpop, modify dmatrix---*;
if (&inoperator) then do;
 isub=j(nobsread,1,0);
 do e=1 to ninelements;
  inpart=(&subpopvar=inelements[e]);
  isub= isub + inpart;
 end;
end;
else isub=(&subpop2);
mvsubpop=(&subpopvar^=.);
isub=isub#mvsubpop;        free &subpopvar mvsubpop;
dmatrix=dmatrix#isub;
 *---if zerowts, modify matrices---*;
if (izerowt) then do;
  izerowt=1;
  keep=loc(w>0);
  delete=loc(w<=0);
  nobsused=ncol(keep);
  nobsdelete=ncol(delete);
   *---modify matrices---*;
  load strata cluster popsize var topvar botvar cvar;
  strata=strata[keep];
  stratadesign=design(strata);
  cluster=cluster[keep];
  popsize=popsize[keep];
  w=w[keep];
  isub=isub[keep];
  dmatrix=dmatrix[keep,];
  if (&ivar) then var=var[keep,];
  else do;
   topvar=topvar[keep,];
   botvar=botvar[keep,];
  end;
  if (&icontrast) then cvar=cvar[keep,];
  store strata cluster popsize var topvar botvar cvar;
  free strata cluster popsize var topvar botvar cvar;
end;
else do;
 izerowt=0;
 nobsused=nobsread;
 keep=0;
 nobsdelete=0;
end;

ndrec=dmatrix[+,];
ndcol=ncol(dmatrix);
store dmatrix; free dmatrix keep;
create _izerowt var{izerowt nobsread nobsused nobsdelete}; append;

*---------------------------*
*  MODULE: M_COMPUTE        *
*---------------------------*;
start m_compute;
*---cross counts---*;
nvar_ndcol=nvar#ndcol;

 *---linearized variables---*;
do v=1 to nvar;
  *---handle any missing---*;
 load var;
 tvar=var[,v]; free var;
 mv=(tvar^=.);
 if (ncol(loc(tvar=.))) then tvar[loc(tvar=.)]=0; load dmatrix;
 tdmatrix=dmatrix#mv;                       free dmatrix;
  *---count statistics---*;
 tnobs=   tdmatrix[+,];
 tnmiss=  ndrec-tnobs;
 tsumwgt= (w#tdmatrix)[+,];
 nobs= nobs || tnobs;
 nmiss= nmiss || tnmiss;
 sumwgt= sumwgt || tsumwgt;

  *---handle any zero columns in tdmatrix---*;
 if (ncol(loc(tdmatrix[+,]=0))) then do;
  _warn101=1;
  warn101= warn101 // _warn101;
  zerocol=loc(tdmatrix[+,]=0);
  nzerocol=ncol(zerocol);
  do z=1 to nzerocol;
   tdmatrix[,zerocol[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
  _warn101=0;
  warn101= warn101 // _warn101;
 end;

  *---variables---*;
 dtotal=   (tvar#w#tdmatrix)[+,];
 dw=       (w#tdmatrix)[+,];      store tdmatrix; free tdmatrix;
 dmean=    dtotal/dw;             load tdmatrix;
 tzdhij_t= (tvar#w)#tdmatrix;     store tzdhij_t; free tzdhij_t;
 tzdhij_m= (repeat(tvar,1,ndcol)-repeat(dmean,nobsused))
            #w#tdmatrix / repeat(dw,nobsused);
                                  store tzdhij_m; free tzdhij_m tdmatrix;

  *---if zero columns in tdmatrix, set columns to missing---*;
 if (_warn101) then dmean[zerocol]=.;

 if v=1 then do;       load tzdhij_t;
  zdhij_t= tzdhij_t;   store zdhij_t;  free zdhij_t tzdhij_t; load tzdhij_m;
  zdhij_m= tzdhij_m;   store zdhij_m;  free zdhij_m tzdhij_m;
  total=   dtotal;
  mean=    dmean;
 end;
 else do;                              load zdhij_t tzdhij_t;
 zdhij_t= zdhij_t || tzdhij_t;
                  store zdhij_t; free zdhij_t tzdhij_t; load zdhij_m tzdhij_m;
 zdhij_m= zdhij_m || tzdhij_m;  store zdhij_m tzdhij_m; free zdhij_m tzdhij_m;
 total=   total || dtotal;
 mean =   mean || dmean;
 end;
end;
                                 free mv tvar;
 *---start stratum calculations---*;
do h=1 to nstrata;
                                 load zdhij_t;
 zdhij_th=   zdhij_t[loc(stratadesign[,h]=1),];  free zdhij_t;
  *---determine number of psus---*;
 load cluster;
 cluster_h=  cluster[loc(stratadesign[,h]=1)];   free cluster;
 if (&icluster) then do;
   uniqclus=unique(cluster_h[,1]);
   nh=ncol(uniqclus);
   end;
 else nh=nrow(zdhij_th);
 npsu=npsu+nh;

  *---compute fpc factor---*;
 load strata popsize;
 popsizeh= max( popsize[loc(strata=uniqstr[h])] ); free strata popsize;
 if popsizeh=0        then fh=1;
 else if popsizeh=-1  then fh=0;
 else                      fh= nh/popsizeh;

  *---err404---*;
 if (nh=1 & popsizeh^=0 & popsizeh^=1) then do;
  ERRSTR=h;
  ERR404=1;
 end;
 else errstr=-99;
 if (ERR404) then goto skip;

 *---compute stratum variance---*;
 if fh=1 then do;
  vtotalh=0; vmeanh=0;
  vtotal= vtotal + vtotalh;
  vmean= vmean + vmeanh;
 end;
 else do;
 *******************
 * total           *
 *******************;
  *---compute psu sums---*;
 if (&icluster) then do;
  zhi_t=j(nh,nvar_ndcol,0);
  do i=1 to nh;
   zhi_t[i,]= (zdhij_th[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_t= t(zhi_t);
  end;
 else do;
  zhi_t=t(zdhij_th);
  end;                                    free zdhij_th;
 *---variance---*;
 zhbar_t= repeat( zhi_t[,+]/nh, 1,nh);
 zdiff_t= zhi_t - zhbar_t;                free zhi_t zhbar_t;
 s2h_t= (zdiff_t#zdiff_t)[,+] / (nh-1);   free zdiff_t;
 vtotalh= (1-fh)#nh#s2h_t;
  *---accumulate stratum variances---*;
 vtotal= vtotal + vtotalh;
 *******************
 * mean            *
 *******************;
 load zdhij_m;
 zdhij_mh=   zdhij_m[loc(stratadesign[,h]=1),];
  *---compute psu sums---*;
 if (&icluster) then do;
  zhi_m=j(nh,nvar_ndcol,0);
  do i=1 to nh;
   zhi_m[i,]= (zdhij_mh[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_m= t(zhi_m);
  end;
 else do;
  zhi_m=t(zdhij_mh);
  end;                                  free zdhij_mh;
 *---variance---*;
 zhbar_m= repeat( zhi_m[,+]/nh, 1,nh);
 zdiff_m= zhi_m - zhbar_m;              free zhi_m zhbar_m;
 s2h_m= (zdiff_m#zdiff_m)[,+] / (nh-1); free zdiff_m;
 vmeanh=  (1-fh)#nh#s2h_m;
  *---accumulate stratum variances---*;
 vmean=  vmean + vmeanh;
 end;

end;

*-----------------------------------*
* Statistics                        *
*-----------------------------------*;
 *---statistics---*;
noutrec=nrow(vmean);
nstatvar=nvar;
ddf=npsu-nstrata;
tddf= repeat(ddf,noutrec);
nobs=t(nobs);
nmiss=t(nmiss);
sumwgt=t(sumwgt);
total= t(total);
mean=t(mean);
if ncol(loc(vtotal<0)) then vtotal[loc(vtotal<0)]=.;
if ncol(loc(vmean<0))  then vmean[loc(vmean<0)]=.;
setotal= sqrt(vtotal);
semean= sqrt(vmean);

 *---handle any zero columns in tdmatrix---*;
do vd=1 to nvar_ndcol;
 if (mean[vd]=.) then do;
  semean[vd]=.; total[vd]=.; setotal[vd]=.;
 end;
end;

 *---other output variables---*;
if (&isubpop) then subpop= repeat("&subpop",nvar_ndcol);
tableindex= repeat(_tableindex,nvar_ndcol);
&outsubgrp1=repeat(t(outs1levels),nvar,1);
&outsubgrp2=repeat(t(outs2levels),nvar,1);
&outsubgrp3=repeat(t(outs3levels),nvar,1);
varlist=t({&var});
do v=1 to nvar;
  variable=   variable // repeat(varlist[v],ndcol);
  varindex= varindex // repeat(v,ndcol);
end;

*-----------------------------------*
* output to datasets                *
*-----------------------------------*;
 *---remaining data summary ---*;
wsum=w[+];
subpopobs=isub[+];
subpopwsum=(w#isub)[+];

 *---create output dataset---*;
create _out var{tddf tableindex varindex variable &outsubpop &grouplist
 nobs nmiss sumwgt mean semean total setotal};
append;

 *---create data summary dataset---*;
if (&isubpop) then subpop="&subpop";
create _summary var{nstatvar &outsubpop nobsused subpopobs nstrata npsu ddf
                    wsum subpopwsum};
append;

skip:
finish m_compute;

*---------------------------*
*  MODULE: MC_COMPUTE       *
*---------------------------*;
start mc_compute;
load cvar;
ncvar=ncol(cvar);              free cvar;
clablist=t({"&clab1" "&clab2" "&clab3" "&clab4" "&clab5"
           "&clab6" "&clab7" "&clab8" "&clab9" "&clab10"});
cvarlist=t({&cvar});
cvarbase= {z_z_z};

 *---iterate contrasts---*;
do cv=1 to ncvar;
if (cv=1) then cvarcompare=cvarlist[1];
else if (cv>1) then do;
 cvarbase=cvarcompare;
 cvarcompare=cvarlist[cv];
end;

if (cvarcompare^=cvarbase) then do;   load cvar;
  *---tcvar---*;
 tcvar=cvar[,cv];                     free cvar;
 mvc=(tcvar^=.);
  chktcvar=tcvar;
  if (ncol(loc(chktcvar=.))) then chktcvar[loc(chktcvar=.)]=99;
  chklevels=unique(chktcvar);
  if (min(chklevels)<=0) then ERR403=1;
  if (ERR403) then goto skip;
 tcvar2=tcvar[loc(tcvar>=1)];
 tcvarlevels= unique(tcvar2);         free chktcvar tcvar2;
 nccol= ncol(tcvarlevels);
 if (ncol(loc(tcvar=.))) then tcvar[loc(tcvar=.)]=-1;
  *---cross counts---*;
 nvar_ndcol_nccol=nvar#ndcol#nccol;
 nvar_ndcol=      nvar#ndcol;
 ndcol_nccol=     ndcol#nccol;
  *---cmatrix---*;
 cmatrix=design(tcvar);
 if (ncol(loc(tcvar=-1))) then cmatrix=cmatrix[,2:nccol+1];
                                      store cmatrix; free cmatrix tcvar;
end;

 *---variance and covariance multipliers---*;
if cv=1 then coef={&coef1}; else if cv=2 then coef={&coef2};
else if cv=3 then coef={&coef3}; else if cv=4 then coef={&coef4};
else if cv=5 then coef={&coef5}; else if cv=6 then coef={&coef6};
else if cv=7 then coef={&coef7}; else if cv=8 then coef={&coef8};
else if cv=9 then coef={&coef9}; else if cv=10 then coef={&coef10};
tcontrastcoef= shape(coef,1,nccol,0);
varmults= tcontrastcoef##2;
covcomb= comb(nccol,2);
covmults= j(covcomb,1,0);
count=0;
do i=1 to nccol;
 do j=2 to nccol;
   if (i<j) then do;
      count= count+1;
      covmults[count]= tcontrastcoef[i]#tcontrastcoef[j];
   end;
 end;
end;
covmults=t(covmults);

 *---itcvar, tndrec (for count statistics)---*;   load cmatrix;
cmatrix_coef= ((cmatrix#tcontrastcoef)^=0);
itcvar=  cmatrix_coef[,+];                   free cmatrix_coef; load dmatrix;
tndrec=  (dmatrix#itcvar)[+,];               store itcvar; free itcvar;

if (cvarcompare^=cvarbase) then do;
  *---dcmatrix---*;
 do d=1 to ndcol;
  dcmatrix= dcmatrix || dmatrix[,d]#cmatrix;
 end;                                         free dmatrix cmatrix;
  *---handle any zero columns in dcmatrix---*;
 if (ncol(loc(dcmatrix[+,]=0))) then do;
  _warn102=1;
  warn102=warn102 // _warn102;
  zerocol2=loc(dcmatrix[+,]=0);
  nzerocol2=ncol(zerocol2);
  do z=1 to nzerocol2;
   dcmatrix[,zerocol2[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
   _warn102=0;
   warn102= warn102 // _warn102;
 end;
 store dcmatrix; free dcmatrix;
end;
else free dmatrix cmatrix;

 *---linearized variables---*;
do v=1 to nvar;
  *---handle any missing---*;                load var;
 tvar=var[,v];                               free var;
 mv=(tvar^=.);                               load dmatrix;
 tdmatrix=dmatrix#mv#mvc;                    free dmatrix;
 if (ncol(loc(tvar=.))) then tvar[loc(tvar=.)]=0;  load dcmatrix;
 tdcmatrix=dcmatrix#mv;                      free dcmatrix;
 *---count statistics---*;                   load itcvar;
 tnobs=   (tdmatrix#itcvar)[+,];
 tnmiss=   tndrec-tnobs;
 tsumwgt= (w#tdmatrix#itcvar)[+,];           free tdmatrix itcvar;
 _nobs= _nobs || tnobs;
 _nmiss= _nmiss || tnmiss;
 _sumwgt= _sumwgt || tsumwgt;

  *---handle any zero columns in tdcmatrix---*;
 if (ncol(loc(tdcmatrix[+,]=0))) then do;
  _warn101=1;
  warn101= warn101 // _warn101;
  zerocol=loc(tdcmatrix[+,]=0);
  nzerocol=ncol(zerocol);
  do z=1 to nzerocol;
   tdcmatrix[,zerocol[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
  _warn101=0;
  warn101= warn101 // _warn101;
 end;

  *---variables---*;
 dctotal=   (tvar#w#tdcmatrix)[+,];
 dcw=       (w#tdcmatrix)[+,];         store tdcmatrix; free tdcmatrix;
 dcmean=    dctotal/dcw;               load tdcmatrix;
 tzdchij_t= (tvar#w)#tdcmatrix;        store tzdchij_t; free tzdchij_t;
 tzdchij_m= (repeat(tvar,1,ndcol_nccol)-repeat(dcmean,nobsused))
             #w#tdcmatrix / repeat(dcw,nobsused);
                                    store tzdchij_m; free tzdchij_m tdcmatrix;
  *---if zero columns in tdcmatrix, set columns to missing---*;
 if (_warn101) then dcmean[zerocol]=.;

  *---if zero columns in dcmatrix, set columns to missing---*;
 if (_warn102) then dcmean[zerocol2]=.;

 if v=1 then do;      load tzdchij_t;
  zdchij_t=tzdchij_t; store zdchij_t; free zdchij_t tzdchij_t; load tzdchij_m;
  zdchij_m=tzdchij_m; store zdchij_m; free zdchij_m tzdchij_m;
  ctotal=  dctotal;
  cmean=   dcmean;
 end;
 else do;    load zdchij_t tzdchij_t;
 zdchij_t= zdchij_t || tzdchij_t;
             store zdchij_t; free zdchij_t tzdchij_t; load zdchij_m tzdchij_m;
 zdchij_m= zdchij_m || tzdchij_m;
             store zdchij_m tzdchij_m; free zdchij_m tzdchij_m;
 ctotal=   ctotal || dctotal;
 cmean =   cmean || dcmean;
 end;
end;         free mv tvar;

 *---start stratum calculations---*;
do h=1 to nstrata;
                                                  load zdchij_t;
 zdchij_th=  zdchij_t[loc(stratadesign[,h]=1),];  free zdchij_t;
  *---determine number of psus---*;
 load cluster;
 cluster_h=   cluster[loc(stratadesign[,h]=1)];   free cluster;
 if (&icluster) then do;
  uniqclus=unique(cluster_h[,1]);
  nh=ncol(uniqclus);
  end;
 else nh=nrow(zdchij_th[,1]);
 if (cv=1) then npsu=npsu+nh;

  *---compute fpc factor---*;
 load strata popsize;
 popsizeh= max( popsize[loc(strata=uniqstr[h])] ); free strata popsize;
 if popsizeh=0        then fh=1;
 else if popsizeh=-1  then fh=0;
 else                      fh= nh/popsizeh;

  *---err404---*;
 if (nh=1 & popsizeh^=0 & popsizeh^=1) then do;
  ERRSTR=h;
  ERR404=1;
 end;
 else errstr=-99;
 if (ERR404) then goto skip;

 *---compute stratum variance---*;
 if fh=1 then do;
  vtotalh=0; vmeanh=0;
  vtotal= vtotal + vtotalh;
  vmean= vmean + vmeanh;
 end;
 else do;
 ********************
 * total            *
 ********************;
  *---compute psu sums---*;
 if (&icluster) then do;
  zhi_t=j(nh,nvar_ndcol_nccol,0);
  do i=1 to nh;
   zhi_t[i,]= (zdchij_th[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_t= t(zhi_t);
  end;
 else do;
  zhi_t=t(zdchij_th);
  end;                              free zdchij_th;

 *---cell variance---*;
 zhbar_t= repeat(zhi_t[,+]/nh,1,nh);
 zdiff_t= zhi_t - zhbar_t;          free zhi_t zhbar_t;
 s2h_t= (zdiff_t#zdiff_t)[,+] / (nh-1);
 vtotaldch= (1-fh)#nh#s2h_t;
 vtotaldch_tr= btran(t(vtotaldch),1,nccol);

 *---start computing covariances---*;

 *---transpose zdiff---*;
 cindex= repeat(t(tcvarlevels),nvar_ndcol);
 cdesign= design(cindex);
 zdiff_t_tr= j(nvar_ndcol,nh#nccol,0);
 do c=1 to nccol;
  zdiff_t_tr[1:nvar_ndcol, nh#(c-1)+1:nh#c]= zdiff_t[loc(cdesign[,c]=1),];
 end;
                                    free zdiff_t;
 *---covariance vector---*;
 c_th=  j(nvar_ndcol,covcomb,0);
 count=0;
 do i=1 to nccol;
  do j=2 to nccol;
   if (i<j) then do;
    count= count+1;
    c_th[,count]=
       (1-fh)#nh#
       (zdiff_t_tr[,nh#(i-1)+1:nh#i]#zdiff_t_tr[,nh#(j-1)+1:nh#j])[,+]/(nh-1);
    end;
  end;
 end;
                                    free zdiff_t_tr;
  *---stratum contrast variance---*;
 vtotalh=  (varmults#vtotaldch_tr)[,+] +  2#(covmults#c_th)[,+];
  *---accumulate stratum contrast variances---*;
 vtotal= vtotal + vtotalh;
 ********************
 * mean             *
 ********************;
 load zdchij_m;
 zdchij_mh=  zdchij_m[loc(stratadesign[,h]=1),];
                                   free zdchij_m;
   *---compute psu sums---*;
 if (&icluster) then do;
  zhi_m=j(nh,nvar_ndcol_nccol,0);
  do i=1 to nh;
   zhi_m[i,]= (zdchij_mh[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_m= t(zhi_m);
  end;
 else do;
  zhi_m=t(zdchij_mh);
  end;
                                    free zdchij_mh;
 *---cell variance---*;
 zhbar_m= repeat(zhi_m[,+]/nh,1,nh);
 zdiff_m= zhi_m - zhbar_m;          free zhi_m zhbar_m;
 s2h_m= (zdiff_m#zdiff_m)[,+] / (nh-1);
 vmeandch=  (1-fh)#nh#s2h_m;
 vmeandch_tr=  btran(t(vmeandch),1,nccol);

 *---start computing covariances---*;

 *---transpose zdiff---*;
 cindex= repeat(t(tcvarlevels),nvar_ndcol);
 cdesign= design(cindex);
 zdiff_m_tr= j(nvar_ndcol,nh#nccol,0);
 do c=1 to nccol;
  zdiff_m_tr[1:nvar_ndcol, nh#(c-1)+1:nh#c]= zdiff_m[loc(cdesign[,c]=1),];
 end;
                                    free zdiff_m;
 *---covariance vector---*;
 c_mh=  j(nvar_ndcol,covcomb,0);
 count=0;
 do i=1 to nccol;
  do j=2 to nccol;
   if (i<j) then do;
    count= count+1;
    c_mh[,count]=
       (1-fh)#nh#
       (zdiff_m_tr[,nh#(i-1)+1:nh#i]#zdiff_m_tr[,nh#(j-1)+1:nh#j])[,+]/(nh-1);
    end;
  end;
 end;
                                    free zdiff_m_tr;
  *---stratum contrast variance---*;
 vmeanh=   (varmults#vmeandch_tr)[,+]  +  2#(covmults#c_mh)[,+];
 *---accumulate stratum contrast variances---*;
 vmean=  vmean + vmeanh;
 end;

 end;

*-----------------------------------*
* Statistics                        *
*-----------------------------------*;
if (cv=1) then do;
 noutrec=nrow(vmean);
 nstatvar=nvar;
 ddf=npsu-nstrata;
 _tddf= repeat(ddf,noutrec);
end;

  *---cell statistics---*;
ctotal= btran(ctotal,1,nccol);
cmean=  btran(cmean,1,nccol);

 *---contrast statistics---*;
_total= (ctotal#tcontrastcoef)[,+];
_mean=  (cmean#tcontrastcoef)[,+];

 *---semean setotal---*;
if ncol(loc(vtotal<0)) then vtotal[loc(vtotal<0)]=.;
if ncol(loc(vmean<0))  then vmean[loc(vmean<0)]=.;
_setotal= sqrt(vtotal);
_semean= sqrt(vmean);

 *---handle any zero columns---*;
do vd=1 to nvar_ndcol;
 if (ncol(loc(cmean[vd,]=.))) then do;
  _mean[vd]=.; _semean[vd]=.; _total[vd]=.; _setotal[vd]=.;
 end;
end;

 *---other table output variables---*;
if (&isubpop) then _subpop= repeat("&subpop",nvar_ndcol);
if cv=1 then do;
 _tableindex= repeat(_tableindex,nvar_ndcol);
 _outsubgrp1=repeat(t(outs1levels),nvar,1);
 _outsubgrp2=repeat(t(outs2levels),nvar,1);
 _outsubgrp3=repeat(t(outs3levels),nvar,1);
 varlist=t({&var});
 do v=1 to nvar;
   _variable= _variable // repeat(varlist[v],ndcol);
   _varindex= _varindex // repeat(v,ndcol);
 end;
end;
_contrastindex=repeat(cv,nvar_ndcol);
_contrastlabel=repeat(clablist[cv],nvar_ndcol);

 *---append contrast-table blocks---*;
if (&isubpop) then subpop= subpop // _subpop;
nobs=       nobs     // t(_nobs);
nmiss=      nmiss    // t(_nmiss);
sumwgt=     sumwgt   // t(_sumwgt);
mean=       mean     // _mean;
semean=     semean   // _semean;
total=      total    // _total;
setotal=    setotal  // _setotal;
outsubgrp1=   outsubgrp1   // _outsubgrp1;
outsubgrp2=   outsubgrp2   // _outsubgrp2;
outsubgrp3=   outsubgrp3   // _outsubgrp3;
tableindex=    tableindex    // _tableindex;
varindex=      varindex      // _varindex;
variable=      variable      // _variable;
contrastindex= contrastindex // _contrastindex;
contrastlabel= contrastlabel // _contrastlabel;
tddf=          tddf // _tddf;

 *---contrast vector variables---*;
_contrastindex2= repeat(cv,nccol,1);
_contrastlabel2= repeat(clablist[cv],nccol,1);

 *--append contrast vectors---*;
contrastindex2= contrastindex2 // _contrastindex2;
contrastlabel2= contrastlabel2 // _contrastlabel2;
cvarlevels=     cvarlevels // t(tcvarlevels);
contrastcoef=   contrastcoef // t(tcontrastcoef);

 *---free matrices for next iteration---*;
free _nobs _nmiss _sumwgt ctotal cmean;
vtotal=0; vmean=0;

end;        * END OF CONTRAST ITERATION;

&outsubgrp1=outsubgrp1;
&outsubgrp2=outsubgrp2;
&outsubgrp3=outsubgrp3;

*-----------------------------------*
* output to datasets                *
*-----------------------------------*;
 *---remaining data summary ---*;
wsum=w[+];
subpopobs=isub[+];
subpopwsum=(w#isub)[+];

 *---create output dataset---*;
create _out var{tddf tableindex varindex variable &outsubpop &grouplist
 nobs nmiss sumwgt mean semean total setotal
 contrastindex contrastlabel};
append;

 *---create data summary dataset---*;
if (&isubpop) then subpop="&subpop";
create _summary var{ncvar nstatvar &outsubpop nobsused subpopobs nstrata npsu
                    ddf wsum subpopwsum};
append;

 *---create contrast vector dataset---*;
create _contrast var{contrastindex2 contrastlabel2 cvarlevels contrastcoef};
append;

skip:
finish mc_compute;

*---------------------------*
*  MODULE: R_COMPUTE        *
*---------------------------*;
start r_compute;
*---cross counts---*;
nratio_ndcol=nratio#ndcol;
 *---linearized variables---*;
do r=1 to nratio;
  *---handle any missing---*;
 load topvar botvar;
 ttopvar=topvar[,r]; free topvar;
 tbotvar=botvar[,r]; free botvar;
 mvt=(ttopvar^=.);
 mvb=(tbotvar^=.);
 mv=mvt#mvb;
 if (ncol(loc(ttopvar=.))) then ttopvar[loc(ttopvar=.)]=0;
 if (ncol(loc(tbotvar=.))) then tbotvar[loc(tbotvar=.)]=0;
 ttopvar=ttopvar#mv;
 tbotvar=tbotvar#mv;                  load dmatrix;
 tdmatrix=dmatrix#mv;                 free dmatrix;
  *---count statistics---*;
 tnobs=   tdmatrix[+,];
 tnmiss=  ndrec-tnobs;
 tsumwgt= (w#tdmatrix)[+,];
 nobs= nobs || tnobs;
 nmiss= nmiss || tnmiss;
 sumwgt= sumwgt || tsumwgt;

  *---handle any zero columns in tdmatrix---*;
 if (ncol(loc(tdmatrix[+,]=0))) then do;
  _warn101=1;
  warn101= warn101 // _warn101;
  zerocol=loc(tdmatrix[+,]=0);
  nzerocol=ncol(zerocol);
  do z=1 to nzerocol;
   tdmatrix[,zerocol[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
  _warn101=0;
  warn101= warn101 // _warn101;
 end;

 *---variables---*;
 if (topequal[r]) then ttopvar=(ttopvar=topvalue[r]);
 if (botequal[r]) then tbotvar=(tbotvar=botvalue[r]);
 dtop=    (ttopvar#w#tdmatrix)[+,];
 dbot=    (tbotvar#w#tdmatrix)[+,];          store tdmatrix; free tdmatrix;
  *---handle any zero columns in dbot---*;
 if (ncol(loc(dbot=0))) then do;
  _warn103=1;
  warn103=warn103 // _warn103;
  zerocol3=loc(dbot=0);
  dbot[zerocol3]=-1;
 end;
 else do;
  _warn103=0;
  warn103=warn103 // _warn103;
 end;

 dratio=   dtop/dbot;                  load tdmatrix;
 tzdhij_r= ( repeat(ttopvar,1,ndcol) -
           repeat(dratio,nobsused)#repeat(tbotvar,1,ndcol) )
           #w#tdmatrix / repeat(dbot,nobsused);
                                       store tzdhij_r; free tzdhij_r tdmatrix;

  *---if zero columns in tdmatrix, set columns to missing---*;
 if (_warn101) then dratio[zerocol]=.;

  *---if zero columns in dbot, set columns to missing---*;
 if (_warn103) then dratio[zerocol3]=.;

 if r=1 then do;
  ratio=     ratio    || dratio;     load tzdhij_r;
  zdhij_r=   zdhij_r  || tzdhij_r;   store zdhij_r; free zdhij_r tzdhij_r;
  numtotal=  numtotal || dtop;
  dentotal=  dentotal || dbot;
 end;
 else do;
  ratio=    ratio    || dratio;      load zdhij_r tzdhij_r;
  zdhij_r=  zdhij_r  || tzdhij_r;    store zdhij_r; free zdhij_r tzdhij_r;
  numtotal= numtotal || dtop;
  dentotal= dentotal || dbot;
 end;

end;                                 free mv mvt mvb ttopvar tbotvar;

 *---start stratum calculations---*;
do h=1 to nstrata;
                                                 load zdhij_r;
 zdhij_rh=   zdhij_r[loc(stratadesign[,h]=1),];  free zdhij_r;
  *---determine number of psus---*;
 load cluster;
 cluster_h=  cluster[loc(stratadesign[,h]=1)];   free cluster;
 if (&icluster) then do;
   uniqclus=unique(cluster_h[,1]);
   nh=ncol(uniqclus);
   end;
 else nh=nrow(zdhij_rh);
 npsu=npsu+nh;

  *---compute fpc factor---*;
 load strata popsize;
 popsizeh= max( popsize[loc(strata=uniqstr[h])] ); free strata popsize;
 if popsizeh=0        then fh=1;
 else if popsizeh=-1  then fh=0;
 else                      fh= nh/popsizeh;

  *---err404---*;
 if (nh=1 & popsizeh^=0 & popsizeh^=1) then do;
  ERRSTR=h;
  ERR404=1;
 end;
 else errstr=-99;
 if (ERR404) then goto skip;

 *---compute stratum variance---*;
 if fh=1 then do;
  vratioh=0;
  vratio= vratio + vratioh;
 end;
 else do;
  *---compute psu sums---*;
 if (&icluster) then do;
  zhi_r=j(nh,nratio_ndcol,0);
  do i=1 to nh;
   zhi_r[i,]= (zdhij_rh[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_r= t(zhi_r);
  end;
 else do;
  zhi_r= t(zdhij_rh);
  end;                                    free zdhij_rh;
  *---variance---*;
 zhbar_r= repeat( zhi_r[,+]/nh, 1,nh);
 zdiff_r= zhi_r - zhbar_r;                free zhi_r zhbar_r;
 s2h_r= (zdiff_r#zdiff_r)[,+] / (nh-1);   free zdiff_r;
 vratioh=  (1-fh)#nh#s2h_r;

  *---accumulate stratum variances---*;
 vratio=  vratio + vratioh;
 end;

end;

*-----------------------------------*
* Statistics                        *
*-----------------------------------*;
 *---statistics---*;
noutrec=nrow(vratio);
nstatvar=nratio;
ddf=npsu-nstrata;
tddf= repeat(ddf,noutrec);
nobs=t(nobs);
nmiss=t(nmiss);
sumwgt=t(sumwgt);
numtotal=t(numtotal);
dentotal=t(dentotal);
ratio=t(ratio);
if ncol(loc(vratio<0)) then vratio[loc(ratio<0)]=.;
seratio= sqrt(vratio);

 *---handle any zero columns in tdmatrix---*;
do rd=1 to nratio_ndcol;
 if (ratio[rd]=.) then do;
  numtotal[rd]=.; dentotal[rd]=.; ratio[rd]=.; seratio[rd]=.;
 end;
end;

 *---other output variables---*;
if (&isubpop) then subpop= repeat("&subpop",nratio_ndcol);
tableindex= repeat(_tableindex,nratio_ndcol);
&outsubgrp1=repeat(t(outs1levels),nratio,1);
&outsubgrp2=repeat(t(outs2levels),nratio,1);
&outsubgrp3=repeat(t(outs3levels),nratio,1);
ratiolist=  t({"&ratio1" "&ratio2" "&ratio3" "&ratio4" "&ratio5"
               "&ratio6" "&ratio7" "&ratio8" "&ratio9" "&ratio10"});
do r=1 to nratio;
 rationame= rationame // repeat(ratiolist[r],ndcol);
 ratioindex= ratioindex // repeat(r,ndcol);
end;

*-----------------------------------*
* output to datasets                *
*-----------------------------------*;
 *---remaining data summary ---*;
wsum=w[+];
subpopobs=isub[+];
subpopwsum=(w#isub)[+];

 *---create output dataset---*;
create _out var{tddf tableindex ratioindex rationame &outsubpop &grouplist
 nobs nmiss sumwgt numtotal dentotal ratio seratio};
append;

 *---create data summary dataset---*;
if (&isubpop) then subpop="&subpop";
create _summary var{nstatvar &outsubpop nobsused subpopobs nstrata npsu ddf
                    wsum subpopwsum};
append;

skip:
finish r_compute;

*---------------------------*
*  MODULE: RC_COMPUTE       *
*---------------------------*;
start rc_compute;
load cvar;
ncvar=ncol(cvar);                          free cvar;
clablist=t({"&clab1" "&clab2" "&clab3" "&clab4" "&clab5"
           "&clab6" "&clab7" "&clab8" "&clab9" "&clab10"});
cvarlist=t({&cvar});
cvarbase={z_z_z};

 *---iterate contrasts---*;
do cv=1 to ncvar;
if (cv=1) then cvarcompare=cvarlist[1];
else if (cv>1) then do;
 cvarbase=cvarcompare;
 cvarcompare=cvarlist[cv];
end;

if (cvarcompare^=cvarbase) then do;        load cvar;
  *---tcvar---*;
 tcvar=cvar[,cv];                          free cvar;
 mvc=(tcvar^=.);
  chktcvar=tcvar;
  if (ncol(loc(chktcvar=.))) then chktcvar[loc(chktcvar=.)]=99;
  chklevels=unique(chktcvar);
  if (min(chklevels)<=0) then ERR403=1;
  if (ERR403) then goto skip;
 tcvar2=tcvar[loc(tcvar>=1)];
 tcvarlevels= unique(tcvar2);              free chktcvar tcvar2;
 nccol= ncol(tcvarlevels);
 if (ncol(loc(tcvar=.))) then tcvar[loc(tcvar=.)]=-1;
  *---cross counts---*;
 nratio_ndcol_nccol=nratio#ndcol#nccol;
 nratio_ndcol=nratio#ndcol;
 ndcol_nccol= ndcol#nccol;
  *---cmatrix---*;
 cmatrix=design(tcvar);
 if (ncol(loc(tcvar=-1))) then cmatrix=cmatrix[,2:nccol+1];
                              store cmatrix; free cmatrix tcvar;
end;

 *---variance and covariance multipliers---*;
if cv=1 then coef={&coef1}; else if cv=2 then coef={&coef2};
else if cv=3 then coef={&coef3}; else if cv=4 then coef={&coef4};
else if cv=5 then coef={&coef5}; else if cv=6 then coef={&coef6};
else if cv=7 then coef={&coef7}; else if cv=8 then coef={&coef8};
else if cv=9 then coef={&coef9}; else if cv=10 then coef={&coef10};
tcontrastcoef= shape(coef,1,nccol,0);
varmults= tcontrastcoef##2;
covcomb= comb(nccol,2);
covmults= j(covcomb,1,0);
count=0;
do i=1 to nccol;
 do j=2 to nccol;
   if (i<j) then do;
      count= count+1;
      covmults[count]= tcontrastcoef[i]#tcontrastcoef[j];
   end;
 end;
end;
covmults=t(covmults);

 *---itcvar, tndrec (for count statistics)---*;  load cmatrix;
cmatrix_coef= ((cmatrix#tcontrastcoef)^=0);
itcvar=  cmatrix_coef[,+];                   free cmatrix_coef; load dmatrix;
tndrec=  (dmatrix#itcvar)[+,];               store itcvar; free itcvar;

if (cvarcompare^=cvarbase) then do;
  *---dcmatrix---*;
 do d=1 to ndcol;
  dcmatrix= dcmatrix || dmatrix[,d]#cmatrix;
 end;                                         free dmatrix cmatrix;
  *---handle any zero columns in dcmatrix---*;
 if (ncol(loc(dcmatrix[+,]=0))) then do;
  _warn102=1;
  warn102= warn102 // _warn102;
  zerocol2=loc(dcmatrix[+,]=0);
  nzerocol2=ncol(zerocol2);
  do z=1 to nzerocol2;
   dcmatrix[,zerocol2[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
  _warn102=0;
  warn102=warn102 // _warn102;
 end;
 store dcmatrix; free dcmatrix;
end;
else free dmatrix cmatrix;

 *---linearized variables---*;
do r=1 to nratio;
  *---handle any missing---*;   load topvar botvar;
 ttopvar=topvar[,r];            free topvar;
 tbotvar=botvar[,r];            free botvar;
 mvt=(ttopvar^=.);
 mvb=(tbotvar^=.);
 mv=mvt#mvb;                                          load dmatrix;
 tdmatrix=dmatrix#mv#mvc;                             free dmatrix;
 if (ncol(loc(ttopvar=.))) then ttopvar[loc(ttopvar=.)]=0;
 if (ncol(loc(tbotvar=.))) then tbotvar[loc(tbotvar=.)]=0;
 ttopvar=ttopvar#mv;
 tbotvar=tbotvar#mv;                                  load dcmatrix;
 tdcmatrix=dcmatrix#mv;                               free dcmatrix;
 *---count statistics---*;                            load itcvar;
 tnobs=   (tdmatrix#itcvar)[+,];
 tnmiss=   tndrec-tnobs;
 tsumwgt= (w#tdmatrix#itcvar)[+,];                    free tdmatrix itcvar;
 _nobs= _nobs || tnobs;
 _nmiss= _nmiss || tnmiss;
 _sumwgt= _sumwgt || tsumwgt;
  *---handle any zero columns in tdcmatrix---*;
 if (ncol(loc(tdcmatrix[+,]=0))) then do;
  _warn101=1;
  warn101=warn101 // _warn101;
  zerocol=loc(tdcmatrix[+,]=0);
  nzerocol=ncol(zerocol);
  do z=1 to nzerocol;
   tdcmatrix[,zerocol[z]]=j(nobsused,1,-1);
  end;
 end;
 else do;
  _warn101=0;
  warn101=warn101 // _warn101;
 end;

  *---variables---*;
 if (topequal[r]) then ttopvar=(ttopvar=topvalue[r]);
 if (botequal[r]) then tbotvar=(tbotvar=botvalue[r]);
 dctop=    (ttopvar#w#tdcmatrix)[+,];
 dcbot=    (tbotvar#w#tdcmatrix)[+,];         store tdcmatrix; free tdcmatrix;
  *---handle any zero columns in dbot---*;
 if (ncol(loc(dcbot=0))) then do;
  _warn103=1;
  warn103=warn103 // _warn103;
  zerocol3=loc(dcbot=0);
  dcbot[zerocol3]=-1;
 end;
 else do;
  _warn103=0;
  warn103=warn103 // _warn103;
 end;

 dcratio=  dctop/dcbot;                       load tdcmatrix;
 tzdchij_r= ( repeat(ttopvar,1,ndcol_nccol) -
            repeat(dcratio,nobsused)#repeat(tbotvar,1,ndcol_nccol) )
            #w#tdcmatrix / repeat(dcbot,nobsused);
                                    store tzdchij_r; free tzdchij_r tdcmatrix;

  *---if zero columns in tdcmatrix, set columns to missing---*;
 if (_warn101) then dcratio[zerocol]=.;

 *---if zero columns in dcmatrix, set columns to missing---*;
 if (_warn102) then; dcratio[zerocol2]=.;

  *---if zero columns in dbot, set columns to missing---*;
 if (_warn103) then dcratio[zerocol3]=.;

 if r=1 then do;
  cratio=    cratio    || dcratio;    load tzdchij_r;
  zdchij_r=  zdchij_r  || tzdchij_r;  store zdchij_r; free zdchij_r tzdchij_r;
  cnumtotal= cnumtotal || dctop;
  cdentotal= cdentotal || dcbot;

 end;
 else do;
  cratio=    cratio    || dcratio;    load zdchij_r tzdchij_r;
  zdchij_r=  zdchij_r  || tzdchij_r;  store zdchij_r; free zdchij_r tzdchij_r;
  cnumtotal= cnumtotal || dctop;
  cdentotal= cdentotal || dcbot;
 end;
end;                                  free mv mvt mvb ttopvar tbotvar;

 *---start stratum calculations---*;
do h=1 to nstrata;
                                                 load zdchij_r;
 zdchij_rh=  zdchij_r[loc(stratadesign[,h]=1),]; free zdchij_r;

 *---determine number of psus---*;
 load cluster;
 cluster_h=   cluster[loc(stratadesign[,h]=1)];  free cluster;
 if (&icluster) then do;
  uniqclus=unique(cluster_h[,1]);
  nh=ncol(uniqclus);
  end;
 else nh=nrow(zdchij_rh[,1]);
 if (cv=1) then npsu=npsu+nh;

  *---compute fpc factor---*;
 load strata popsize;
 popsizeh= max( popsize[loc(strata=uniqstr[h])] ); free strata popsize;
 if popsizeh=0        then fh=1;
 else if popsizeh=-1  then fh=0;
 else                      fh= nh/popsizeh;

  *---err404---*;
 if (nh=1 & popsizeh^=0 & popsizeh^=1) then do;
  ERRSTR=h;
  ERR404=1;
 end;
 else errstr=-99;
 if (ERR404) then goto skip;

 *---compute stratum variance---*;
 if fh=1 then do;
  vratioh=0;
  vratio= vratio + vratioh;
 end;
 else do;
  *---compute psu sums---*;
 if (&icluster) then do;
  zhi_r=j(nh,nratio_ndcol_nccol,0);
  do i=1 to nh;
   zhi_r[i,]= (zdchij_rh[loc(cluster_h=uniqclus[i]),]) [+,];
  end;
  zhi_r= t(zhi_r);
  end;
 else do;
  zhi_r=t(zdchij_rh);
  end;                                    free zdchij_rh;
 *---cell variance---*;
 zhbar_r= repeat( zhi_r[,+]/nh, 1,nh);
 zdiff_r= zhi_r - zhbar_r;                free zhi_r zhbar_r;
 s2h_r= (zdiff_r#zdiff_r)[,+] / (nh-1);
 vratiodch=  (1-fh)#nh#s2h_r;
 vratiodch_tr=  btran(t(vratiodch),1,nccol);

 *---start computing covariances---*;

 *---transpose zdiff---*;
 cindex= repeat(t(tcvarlevels),nratio_ndcol);
 cdesign= design(cindex);
 zdiff_r_tr= j(nratio_ndcol,nh#nccol,0);
 do c=1 to nccol;
  zdiff_r_tr[1:nratio_ndcol, nh#(c-1)+1:nh#c]= zdiff_r[loc(cdesign[,c]=1),];
 end;
                                          free zdiff_r;
 *---covariance vector---*;
 c_rh=  j(nratio_ndcol,covcomb,0);
 count=0;
 do i=1 to nccol;
  do j=2 to nccol;
   if (i<j) then do;
    count= count+1;
     c_rh[,count]=
      (1-fh)#nh#
      (zdiff_r_tr[,nh#(i-1)+1:nh#i]#zdiff_r_tr[,nh#(j-1)+1:nh#j])[,+] /(nh-1);
    end;
  end;
 end;                                     free zdiff_r_tr;
  *---stratum contrast variance---*;
 vratioh=   (varmults#vratiodch_tr)[,+]  +  2#(covmults#c_rh)[,+];
  *---accumulate stratum contrast variances---*;
 vratio=  vratio + vratioh;
 end;

end;

*-----------------------------------*
* Statistics                        *
*-----------------------------------*;
if (cv=1) then do;
 noutrec=nrow(vratio);
 nstatvar=nratio;
 ddf=npsu-nstrata;
 _tddf= repeat(ddf,noutrec);
end;

  *---cell statistics---*;
cnumtotal=  btran(cnumtotal,1,nccol);
cdentotal=  btran(cdentotal,1,nccol);
cratio=  btran(cratio,1,nccol);

 *---contrast statistics---*;
_ratio=     (cratio#tcontrastcoef)[,+];
_numtotal=  (cnumtotal#tcontrastcoef)[,+];
_dentotal=  (cdentotal#tcontrastcoef)[,+];

 *---seratio---*;
if ncol(loc(vratio<0))  then vratio[loc(vratio<0)]=.;
_seratio= sqrt(vratio);

 *---handle any zero columns---*;
do rd=1 to nratio_ndcol;
 if (ncol(loc(cratio[rd,]=.))) then do;
  _numtotal[rd]=.; _dentotal[rd]=.; _ratio[rd]=.; _seratio[rd]=.;
 end;
end;

 *---other output variables---*;
if (&isubpop) then _subpop= repeat("&subpop",nratio_ndcol);
if cv=1 then do;
 _tableindex= repeat(_tableindex,nratio_ndcol);
 _outsubgrp1=repeat(t(outs1levels),nratio,1);
 _outsubgrp2=repeat(t(outs2levels),nratio,1);
 _outsubgrp3=repeat(t(outs3levels),nratio,1);
 ratiolist=  t({"&ratio1" "&ratio2" "&ratio3" "&ratio4" "&ratio5"
                "&ratio6" "&ratio7" "&ratio8" "&ratio9" "&ratio10"});
 do r=1 to nratio;
  _rationame=  _rationame // repeat(ratiolist[r],ndcol);
  _ratioindex= _ratioindex // repeat(r,ndcol);
 end;
end;
_contrastindex=repeat(cv,nratio_ndcol);
_contrastlabel=repeat(clablist[cv],nratio_ndcol);

 *---append contrast-table blocks---*;
if (&isubpop) then subpop= subpop // _subpop;
nobs=       nobs       //  t(_nobs);
nmiss=      nmiss      //  t(_nmiss);
sumwgt=     sumwgt     //  t(_sumwgt);
numtotal=   numtotal   //  _numtotal;
dentotal=   dentotal   //  _dentotal;
ratio=      ratio      //  _ratio;
seratio=    seratio    //  _seratio;
outsubgrp1= outsubgrp1 //  _outsubgrp1;
outsubgrp2= outsubgrp2 //  _outsubgrp2;
outsubgrp3= outsubgrp3 //  _outsubgrp3;
tableindex= tableindex //  _tableindex;
ratioindex= ratioindex //  _ratioindex;
rationame=  rationame  //  _rationame;
contrastindex= contrastindex // _contrastindex;
contrastlabel= contrastlabel // _contrastlabel;
tddf=          tddf // _tddf;

 *---contrast vector variables---*;
_contrastindex= repeat(cv,nccol,1);
_contrastlabel= repeat(clablist[cv],nccol,1);

 *--append contrast vectors---*;
contrastindex2= contrastindex2 // _contrastindex;
contrastlabel2= contrastlabel2 // _contrastlabel;
cvarlevels=     cvarlevels     // t(tcvarlevels);
contrastcoef=   contrastcoef   // t(tcontrastcoef);

 *---free matrices for next iteration---*;
free _nobs _nmiss _sumwgt cnumtotal cdentotal cratio;
vratio=0;

end;        * END OF CONTRAST ITERATION;

&outsubgrp1=outsubgrp1;
&outsubgrp2=outsubgrp2;
&outsubgrp3=outsubgrp3;

*-----------------------------------*
* output to datasets                *
*-----------------------------------*;
 *---remaining data summary ---*;
wsum=w[+];
subpopobs=isub[+];
subpopwsum=(w#isub)[+];

 *---create output dataset---*;
create _out var{tddf tableindex ratioindex rationame &outsubpop &grouplist
 nobs nmiss sumwgt numtotal dentotal ratio seratio
 contrastindex contrastlabel};
append;

 *---create data summary dataset---*;
if (&isubpop) then subpop="&subpop";
create _summary var{ncvar nstatvar &outsubpop nobsused subpopobs nstrata npsu
                    ddf wsum subpopwsum};
append;

 *---create contrast vector dataset---*;
create _contrast var{contrastindex2 contrastlabel2 cvarlevels contrastcoef};
append;

skip:
finish rc_compute;

*---------------------------*
*  RUN MODULES              *
*---------------------------*;
if (ERR401=0 & ERR402=0) then do;
 if (&ivar=1 & &icontrast=0)      then run m_compute;
 else if (&ivar=1 & &icontrast=1) then run mc_compute;
 else if (&ivar=0 & &icontrast=0) then run r_compute;
 else if (&ivar=0 & &icontrast=1) then run rc_compute;

  *---print any necessary warning messages---*;
 if (warn101[+]) then do;
  file log; put;
  put 'WARNING: At least one subgroup cell is not represented due to missing';
  put '         values or a subpop request.  Some statistics were set to';
  put '         missing.';
  put;
 end;
 if (warn102[+]) then do;
  file log; put;
  put 'WARNING: At least one contrast level in a subgroup cell is not';
  put '         represented.  Some statistics were set to missing.';
  put;
 end;
 if (warn103[+]) then do;
  file log; put;
  put 'WARNING: At least one subgroup cell contains a zero denominator in';
  put '         the ratio.  Some statistics were set to missing.';
  put;
 end;
end;

skip:
  *---output error indicators to dataset---*;
 create _err var{ERR401 ERR402 ERR403 ERR404 ERR405 ERR406 ERR407 ERRSTR};
 append;

free / ;
finish compute;
run compute;
quit;

%mend compute;

*--------------------------------------------------------*
*                                                        *
*    OUTPUT                                              *
*                                                        *
*--------------------------------------------------------*;
%macro output;
 *---process alpha option---*;
%if (%length(&alpha)) %then %let cipct=%sysevalf(100*(1-&alpha))%;
%else %let cipct=95%;

 *---label for contrastlabel---*;
%if (&icontrast) %then %do;
 %let labelcontrastlabel=%str(contrastlabel='Contrast Label');
 %let labelcontrastindex=%str(contrastindex='Contrast Index');
%end;
%else %do;
 %let labelcontrastlabel=;
 %let labelcontrastindex=;
%end;

 *---summary output---*;
data _null_;
set _izerowt;
call symput('izerowt',izerowt);
data _null_;

%if (&izerowt) %then %do;
 %let outnobsused=nobsused;
 %let outnobsdelete=nobsdelete;
 %let nobsusedlabel=nobsused='Number of Observations Used';
 %let nobsdeletelabel=nobsdelete='Number of Observations Deleted';
%end;
%else %do;
 %let outnobsused=;
 %let outnobsdelete=;
 %let nobsusedlabel=;
 %let nobsdeletelabel=;
%end;

%if (&tableindex=1) %then %do;
 data &outsummary(keep=&outsubpop nobsread &outnobsused &outnobsdelete
                       &outsubpopobs nstrata npsu ddf wsum &outsubpopwsum);
 merge _summary _izerowt;
 label
 &nobsusedlabel
 &nobsdeletelabel
 &subpoplabel
 &subpopobslabel
 nobsread='Number of Observations Read'
 nstrata='Number of Strata'
 npsu='Number of PSUs'
 ddf='Denominator Degrees of Freedom'
 wsum='Weight Sum'
 &subpopwsumlabel
 ;

 *---contrast output---*;
 %if  (&icontrast) %then %do;
  data &outcontrast(keep=contrastlabel contrastvar values coefficients);
  length CONTRASTVAR $ &maxlength;
  set _contrast;
  by contrastindex2;
  if contrastindex2=1 then CONTRASTVAR="&cvar1";
  else if contrastindex2=2 then CONTRASTVAR="&cvar2";
  else if contrastindex2=3 then CONTRASTVAR="&cvar3";
  else if contrastindex2=4 then CONTRASTVAR="&cvar4";
  else if contrastindex2=5 then CONTRASTVAR="&cvar5";
  else if contrastindex2=6 then CONTRASTVAR="&cvar6";
  else if contrastindex2=7 then CONTRASTVAR="&cvar7";
  else if contrastindex2=8 then CONTRASTVAR="&cvar8";
  else if contrastindex2=9 then CONTRASTVAR="&cvar9";
  else if contrastindex2=10 then CONTRASTVAR="&cvar10";
  CONTRASTLABEL=contrastlabel2;
  VALUES=cvarlevels;
  COEFFICIENTS=contrastcoef;
  if (not first.contrastindex2) then do;
   contrastvar=' ';
   contrastlabel=' ';
  end;
  label
  contrastlabel='Contrast Label'
  contrastvar='Contrast Variable'
  values='Contrast Variable Values'
  coefficients='Contrast Coefficients'
  ;
 %end;

%end;

 *---table output---*;
%if (&ivar) %then %do;
data _toutput(keep=tableindex &indextype &outsubpop &statvar &grouplist
                  &finalstats_both &outcontrastlabel &outcontrastindex);
 set _out;
  *---create variables for printing ---*;
 nmiss_p=nmiss;
 sumwgt_p=sumwgt;
 mean_p=mean;
  *---clmean---*;
 if (mean=.) then do;
  LCLMEAN=.;
  UCLMEAN=.;
 end;
 else if (tddf=0) then do;
  LCLMEAN=mean;
  UCLMEAN=mean;
 end;
 else do;
  LCLMEAN= mean - semean*tinv(1-&alpha/2,tddf);
  UCLMEAN= mean + semean*tinv(1-&alpha/2,tddf);
 end;
  *---cvmean---*;
 if (mean in (. 0)) then CVMEAN=.;
 else                    CVMEAN=semean/mean;
  *---tmean,pmean---*;
 if (semean in (. 0)) then do;
  TMEAN=.;
  PMEAN=.;
 end;
 else do;
  TMEAN=  mean/semean;
  PMEAN=  sdf('f',tmean**2,1,tddf);
 end;
  *---cltotal---*;
 if (total=.) then do;
  LCLTOTAL=.;
  UCLTOTAL=.;
 end;
 else if (tddf=0) then do;
  LCLTOTAL=total;
  UCLTOTAL=total;
 end;
 else do;
  LCLTOTAL= total - setotal*tinv(1-&alpha/2,tddf);
  UCLTOTAL= total + setotal*tinv(1-&alpha/2,tddf);
 end;
  *---cvtotal---*;
 if (total in (. 0)) then CVTOTAL=.;
 else                     CVTOTAL=setotal/total;
  *---labels---*;
 label
 tableindex='Table Index'
 varindex='Variable Index'
 variable='Variable'
 nobs='N'
 nmiss='N Miss'
 nmiss_p='N Miss'
 sumwgt='Sum of Weights'
 sumwgt_p='Sum of Weights'
 mean='Mean'
 mean_p='Mean'
 semean='Std Error  of Mean'
 lclmean="Lower" &cipct"  CL for Mean"
 uclmean="Upper" &cipct"  CL for Mean"
 cvmean='CV for   Mean'
 tmean='t Value Mean=0'
 pmean='Pr > |t| Mean=0'
 total='Total'
 setotal='Std Error  of Total'
 lcltotal="Lower" &cipct"  CL for Total"
 ucltotal="Upper" &cipct"  CL for Total"
 cvtotal='CV for  Total'
 &subpoplabel
 &labelcontrastindex
 &labelcontrastlabel
 ;
%end;
%else %do;
data _toutput(keep=tableindex &indextype &outsubpop &statvar &grouplist
                  &finalstats_both &outcontrastlabel &outcontrastindex);
 set _out;
  *---create variables for printing ---*;
 nmiss_p=nmiss;
 sumwgt_p=sumwgt;
  *---clratio---*;
 if (ratio=.) then do;
  LCLRATIO=.;
  UCLRATIO=.;
  end;
 else do;
  LCLRATIO= ratio - seratio*tinv(1-&alpha/2,tddf);
  UCLRATIO= ratio + seratio*tinv(1-&alpha/2,tddf);
  end;
  *---cvratio---*;
 if (ratio in (. 0)) then CVRATIO=.;
 else                     CVRATIO=seratio/ratio;
  *---tratio,pratio---*;
 if (seratio in (. 0)) then do;
  TRATIO=.;
  PRATIO=.;
 end;
 else do;
  TRATIO=  ratio/seratio;
  PRATIO=  sdf('f',tratio**2,1,tddf);
 end;
  *---labels---*;
 label
 tableindex='Table Index'
 ratioindex='Ratio Index'
 rationame='Ratio Name'
 nobs='N'
 nmiss='N Miss'
 nmiss_p='N Miss'
 sumwgt='Sum of Weights'
 sumwgt_p='Sum of Weights'
 numtotal='Numerator Total'
 dentotal='Denominator Total'
 ratio='Ratio'
 seratio='Std Error  of Ratio'
 lclratio="Lower" &cipct"  CL for Ratio"
 uclratio="Upper" &cipct"  CL for Ratio"
 cvratio='CV for  Ratio'
 tratio='t Value Ratio=0'
 pratio='Pr > |t| Ratio=0'
 &subpoplabel
 &labelcontrastindex
 &labelcontrastlabel
 ;
%end;
 *---append table output datasets---*;
%if (&icontrast and &ivar) %then %do;
 proc sort data=_toutput; by varindex &grouplist contrastindex; run;
%end;
%if (&icontrast and &ivar=0) %then %do;
 proc sort data=_toutput; by ratioindex &grouplist contrastindex; run;
%end;
%if (&tableindex=1) %then %do;
 data &outtable;
 set _toutput(keep=tableindex &indextype &outsubpop &statvar &grouplist
                  &finalstats_out &outcontrastlabel &outcontrastindex);
%end;
%else %do;
 data &outtable;
 set &outtable _toutput(keep=tableindex &indextype &outsubpop &statvar
              &grouplist &finalstats_out &outcontrastlabel &outcontrastindex);
%end;
%mend output;

*--------------------------------------------------------*
*                                                        *
*    PRINT                                               *
*                                                        *
*--------------------------------------------------------*;
%macro print;
*------------------*
* tables: prep     *
*------------------*;
 *---nstatvar---*;
data _null_;
set _summary;
call symput('nstatvar',nstatvar);
data _null_;
 *---proc format---*;
proc format;
value overall
0='Total'
;
value &s1fmt
&s1fmtcodes
0='Total'
99e15='Missing'
;
value &s2fmt
&s2fmtcodes
0='Total'
99e15='Missing'
;
value &s3fmt
&s3fmtcodes
0='Total'
99e15='Missing'
;
run;
%if (&tabletype=0) %then %do;
 %let printformat=_overall_ overall.;
 %let label=label _overall_='Overall';
 %let box=By: Overall;
%end;
%else %if (&tabletype=1) %then %do;
 %let printformat=&subgrp1 &s1fmt..;
 %let label=;
 %let box=By: &subgrp1;
%end;
%else %if (&tabletype=2) %then %do;
 %let printformat=&subgrp1 &s1fmt.. &subgrp2 &s2fmt..;
 %let label=;
 %let box=By: &subgrp1, &subgrp2;
%end;
%else %if (&tabletype=3) %then %do;
 %let printformat=&subgrp3 &s3fmt.. &subgrp1 &s1fmt.. &subgrp2 &s2fmt..;
 %let label=;
 %let box=By: &subgrp1, &subgrp2;
%end;

 *---if contrast, add to box---*;
%if (&icontrast) %then %let box=&box, Contrast;

%if (&style=BOX) %then %do;
*------------------*
* summary          *
*------------------*;
%if (&tableindex=1) %then %do;
 *---macro variables---*;
data _null_;
set _summary;
call symput('nobsused',nobsused);
call symput('subpopobs',subpopobs);
call symput('nstrata',nstrata);
call symput('npsu',npsu);
call symput('ddf',ddf);
call symput('wsum',round(wsum));
call symput('subpopwsum',round(subpopwsum));
data _null_;

data _temp(keep=noteclupop);
set _izerowt;
noteclupop=&iclupop;
call symput('izerowt',izerowt);
call symput('nobsread',nobsread);
call symput('nobsused',nobsused);
call symput('nobsdelete',nobsdelete);
run;
%if (&isubpop) %then %let subtitle=Subpopulation:   &subpop;
%if (&izerowt) %then
 %let nobsrow= put "Number of Observations Read:    &nobsread"    /
                   "Number of Observations Used:    &nobsused"
                      @51 "Weight Sum: &wsum"                     /
                   "Number of Observations Deleted: &nobsdelete"  /
                   " (Nonpositive Weights)";
%else %let nobsrow= put "Number of Observations Read:    &nobsread"
                   @51 "Weight Sum: &wsum";

%if (&isubpop=0) %then %do;
 data _null_;
 set _temp;
 file print;
 put "&usertitle";
 put;
 put "The SURVEYMEANS SUBGROUP Macro";
 put;
 put;
 put "            Data Summary";
 put;
 &nobsrow;
 put;
 put "Number of Strata:               &nstrata";
 put "Number of PSUs:                 &npsu";
 put "Denominator Degrees of Freedom: &ddf";
 if (noteclupop) then do;
 put;
 put"NOTE:  Finite population correction to sample PSUs (&cluster)";
 put"assumes inclusion of all listing units within sample PSUs.";
 end;
 run;
 %end;
%else %do;
 data _null_;
 set _temp;
 file print;
 put "&usertitle";
 put;
 put "The SURVEYMEANS SUBGROUP Macro";
 put;
 put;
 put "            Data Summary";
 put;
 put "&subtitle";
 put;
 &nobsrow;
 put "Subpopulation Observations:     &subpopobs"
         @51 "Weight Sum: &subpopwsum";
 put;
 put "Number of Strata:               &nstrata";
 put "Number of PSUs:                 &npsu";
 put "Denominator Degrees of Freedom: &ddf";
 if (noteclupop) then do;
 put;
 put"NOTE:  Finite population correction to sample PSUs (&cluster)";
 put"assumes inclusion of all listing units within sample PSUs.";
 end;
 run;
 %end;
*------------------*
* contrast         *
*------------------*;
%if (&icontrast) %then %do;
 proc print data=&outcontrast noobs label;
 var contrastlabel contrastvar values Coefficients;
 title'            Contrast Vector(s)';
 label
 contrastlabel='Label'
 contrastvar='Variable'
 values='Values'
 coefficients='Coefficients'
 ;
 run;
%end;
%end;

*------------------*
* tables           *
*------------------*;
%do t=1 %to &nstatvar;
 %if (&t=1) %then %do;
    %let tablerequest=%sysfunc(tranwrd(&tablerequest,#,*));
    %if (not &icontrast) %then
          %let tabletitle=title"            Table:  &tablerequest";
    %else %let tabletitle=title"            Table:  &tablerequest *Contrast";
    %end;
 %else %let tabletitle=title;
 %if (&ivar) %then %do;
   %let analvar=%scan(&var,&t);
   %let statvartitle=title4"Analysis Variable:  &analvar";
   %end;
 %else %do;
   proc sort data=_toutput out=_fname nodupkey; by ratioindex; run;
   data _null_;
   set _fname;
   if (_n_=&t) then call symput('analratio',rationame);
   run;
   %let statvartitle=title4"Ratio:              &analratio";
   %end;
 *---proc tabulate---*;
 %if not (&icontrast) %then %do;
  proc tabulate data=_toutput formchar='|----|+|---' noseps;
  class &grouplist;
  var &finalstats_print;
  table &printtablerequest, (&finalstats_print2)*max=''/
    box="&box" rts=&colwidth indent=4 nocontinued condense;
  format &printformat;
  &label;
  &tabletitle;
  &statvartitle;
  where &indextype=&t;
  run;
  %end;
  %else %if (&icontrast) %then %do;
  proc tabulate data=_toutput formchar='|----|+|---' noseps;
  class &grouplist contrastlabel;
  var &finalstats_print;
  table &printtablerequest*contrastlabel, (&finalstats_print2)*max=''/
    box="&box" rts=&colwidth indent=4 nocontinued condense;
  format &printformat;
  &label;
  &tabletitle;
  &statvartitle;
  where &indextype=&t;
  run;
  %end;
 %end;
%end;

%else %if (&style=WORDREADY) %then %do;
*------------------*
* summary          *
*------------------*;
%if (&tableindex=1) %then %do;

data _temp(keep=nobsdelete noteclupop);
set _izerowt;
noteclupop=&iclupop;
call symput('izerowt',izerowt);
run;
%if (&isubpop) %then %let subtitle=Subpopulation:  &subpop;
%if (&izerowt) %then
  %let nobsrow= put 'Number of Observations Read:#' nobsread            /
          'Number of Observations Used:#' nobsused '#Weight Sum:#' wsum /
          'Number of Observations Deleted:#' nobsdelete                 /
          '(Nonpositive Weights)';
%else %let nobsrow= put 'Number of Observations Read:#' nobsread
                        '#Weight Sum:#' wsum;

%if (&isubpop=0) %then %do;
data _null_;
merge &outsummary _temp;
file print;
 put "&usertitle";
 put;
 put 'The SURVEYMEANS SUBGROUP Macro';
 put;
 put;
 put '            Data Summary';
 put;
 &nobsrow;
 put '#';
 put 'Number of Strata:#' nstrata;
 put 'Number of PSUs:#' npsu;
 put 'Denominator Degrees of Freedom:#' ddf;
 format wsum 12.0;
 if (noteclupop) then do;
 put;
 put"NOTE:  Finite population correction to sample PSUs (&cluster)";
 put"assumes inclusion of all listing units within sample PSUs.";
 end;
 run;
%end;
%else %do;
data _null_;
merge &outsummary _temp;
file print;
 put "&usertitle";
 put;
 put 'The SURVEYMEANS SUBGROUP Macro';
 put;
 put;
 put '            Data Summary';
 put;
 put "&subtitle #";
 put '#';
 &nobsrow;
 put 'Subpopulation Observations:#' subpopobs '#Weight Sum:#' subpopwsum;
 put '#';
 put 'Number of Strata:#' nstrata;
 put 'Number of PSUs:#' npsu;
 put 'Denominator Degrees of Freedom:#' ddf;
 format wsum subpopwsum 12.0;
 if (noteclupop) then do;
 put;
 put"NOTE:  Finite population correction to sample PSUs (&cluster)";
 put"assumes inclusion of all listing units within sample PSUs.";
 end;
 run;
%end;
*------------------*
* contrast         *
*------------------*;
%if (&icontrast) %then %do;
data _null_;
set &outcontrast;
file print;
if _n_=1 then do;
 put '            Contrast Vector(s)';
 put;
 put 'Label#Variable#Values#Coefficients';
end;
put contrastlabel '#' contrastvar '#' values '#' coefficients;
run;
%end;
%end;
*------------------*
* tables           *
*------------------*;
 *---putlist,outlist---*;
%let list= &grouplist &outcontrastlabel &finalstats;
%let outlist=;
%let putlist=put;
%let outindex=1;
%let outputvar=%scan(&list,&outindex);
%do %while (%length(&outputvar));
 %let outlist= &outlist &outputvar #;
 %let putlist= &putlist &outputvar '#';
 %let outindex=%eval(&outindex+1);
 %let outputvar=%scan(&list,&outindex);
%end;
%let length1=%length(&outlist);
%let length2=%length(&putlist);
%let outlist=%substr(&outlist,1,&length1-1);
%let putlist=%substr(&putlist,1,&length2-3);
%let outlist=%sysfunc(compress(&outlist));

 *---print tables---*;
%do t=1 %to &nstatvar;
 %if (&t=1) %then %do;
  %let tablerequest=%sysfunc(tranwrd(&tablerequest,#,*));
  %if (not &icontrast) %then
        %let tabletitle=Table:  &tablerequest;
  %else %let tabletitle=Table:  &tablerequest *Contrast;
  data _null_;
  file print;
  put "&tabletitle";
  run;
 %end;
 %if (&ivar) %then %do;
  %let analvar=%scan(&var,&t);
  %let statvartitle=Analysis Variable:  &analvar;
  %end;
 %else %do;
  proc sort data=_toutput out=_fname nodupkey; by ratioindex; run;
  data _null_;
  set _fname;
  if (_n_=&t) then call symput('analratio',rationame);
  run;
  %let statvartitle=Ratio:              &analratio;
 %end;
 %if not (&icontrast) %then %do;
  data _null_;
  set _toutput;
  file print;
  if _n_=1 then do;
   put "&statvartitle";
   put;
   put "&outlist";
  end;
  &putlist;
  format &printformat &wordfmt;
  where &indextype=&t;
 %end;
 %else %if (&icontrast) %then %do;
  data _null_;
  set _toutput;
  file print;
  if _n_=1 then do;
   put "&statvartitle";
   put;
   put "&outlist";
  end;
  &putlist;
  format &printformat &wordfmt;
  where &indextype=&t;
 %end;
%end;
%end;

%mend print;

*--------------------------------------------------------------*
*                                                              *
*  TABLE                                                       *
*                                                              *
*--------------------------------------------------------------*;
%macro table;
%let star1=%index(&tablerequest,#);
%if (&star1) %then %do;
 %let temp=%substr(&tablerequest,&star1+1);
 %let star2=%index(&temp,#);
%end;
%if (&tablerequest=_overall_) %then %do;
 %let tabletype=0;
 %let numsubvars=1;
 %let _overall_=1;
 %let subgrp1=_overall_;
 %let subgrp2=-99;
 %let subgrp3=-99;
 %let grouplist=_OVERALL_;
 %let printtablerequest=_overall_;
 %let outsubgrp1=_overall_;
 %let outsubgrp2=ignore;
 %let outsubgrp3=ignore;
%end;
%else %if (not &star1) %then %do;
 %let tabletype=1;
 %let numsubvars=1;
 %let _overall_=0;
 %let subgrp1=&tablerequest;
 %if (&subgrp1=_overall_) %then %let _overall_=1;
 %let subgrp2=-99;
 %let subgrp3=-99;
 %let grouplist=&subgrp1;
 %let printtablerequest=&subgrp1;
 %let outsubgrp1=&subgrp1;
 %let outsubgrp2=ignore;
 %let outsubgrp3=ignore;
%end;
%else %if (not &star2) %then %do;
 %let tabletype=2;
 %let numsubvars=2;
 %let _overall_=0;
 %let subgrp1=%substr(&tablerequest,1,&star1-1);
 %let subgrp2=%substr(&tablerequest,&star1+1);
 %let subgrp3=-99;
 %let grouplist=&subgrp1 &subgrp2;
 %let printtablerequest=&subgrp1*&subgrp2;
 %let outsubgrp1=&subgrp1;
 %let outsubgrp2=&subgrp2;
 %let outsubgrp3=ignore;
%end;
%else %if (&star2) %then %do;
 %let tabletype=3;
 %let numsubvars=3;
 %let _overall_=0;
 %let subgrp3=%substr(&tablerequest,1,&star1-1);
 %let subgrp1=%substr(&temp,1,&star2-1);
 %let subgrp2=%substr(&temp,&star2+1);
 %let grouplist= &subgrp3 &subgrp1 &subgrp2;
 %let printtablerequest=&subgrp3, &subgrp1*&subgrp2;
 %let outsubgrp1=&subgrp1;
 %let outsubgrp2=&subgrp2;
 %let outsubgrp3=&subgrp3;
%end;
 *---formats---*;
%let s1fmt=_commona; %let s2fmt=_commonb; %let s3fmt=_commonc;
%let s1fmtcodes=; %let s2fmtcodes=; %let s3fmtcodes=;
%if %length(&format) %then %do;
 %let i=1;
 %let tformat= %scan(&format,1,&);
 %do %while (%length(&tformat));
  %let tsubgrp=%upcase(%scan(&tformat,1));
  %let subgrplength=%length(&tsubgrp);
  %let tfmtsubgrp=%substr(&tformat,%eval(&subgrplength+1));
  %if (&tsubgrp=%upcase(&subgrp1)) %then %do;
   %let s1fmt=&tsubgrp._;
   %let s1fmtcodes=&tfmtsubgrp;
  %end;
  %else %if (&tsubgrp=%upcase(&subgrp2)) %then %do;
   %let s2fmt=&tsubgrp._;
   %let s2fmtcodes=&tfmtsubgrp;
  %end;
  %else %if (&tsubgrp=%upcase(&subgrp3)) %then %do;
   %let s3fmt=&tsubgrp._;
   %let s3fmtcodes=&tfmtsubgrp;
  %end;
  %let i=%eval(&i+1);
  %let tformat=%scan(&format,&i,&);
 %end;
%end;
%mend table;

*--------------------------------------------------------------*
*                                                              *
*  RATIO                                                       *
*                                                              *
*--------------------------------------------------------------*;
%macro ratio;
 *---initialize---*;
%let topvar=; %let topequal=;  %let topvalue=;
%let botvar=;  %let botequal=;  %let botvalue=;
%let ratio1=; %let ratio2=; %let ratio3=; %let ratio4=; %let ratio5=;
%let ratio6=; %let ratio7=; %let ratio8=; %let ratio9=; %let ratio10=;

 *---process---*;
%let i=1;
%let tratio=%scan(&ratio,1,&);
%do %while (%length(&tratio) and (&i<=10));
 %let rslash=%index(&tratio,/);
 %let leftside=%substr(&tratio,1,&rslash-1);
 %let rightside=%substr(&tratio,&rslash+1);
 %let ttopequal= %index(&leftside,=);
 %let tbotequal= %index(&rightside,=);
 %if (&ttopequal) %then %do;
  %let ttopvar= %substr(&leftside,1,&ttopequal-1);
  %let ttopvalue=%substr(&leftside,&ttopequal+1);
 %end;
 %else %do;
  %let ttopvar=&leftside;
  %let ttopvalue=0;
 %end;
 %if (&tbotequal) %then %do;
  %let tbotvar= %substr(&rightside,1,&tbotequal-1);
  %let tbotvalue=%substr(&rightside,&tbotequal+1);
 %end;
 %else %do;
  %let tbotvar=&rightside;
  %let tbotvalue=0;
 %end;
 %if (&i=1) %then       %let ratio1=&tratio;
 %else %if (&i=2) %then %let ratio2=&tratio;
 %else %if (&i=3) %then %let ratio3=&tratio;
 %else %if (&i=4) %then %let ratio4=&tratio;
 %else %if (&i=5) %then %let ratio5=&tratio;
 %else %if (&i=6) %then %let ratio6=&tratio;
 %else %if (&i=7) %then %let ratio7=&tratio;
 %else %if (&i=8) %then %let ratio8=&tratio;
 %else %if (&i=9) %then %let ratio9=&tratio;
 %else %if (&i=10) %then %let ratio10=&tratio;
 %let topvar= &topvar &ttopvar;
 %let topequal= &topequal &ttopequal;
 %let topvalue= &topvalue &ttopvalue;
 %let botvar= &botvar &tbotvar;
 %let botequal= &botequal &tbotequal;
 %let botvalue= &botvalue &tbotvalue;
 %let i=%eval(&i+1);
 %let tratio=%scan(&ratio,&i,&);
%end;
%mend ratio;

*--------------------------------------------------------------*
*                                                              *
*  CONTRAST                                                    *
*                                                              *
*--------------------------------------------------------------*;
%macro contrast;
 *---initialize---*;
%let cvar1=;%let cvar2=; %let cvar3=; %let cvar4=; %let cvar5=;
%let cvar6=;%let cvar7=; %let cvar8=; %let cvar9=; %let cvar10=;
%let clab1=;%let clab2=; %let clab3=; %let clab4=; %let clab5=;
%let clab6=;%let clab7=; %let clab8=; %let clab9=; %let clab10=;
 *---process---*;
%let tcontrast=%scan(&contrast,1,&);
%let i=1; %let maxlength=0;
%do %while(%length(&tcontrast) and (&i<=10));
 %let cont2=%sysfunc(translate(&tcontrast,"$","'"));
 %let cont2=%substr(&cont2,2);
 %let dollar=%index(&cont2,$);
 %let clab=%substr(&cont2,1,&dollar-1);
 %let cont3=%substr(&cont2,&dollar+1);
 %let cvar=%scan(&cont3,1);
 %let firstcoef=%scan(&cont3,2,%str( ));
 %let pfirstcoef=%index(&cont3,&firstcoef);
 %let coef=%substr(&cont3,&pfirstcoef);
 %if &i=1 %then %do;
   %let cvar1=&cvar; %let clab1=&clab; %let coef1=&coef; %end;
 %else %if &i=2 %then %do;
   %let cvar2=&cvar; %let clab2=&clab; %let coef2=&coef; %end;
 %else %if &i=3 %then %do;
   %let cvar3=&cvar; %let clab3=&clab; %let coef3=&coef; %end;
 %else %if &i=4 %then %do;
   %let cvar4=&cvar; %let clab4=&clab; %let coef4=&coef; %end;
 %else %if &i=5 %then %do;
   %let cvar5=&cvar; %let clab5=&clab; %let coef5=&coef; %end;
 %else %if &i=6 %then %do;
   %let cvar6=&cvar; %let clab6=&clab; %let coef6=&coef; %end;
 %else %if &i=7 %then %do;
   %let cvar7=&cvar; %let clab7=&clab; %let coef7=&coef; %end;
 %else %if &i=8 %then %do;
   %let cvar8=&cvar; %let clab8=&clab; %let coef8=&coef; %end;
 %else %if &i=9 %then %do;
   %let cvar9=&cvar; %let clab9=&clab; %let coef9=&coef; %end;
 %else %if &i=10 %then %do;
   %let cvar10=&cvar; %let clab10=&clab; %let coef10=&coef; %end;
 %let tlength=%length(&cvar);
 %let maxlength=%sysfunc(max(&tlength,&maxlength));
 %let i=%eval(&i+1);
 %let tcontrast=%scan(&contrast,&i,&);
%end;
%let cvar=&cvar1 &cvar2 &cvar3 &cvar4 &cvar5 &cvar6 &cvar7 &cvar8
          &cvar9 &cvar10;
%mend contrast;

*--------------------------------------------------------------*
*                                                              *
*  SUBPOP                                                      *
*                                                              *
*--------------------------------------------------------------*;
%macro subpop;
 %let subpop2=%upcase(&subpop);
 %if %sysfunc(indexc(&subpop2, '(' )) and %sysfunc(indexc(&subpop2,')'))
   %then %do;
    %let inoperator=1;
    %let subpoplabel=Subpop='Subpopulation';
    %let subpopvar=%scan(&subpop2,1);
    %let leftpar= %sysfunc(indexc(&subpop,'('));
    %let rightpar=%sysfunc(indexc(&subpop,')'));
    %let numpos=%eval(&rightpar-&leftpar-1);
    %let inelements=%substr(&subpop2,&leftpar+1,&numpos);
    %let subpop2=ignore=1;
    %end;
 %else %do;
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str (=), %str( = )  ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str (^), %str( ^ )  ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str (>), %str( > )  ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str (<), %str( < )  ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str (~), %str( ^ )  ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( >  = ), %str( >= ) ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( <  = ), %str( <= ) ));
  %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( ^  = ), %str( ^= ) ));
  %let word2=%scan(&subpop,2);
  %if (&word2=EQ) or (&word2=NE) or (&word2=GT) or (&word2=LT) or
      (&word2=GE) or (&word2=LE) %then %do;
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( EQ ), %str( = )  ));
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( NE ), %str( ^= ) ));
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( GT ), %str( > )  ));
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( LT ), %str( < )  ));
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( GE ), %str( >= ) ));
    %let subpop2=%sysfunc(tranwrd(&subpop2,%str ( LE ), %str( <= ) ));
    %let subpoplabel=Subpop='Subpopulation';
    %let subpopvar=%scan(&subpop2,1);
    %end;
  %if  (%sysfunc(indexc(&subpop2,'~=^><'))=0) %then %let ERR301=1;
 %end;
%mend subpop;

*--------------------------------------------------------------*
*                                                              *
*  STATISTICS                                                  *
*                                                              *
*--------------------------------------------------------------*;
%macro statistics;
*---------------------------*
* process parameter         *
*---------------------------*;
 *---initialize---*;
%let fnobs=12.0; %let fnmiss=12.0; %let fsumwgt=12.1; %let fmean=12.4;
%let fsemean=12.4; %let fclmean=12.4; %let fcvmean=12.4; %let ftotal=12.1;
%let fsetotal=12.1; %let fcltotal=12.1; %let fcvtotal=12.4;
%let fnumtotal=12.1; %let fdentotal=12.1; %let fratio=12.4;
%let fseratio=12.4; %let fclratio=12.4; %let fcvratio=12.4;
%let ftmean=12.2; %let ftratio=12.2; %let fpmean=12.4; %let fpratio=12.4;
%let noprint=0;

 *---options---*;
%let tstatistics=%str( ) %upcase(&statistics) %str( );
%let slash=%index(&tstatistics,/);
%if (&slash) %then %do;
 %let statoptions= %substr(&tstatistics,&slash+1);
 %let tstatistics= %substr(&tstatistics,1,&slash-1);
 %let noprint= %index(&statoptions,NOPRINT);
 %let statoptions=%sysfunc(tranwrd(&statoptions,NOPRINT,%str( )));
 %let equal= %index(&statoptions,=);
 %if (&equal) %then %do;
   *---first format---*;
  %let leftside=%upcase(%substr(&statoptions,1,&equal-1));
  %let rest=%substr(&statoptions,&equal+1);
  %let rightside=%scan(&rest,1,%str( ) );
  %if (&ivar) %then %varfmt; %else %ratiofmt;
   *---remaining formats---*;
  %let rest=&rest z_z_z;
  %let nextword=%scan(&rest,2,%str( ) );
  %do %while (&nextword ne z_z_z);
   %let pnextword=%index(&rest,&nextword);
   %let truncate=%substr(&rest,&pnextword);
   %let equal= %index(&truncate,=);
   %let leftside=%upcase(%substr(&truncate,1,&equal-1));
   %let rest=%substr(&truncate,&equal+1);
   %let rightside=%scan(&rest,1,%str( ) );
   %if (&ivar) %then %varfmt; %else %ratiofmt;
   %let nextword=%scan(&rest,2,%str( ) );
  %end;
   *---remove z_z_z from last format---*;
  %let pz_z_z=%index(&rest,z_z_z);
  %let rightside=%substr(&rest,1,&pz_z_z-1);
  %if (&ivar) %then %varfmt; %else %ratiofmt;
 %end;
%end;
*---------------------------*
* check parameter           *
*---------------------------*;
%let unknown=;
%if %length(&tstatistics) and (&tstatistics ne ALL) %then %do;
 %let i=1;
 %let keyword=%scan(&tstatistics,&i);
 %do %until (%length(&keyword)=0);
 %if (&keyword ne NOBS) and (&keyword ne NMISS) and (&keyword ne SUMWGT)and
   (&keyword ne MEAN) and (&keyword ne SEMEAN)and (&keyword ne CLMEAN) and
   (&keyword ne CVMEAN) and (&keyword ne TOTAL) and
   (&keyword ne SETOTAL) and (&keyword ne CLTOTAL) and
   (&keyword ne CVTOTAL) and (&keyword ne NUMTOTAL) and
   (&keyword ne DENTOTAL) and (&keyword ne RATIO) and
   (&keyword ne SERATIO) and (&keyword ne CLRATIO) and
   (&keyword ne CVRATIO) and (&keyword ne TMEAN) and (&keyword ne PMEAN) and
   (&keyword ne TRATIO) and (&keyword ne PRATIO)
   %then %do;
     %let ERR203=1;
     %let unknown=&unknown %scan(&tstatistics,&i);
     %end;
  %if (&ivar) %then %do;
     %if (&keyword=RATIO) or (&keyword=SERATIO) or (&keyword=CLRATIO) or
     (&keyword=CVRATIO) or (&keyword=NUMTOTAL) or (&keyword=DENTOTAL) or
     (&keyword=TRATIO) or (&keyword=PRATIO)
        %then %let  ERR201=1;
  %end;
  %else %if (&iratio) %then %do;
     %if (&keyword=MEAN) or (&keyword=SEMEAN) or (&keyword=CLMEAN) or
     (&keyword=CVMEAN) or (&keyword=TOTAL) or (&keyword=SETOTAL) or
     (&keyword=CLTOTAL) or (&keyword=CVTOTAL) or (&keyword=TMEAN) or
     (&keyword=PMEAN)
       %then %let ERR202=1;
  %end;
  %if (&icontrast=0) %then %do;
     %if (&keyword=TMEAN) or (&keyword=PMEAN) or (&keyword=TRATIO) or
         (&keyword=PRATIO)

        %then %let  ERR204=1;
  %end;
  %let i=%eval(&i+1);
  %let keyword=%scan(&tstatistics,&i);
  %end;
%end;

*---------------------------*
* create finalstats         *
*---------------------------*;
 *---statistics=all---*;
%if (&tstatistics=ALL) %then %do;
 %if (&ivar=1 and &icontrast=0) %then
   %let tstatistics=%str( NOBS NMISS SUMWGT MEAN SEMEAN CLMEAN CVMEAN
TOTAL SETOTAL CLTOTAL CVTOTAL );
 %else %if (&ivar=1 and &icontrast=1) %then
   %let tstatistics=%str( NOBS NMISS SUMWGT MEAN SEMEAN CLMEAN CVMEAN TMEAN
PMEAN TOTAL SETOTAL CLTOTAL CVTOTAL );
 %else %if (&iratio=1 and &icontrast=0) %then
   %let tstatistics=%str( NOBS NMISS SUMWGT NUMTOTAL DENTOTAL RATIO SERATIO
CLRATIO CVRATIO );
 %else %if (&iratio=1 and &icontrast=1) %then
   %let tstatistics=%str( NOBS NMISS SUMWGT NUMTOTAL DENTOTAL RATIO SERATIO
CLRATIO CVRATIO TRATIO PRATIO );
 %end;

 *---statistics=blank---*;
%if (%length(&tstatistics)=0) %then %do;
 %if (&ivar=1 and &icontrast=0) %then
       %let tstatistics= %str( ) NOBS MEAN SEMEAN CLMEAN %str( );
 %else %if (&ivar=1 and &icontrast=1) %then
       %let tstatistics=%str( ) NOBS MEAN SEMEAN TMEAN PMEAN %str( );
 %else %if (&iratio=1 and &icontrast=0) %then
       %let tstatistics=%str( ) NOBS RATIO SERATIO CLRATIO %str( );
 %else %if (&iratio=1 and &icontrast=1) %then
       %let tstatistics=%str( ) NOBS RATIO SERATIO TRATIO PRATIO %str( );
%end;

 *---process clmean,cltotal,clratio---*;
%let finalstats=%str( )
  %sysfunc(tranwrd(&tstatistics,
      %str( CLMEAN ),%str( LCLMEAN UCLMEAN ))) %str( );
%let finalstats=%str( )
  %sysfunc(tranwrd(&finalstats,
     %str( CLTOTAL ),%str( LCLTOTAL UCLTOTAL ))) %str( );
%let finalstats=%str( )
  %sysfunc(tranwrd(&finalstats,
     %str( CLRATIO ),%str( LCLRATIO UCLRATIO ))) %str( );

 *--- finalstats_out ---*;
%let finalstats_out=&finalstats;

 *--- finalstats_print ---*;
%let finalstats_print=%str( ) &finalstats %str( );
%let finalstats_print=%str( ) %sysfunc(tranwrd(&finalstats_print,
       %str( NMISS ),%str( NMISS_P ))) %str( );
%let finalstats_print=%str( ) %sysfunc(tranwrd(&finalstats_print,
       %str( SUMWGT ),%str( SUMWGT_P ))) %str( );
%let finalstats_print=%str( ) %sysfunc(tranwrd(&finalstats_print,
       %str( MEAN ),%str( MEAN_P ))) %str( );

 *--- finalstats_both ---*;
%let finalstats_both=%str( ) &finalstats %str( );
%let finalstats_both=%str( ) %sysfunc(tranwrd(&finalstats_both,%str( NMISS ),
       %str( NMISS NMISS_P ))) %str( );
%let finalstats_both=%str( ) %sysfunc(tranwrd(&finalstats_both,%str( SUMWGT ),
       %str( SUMWGT SUMWGT_P ))) %str( );
%let finalstats_both=%str( ) %sysfunc(tranwrd(&finalstats_both,%str( MEAN ),
       %str( MEAN MEAN_P ))) %str( );

*---------------------------*
* print formats             *
*---------------------------*;
%if (&noprint=1) %then %let finalstats_print2=ignore;
%else %do;
 *---finalstats_print2---*;
 %let finalstats_print2=%str( ) &finalstats_print %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
     %str( NOBS ),%str( NOBS*f=&fnobs ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
     %str( NMISS_P ),%str( NMISS_P*f=&fnmiss ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
     %str( SUMWGT_P ),%str( SUMWGT_P*f=&fsumwgt ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( MEAN_P ),%str( MEAN_P*f=&fmean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( SEMEAN ),%str( SEMEAN*f=&fsemean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( LCLMEAN ),%str( LCLMEAN*f=&fclmean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( UCLMEAN ),%str( UCLMEAN*f=&fclmean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( CVMEAN ),%str( CVMEAN*f=&fcvmean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( TOTAL ),%str( TOTAL*f=&ftotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( SETOTAL ),%str( SETOTAL*f=&fsetotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( LCLTOTAL ),%str( LCLTOTAL*f=&fcltotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( UCLTOTAL ),%str( UCLTOTAL*f=&fcltotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( CVTOTAL ),%str( CVTOTAL*f=&fcvtotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( NUMTOTAL ),%str( NUMTOTAL*f=&fnumtotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( DENTOTAL ),%str( DENTOTAL*f=&fdentotal ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( RATIO ),%str( RATIO*f=&fratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( SERATIO ),%str( SERATIO*f=&fseratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( LCLRATIO ),%str( LCLRATIO*f=&fclratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( UCLRATIO ),%str( UCLRATIO*f=&fclratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( CVRATIO ),%str( CVRATIO*f=&fcvratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( TMEAN ),%str( TMEAN*f=&ftmean ))) %str( );
  %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( PMEAN ),%str( PMEAN*f=&fpmean ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( TRATIO ),%str( TRATIO*f=&ftratio ))) %str( );
 %let finalstats_print2= %str( ) %sysfunc(tranwrd(&finalstats_print2,
      %str( PRATIO ),%str( PRATIO*f=&fpratio ))) %str( );

  *---wordfmt---*;
 %let wordfmt=%str( ) &finalstats %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( NOBS ),%str( NOBS &fnobs ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( NMISS ),%str( NMISS &fnmiss ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( SUMWGT ),%str( SUMWGT &fsumwgt ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( MEAN ),%str( MEAN &fmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( SEMEAN ),%str( SEMEAN &fsemean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( LCLMEAN ),%str( LCLMEAN &fclmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( UCLMEAN ),%str( UCLMEAN &fclmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( CVMEAN ),%str( CVMEAN &fcvmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( TOTAL ),%str( TOTAL &ftotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( SETOTAL ),%str( SETOTAL &fsetotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( LCLTOTAL ),%str( LCLTOTAL &fcltotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( UCLTOTAL ),%str( UCLTOTAL &fcltotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( CVTOTAL ),%str( CVTOTAL &fcvtotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( NUMTOTAL ),%str( NUMTOTAL &fnumtotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( DENTOTAL ),%str( DENTOTAL &fdentotal ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( RATIO ),%str( RATIO &fratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( SERATIO ),%str( SERATIO &fseratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( LCLRATIO ),%str( LCLRATIO &fclratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( UCLRATIO ),%str( UCLRATIO &fclratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( CVRATIO ),%str( CVRATIO &fcvratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( TMEAN ),%str( TMEAN &ftmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( PMEAN ),%str( PMEAN &fpmean ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( TRATIO ),%str( TRATIO &ftratio ))) %str( );
 %let wordfmt=%sysfunc(tranwrd(&wordfmt,
      %str( PRATIO ),%str( PRATIO &fpratio ))) %str( );
%end;
%mend statistics;

 *---VARFMT---*;
%macro varfmt;
%if (&leftside=FNOBS)          %then %let fnobs=&rightside;
%else %if (&leftside=FNMISS)   %then %let fnmiss=&rightside;
%else %if (&leftside=FSUMWGT)  %then %let fsumwgt=&rightside;
%else %if (&leftside=FMEAN)    %then %let fmean=&rightside;
%else %if (&leftside=FSEMEAN)  %then %let fsemean=&rightside;
%else %if (&leftside=FTOTAL)   %then %let ftotal=&rightside;
%else %if (&leftside=FSETOTAL) %then %let fsetotal=&rightside;
%else %if (&leftside=FCLMEAN)  %then %let fclmean=&rightside;
%else %if (&leftside=FCVMEAN)  %then %let fcvmean=&rightside;
%else %if (&leftside=FCLTOTAL) %then %let fcltotal=&rightside;
%else %if (&leftside=FCVTOTAL) %then %let fcvtotal=&rightside;
%else %if (&leftside=FTMEAN)   %then %let ftmean=&rightside;
%else %if (&leftside=FPMEAN)   %then %let fpmean=&rightside;
%mend varfmt;
 *---RATIOFMT---*;
%macro ratiofmt;
%if (&leftside=FNOBS)           %then %let fnobs=&rightside;
%else %if (&leftside=FNMISS)    %then %let fnmiss=&rightside;
%else %if (&leftside=FSUMWGT)   %then %let fsumwgt=&rightside;
%else %if (&leftside=FNUMTOTAL) %then %let fnumtotal=&rightside;
%else %if (&leftside=FDENTOTAL) %then %let fdentotal=&rightside;
%else %if (&leftside=FRATIO)    %then %let fratio=&rightside;
%else %if (&leftside=FSERATIO)  %then %let fseratio=&rightside;
%else %if (&leftside=FCLRATIO)  %then %let fclratio=&rightside;
%else %if (&leftside=FCVRATIO)  %then %let fcvratio=&rightside;
%else %if (&leftside=FTRATIO)   %then %let ftratio=&rightside;
%else %if (&leftside=FPRATIO)   %then %let fpratio=&rightside;
%mend ratiofmt;

*--------------------------------------------------------------*
*                                                              *
*  GLOBAL                                                      *
*                                                              *
*--------------------------------------------------------------*;
%global
finalstats finalstats_out finalstats_print finalstats_both finalstats_print2
wordfmt noprint unknown topvar topequal topvalue botvar botequal botvalue
ratio1 ratio2 ratio3 ratio4 ratio5 ratio6 ratio7 ratio8 ratio9 ratio10
cvar cvar1 cvar2 cvar3 cvar4 cvar5 cvar6 cvar7 cvar8 cvar9 cvar10
clab1 clab2 clab3 clab4 clab5 clab6 clab7 clab8 clab9 clab10
coef1 coef2 coef3 coef4 coef5 coef6 coef7 coef8 coef9 coef10 maxlength
subpoplabel subpopvar subpop2
tabletype numsubvars grouplist printtablerequest
outsubgrp1 outsubgrp2 outsubgrp3 _overall_ subgrp1 subgrp2 subgrp3
s1fmt s2fmt s3fmt s1fmtcodes s2fmtcodes s3fmtcodes;