*****************************************************************************************;
* Macro Name: WtdCatmod (File, NRepWts, NRMeas, RespVars, NCov, Covs, Effects)          *;
*                                                                                       *;
* Input:     File = Name of input SAS dataset which includes:                           *;
*                       the categorical response variables, coded as 1,2,...Ncat        *;
*                       the covariates, coded as 1 or 2                                 *;
*                       the weights, named wt0, wt1, wt2, ..., wtNRepWts                *;
*                       the full sample weight is wt0, the others are replicate weights *;
*                       the value of Ncat is determined by the program                  *;
*                       each response variable must have Ncat categories                *;
*         NRepWts = Number of sets of replicate weights                                 *;
*         NCov    = Number of binary covariates                                         *;
*         Covs    = List of names of the covariates, separated by spaces                *;
*         NRMeas   = Number of repeated measurements                                    *;
*         RespVars = List of response variable names, separated by spaces               *;
*         Effects  = List of effects to be included in the model, separated by spaces   *;
*                                                                                       *;
* Purpose:    Runs PROC CATMOD taking the survey design (delete-cluster jackknife)into  *;
*             account, by automating the construction of the direct input data set and  *;
*             the user-supplied design matrix                                           *;
*****************************************************************************************;

%macro WtdCatmod (File, NRepWts, NRMeas, RespVars, NCov, Covs, Effects);

*---------------------------------------------------------------------------------------*;
* Make macro variables of the response variables:  ResponseT1, ... ResponseT&NRMeas     *;
*---------------------------------------------------------------------------------------*;

data _null;
%do i=1 %to &NRMeas;
   call symput("ResponseT&i", scanq("&RespVars", &i));
   %end;
run;

*---------------------------------------------------------------------------------------*;
* Find out how many categories there are in first response variable                     *;
*---------------------------------------------------------------------------------------*;

proc freq data=&File noprint;
tables &ResponseT1 / out = NoOfLevels;
run;

data _null_;
set NoOfLevels nobs=TotalObs;
if (_n_ = 1) then call symput("NCatsResponse", trim(left(TotalObs)));
run;

*---------------------------------------------------------------------------------------*;
* Assign constants to macro variables that are used throughout the program              *;
*---------------------------------------------------------------------------------------*;

%let NIndVars = %eval(2 * (&NCatsResponse - 1));    /* Number of indicator variables    */
%let CovMatxSize = %eval(&NIndVars ** 2);        /* Dimension of cov matrix, within pop */
%let NPops = %eval(2 ** &NCov);                     /* Number of populations            */
%let Dim = %eval(&NPops * &NIndVars);            /* Overall covariance matrix dimension */

*---------------------------------------------------------------------------------------*;
* Create indicator variables for each response variable                                 *;
*---------------------------------------------------------------------------------------*;

data FmtdData;
   set &File;
   array IndVar (&NIndVars);
   %do i=1 %to &NRMeas;
      %do j=1 %to &NCatsResponse-1;
         IndVar((&NCatsResponse-1) * (&i - 1) + &j) = 0;
         if (&&ResponseT&i = &j) then
            IndVar((&NCatsResponse-1) * (&i - 1) + &j) = 1;
      %end;
   %end;
run;

*---------------------------------------------------------------------------------------*;
* Calculate the mean of each response variable using the full sample weights and each   *;
* set of replicate weights                                                              *;
*---------------------------------------------------------------------------------------*;

proc sort data = FmtdData;
by &Covs;
run;

%do I = 0 %to &NRepWts;

   proc means data = FmtdData noprint;
   var IndVar1 - IndVar&NIndVars;
   weight Wt&i;
   by &Covs;
   output
      out = MeansRep&i (drop = _type_ _freq_)
      mean = MeanIndVar1 - MeanIndVar&NIndVars;
   run;

*---------------------------------------------------------------------------------------*;
* Create variable Rep to indicate the replicate weight                                  *;
*---------------------------------------------------------------------------------------*;

   %if (&i ^= 0) %then %do;
      data MeansRep&i;
      set MeansRep&i;
      Rep = &i;
      run;
   %end;

%end;

*---------------------------------------------------------------------------------------*;
* Append datasets of the means calculated from each set of replicate weights            *;
*---------------------------------------------------------------------------------------*;

data MeansAllReps;
set %do i = 1 %to &NRepWts;  MeansRep&I  %end;;
by &Covs;
run;

*---------------------------------------------------------------------------------------*;
* Clean up extraneous datasets                                                          *;
*---------------------------------------------------------------------------------------*;

proc datasets nodetails nolist;
delete %do i = 1 %to &NRepWts;  MeansRep&i  %end;;
run;

*---------------------------------------------------------------------------------------*;
* Evaluate the difference (Delta) between the means using each set of replicate weights *;
* and the mean using the full sample weight, to subsequently evaluate the covariance    *; 
* matrix within each population (CovWp) for each set of replicate weights               *;
*---------------------------------------------------------------------------------------*;

data VectorProd;
merge
   MeansAllReps
   MeansRep0 (rename = (%do i=1 %to &NIndVars; MeanIndVar&i = SampMean&i  %end;));
by &Covs;
array MeanIndVar [&NIndVars];
array SampMean [&NIndVars];
array Delta [&NIndVars];
VarMult =  (&NRepWts - 1) / &NRepWts;
do i=1 to &NIndVars;
   Delta[i] = MeanIndVar[i] - SampMean[i];
   end;

array CovWp[&NIndVars, &NIndVars] CovWp1-CovWp&CovMatxSize;
do i = 1 to &NIndVars;
   do j = 1 to &NIndVars;
      CovWp[i,j] = Delta[i] * Delta[j] * VarMult;
      end;
   end;
drop i j Rep MeanIndVar1-MeanIndVar&NIndVars SampMean1-SampMean&NindVars
   Delta1-Delta&NIndVars;
run;

*---------------------------------------------------------------------------------------*;
* Sum the covariance matrices within the populations over the sets of replicate weights *;
*---------------------------------------------------------------------------------------*;

proc means   data = VectorProd   noprint;
var CovWp1 - CovWp&CovMatxSize;
by &Covs;
output
   out = cov(drop = _type_ _freq_)
   sum = Cov1 - Cov&CovMatxSize;
run;

*---------------------------------------------------------------------------------------*;
* Creates the cov rows for the dataset by distributing the covariance matrices within   *;
* the populations appropriately into the overall covariance matrix                      *;
*---------------------------------------------------------------------------------------*;

data CovMatrix;
length b1-b&Dim 8. _type_ _name_ $5.;
set cov;
_type_ = 'cov';
%do i = 1 %to &Dim;  b&i = 0;  %end;
%do k = 1 %to &NPops;
   if (_n_ = &k) then do;
   %do i=1 %to &NIndVars;
      %do j=1 %to &NIndVars;
         b%eval(((&k-1)*&NIndVars) + &j) = cov%eval(((&i-1)*&NIndVars) + &j);
      %end;
      _name_ = "b%eval(((&k-1)*&NIndVars) + &i)";
      output;
   %end;
   end;
%end;
drop Cov1-Cov&CovMatxSize;
run;

*---------------------------------------------------------------------------------------*;
* Creates the parms rows for the dataset by distributing the response function          *;
* estimates within the populations appropriately                                        *;
*---------------------------------------------------------------------------------------*;

data SampMeansAcross;
   set MeansRep0 nobs=TotalObs;
   array b[&Dim];
   retain b1--b&Dim;
   array MeanIndVar[&NIndVars];
   do i=1 to &NIndVars;
      b[((_n_ - 1) * &NIndVars) + i] = MeanIndVar[i];
      end;
   _type_ = 'parms';
   _name_ = '   ';
   if (_n_= TotalObs);
   keep b1-b&Dim _type_ _name_;
run;

*---------------------------------------------------------------------------------------*;
* Combines the parms row with the cov rows, saves the dataset for later input to        *;
* PROC CATMOD                                                                           *;
*---------------------------------------------------------------------------------------*;

data Matrix4Catmod;
set SampMeansAcross CovMatrix;
drop &Covs;
run;

*---------------------------------------------------------------------------------------*;
* Run PROC CATMOD with unweighted data to obtain the design matrix and hypothesis test  *;
* specifications for the MODEL statement in final CATMOD                                *;
*---------------------------------------------------------------------------------------*;

ods output   ResponseDesign = DesMatrix;
ods output   Estimates = Est;
ods listing close;

proc catmod data=&File;
response marginals;
population &Covs;
model &ResponseT1  %do i=2 %to &NRMeas; * &&ResponseT&I  %end;
   = _response_ &Effects / design;
Repeated Response &NRMeas;
run;

ods output close;
ods listing;
quit;

*---------------------------------------------------------------------------------------*;
* Use the DesMatrix dataset to build macro variables, where each variable is one row of *;
* the design  matrix portion of the model statement, initialize the macro variable      *;
* indicating the number of rows in the design matrix                                    *;
*---------------------------------------------------------------------------------------*;

proc contents data=DesMatrix noprint out=VarsDesMatrix;
run;

data _null_;
   set VarsDesMatrix end=eof;
   retain count 0;
   if (upcase(substr(name, 1, 1)) = 'D') then
      if (0 <= substr(name, 2, 1) <= 9) then count + 1;
   if (eof) then call symput("DimDesMatrix", trim(left(count)));
run;

*---------------------------------------------------------------------------------------*;
* Notice that the "length" of DesMatrix must be at least three times the number of      *;
* columns in the design matrix.  If there are more than 100 columns (model parameters)  *;
* in your design matrix, simply increase the parameter in the length statement below.   *;
*---------------------------------------------------------------------------------------*;

data _null_;
set DesMatrix end=eof;
length DesMatrix $300.;
DesMatrix = trim(left(d1))
%do I = 2  %to &DimDesMatrix; || put(d&i, 4.)   %end;;
call symput("DesMatrix" || left(trim( _n_ )), DesMatrix);
if (eof) then call symput("NoDesMatrix", left(trim( _n_ )));
run;

*---------------------------------------------------------------------------------------*;
* Use the EST dataset to build macro variables, where each variable is a hypothesis     *;
* test specification for the model statement                                            *;
*---------------------------------------------------------------------------------------*;

proc sort data=Est;
by Effect;
run;

proc transpose
data=Est  out=EaCov (drop=_name_)  prefix=LVL;
var Parameter;
by Effect;
run;

proc contents data=EaCov noprint out=vars;
run;

data _null_;
set vars end=eof;
retain MaxCol 0;
if (substr(name, 1, 3) = 'LVL') then MaxCol + 1;
if (eof) then  call symput("MaxCol", trim(left(MaxCol)));
run;

data _null_;
set EaCov;
call symput("col" || trim(left(_n_)), trim(left(LVL1))
      %do i=2 %to &MaxCol;
      || ' ' || trim(left(LVL&i))
      %end;
   || ' = "' || trim(left(Effect)) || '"' );
call symput ("NoCol", trim(left(_n_)));
run;

*---------------------------------------------------------------------------------------*;
* Perform weighted CATMOD, using macro variables DESMATRIX1, DESMATRIX2,... to specify  *;
* the design matrix and macro variables COL1, COL2,... to specify the hypothesis tests  *;
* in the MODEL statement                                                                *;
*---------------------------------------------------------------------------------------*;

proc catmod   data = Matrix4Catmod;
response read b1-b&Dim;
model _f_ = (%do i=1 %to &NoDesMatrix-1;  &&DesMatrix&i,  %end;
   &&DesMatrix&NoDesMatrix)
   (%do i=1 %to &NoCol-1;  &&Col&i,  %end;  &&Col&i)
   / nodesign noparm;
run;
quit;

%mend WtdCatmod;

*---------------------------------------------------------------------------------------*;
* Invoke macro                                                                          *;
*---------------------------------------------------------------------------------------*;

%WtdCatmod (
   File = intensity,
   NRepWts = 84,
   NRMeas = 2,
   RespVars = Resp1 Resp2,
   NCov = 4,
   Covs = Age Educ Sex Quit,
   Effects = Age Educ Sex|Quit Quit|_response_);

*****************************************************************************************;