%macro par(bdata=, pdata=, vpdata=, vpvarname=_pv_,
	   n_or_p=p, n_or_pname=prev, 
	   modvar=, fixedvar=,
	   fullpar=F, partialpar=T,
           citransf=fisherz,
	   vprevopt=FIXED,
	   sampsz=,
           missvarlist=NONE, /** missing indicator list **/
	   notes=nonotes
    );
   

   /**********************************************
    varname&i  all names in v-c matrix for coeffs
    b_varname&i coeff of ith variable
    cvb_&i_&j  i,j th element of v-c matrix for coeffs
    **************************************************/

/** RL: sort it by _NAME_ so that when the estimate is sorted, the var-covar will change accordingly **/
%if %length(&bdata) ne 0 %then %do;
proc sort data=&bdata; by _NAME_;
%end;

proc optsave  out=_tmpoptsfromrun_;  run;
   options &notes;
   
%local i j k l q _nbet _nmodvar _nfixedvar _incomp _fixedvarp
   _nlev _overlap_ _varnameprb_ _rrprob_
   _varfull_ _varpart_ _fullpar_ _partlpar_ _errmv_
   ;
%global _parcallno_;
%if &_parcallno_ eq %then %do;  %let _parcallno_=1;  %end;
%else %do;  %let _parcallno_=%eval(&_parcallno_ + 1);  %end;
/*
title4 "PAR macro call number &_parcallno_";
*/
   
   /* changing all user macro vbls to upcase */
%let vprevopt=%upcase(&vprevopt);
%let n_or_p = %upcase(&n_or_p);
%if %length(&modvar) ne 0 %then %do;
  %let modvar=%upcase(&modvar);
%end;
%if %length(&fixedvar) ne 0 %then %do;
  %let fixedvar=%upcase(&fixedvar);
%end;
%let n_or_pname=%upcase(&n_or_pname);
%let fullpar=%upcase(&fullpar);
%let partialpar=%upcase(&partialpar);
%let citransf=%upcase(&citransf);
%if %length(&missvarlist) eq 0 %then %do;
  data _null_;  errmv=1;  call symput('_errmv_', errmv);  run;
  data _null_;
    put 'ERROR: in PAR macro call:  The parameter MISSVARLIST is empty.';
    put '  If your model has no missing indicators, you must ';
    put '  set missvarlist=NONE ';
    file print;
    put 'ERROR: in PAR macro call:  The parameter MISSVARLIST is empty.';
    put '  If your model has no missing indicators, you must ';
    put '  set missvarlist=NONE ';
  run;
%end;
%if &_errmv_ eq 1 %then %goto errend;

%let missvarlist=%upcase(&missvarlist);
%if &missvarlist eq NONE %then %do;  %let missvarlist= ;  %end;  /* now missvarlist can be empty */
%let _nmissvar=%numargs(&missvarlist);
%if &_nmissvar ne 0 %then %do;
  %do i=1 %to &_nmissvar;  %let _missvar&i = %scan(&missvarlist, &i, %str( ));  %end;
%end;

   
%if %length(&fixedvar) le 200 %then %do;  %let _fixedvarp=&fixedvar;  %end;
%else %do;  %let _fixedvarp=%substr(&fixedvar, 1, 200) and other variables ;  %end;
   

 /***** determine whether there are non-modifiable vbls.
  if so, then set partialpar=T.  ****/
%let _nfixedvar=%numargs(&fixedvar);
/* count fixed variables */
%if &_nfixedvar ne 0 %then %do;
   %let partialpar=T;
   data _fvnames_;
   length name $32 ;
   %do i=1 %to &_nfixedvar;
/* fixed vbls indexed by norder in list */
      %let _fvn_&i=%scan(&fixedvar, &i, %str( ));
      %let _fix_&&_fvn_&i = 1;
      *name="%scan(&fixedvar, &i, %str( ))";
      name="&&_fvn_&i";
      output;
      %end;
   run;
%end; 
%else %do;
   data _fvnames_;
       length name $32;  name=' ';
       output;
        run;
   %let partialpar=F;
  %end;

/* count modifiable vbls */
%let _nmodvar=%numargs(&modvar);
%if &_nmodvar eq 0 %then %do;
   %let partialpar=F;  /* no actual modifiable vbls */
   %let fullpar=T ;
   data _vvnames_;  length name $32 ;  name=' ';  run;
%end;
%else %do;
   data _vvnames_;
   length name $32 ;
   %do i=1 %to &_nmodvar;
/* modifiable vbls indexed by order in list */
      %let _vvn_&i = %scan(&modvar, &i, %str( ));
      %let _fix_&&_vvn_&i = 0;
      name="&&_vvn_&i";
      output;
      %end;
   run;
%end;

/****** check for overlap between fixedvar and modvar lists  ****/
proc sort data=_fvnames_;  by name;  run;
proc sort data=_vvnames_;  by name;  run;

data _trynames_;
   merge _fvnames_ (in=inf) _vvnames_ (in=inv);
   by name;
   retain _overlap_ 0;
   if inf and inv then do;
       _overlap_=max(1, _overlap_);
       output;
       end;
   call symput ('_overlap_', _overlap_);
   run;
 
%if &_overlap_ eq 1 %then %do;
   data _null_;
   put 'ERROR: in PAR macro call:  The FIXEDVAR list overlaps with the MODVAR list.';
   put 'The PAR macro will stop.';

   file print;
   put 'ERROR: in PAR macro call:  The FIXEDVAR list overlaps with the MODVAR list.';
   put 'The PAR macro will stop.';
   run;

   ods listing;
   options notes;
   proc print data=_trynames_;
   title5 'Variable names that are in both FIXEDVAR and MODVAR';
   run;
   title5;
   
options &notes;
%end;


%if &_overlap_ eq 1 %then %goto errend;

/***** check that there is something to do ******/   
%if &_nmodvar eq 0 and &_nfixedvar eq 0 %then %do;
   data _null_;
   put 'ERROR: in PAR macro call:  You did not give any values for MODVAR or FIXEDVAR';
   put 'The PAR macro will stop.';
   file print;
   put 'ERROR: in PAR macro call:  You did not give any values for MODVAR or FIXEDVAR';
   put 'The PAR macro will stop.';
   run;
%end;
%if &_nmodvar eq 0 and &_nfixedvar eq 0 %then %goto errend;
   
/******* modify the parameters depending on FIXEDVAR and MODVAR ******/
%if partialpar eq T %then %do;
   %if &_nmodvar eq 0 %then %do;
      %let fullpar=T;
      data _null_;
      put 'WARNING for PAR macro run:  You did not specify MODVAR.';
      put '   Only full PAR will be computed.';
      file print;
      put 'WARNING for PAR macro run:  You did not specify MODVAR.';
      put '   Only full PAR will be computed.';
      run;
      %end;
   %if &_nfixedvar eq 0 %then %do;
      %let fullpar=T;
      data _null_;
      put 'WARNING for PAR macro run:  You did not specify FIXEDVAR.';
      put '   Only full PAR will be computed.';
      file print;
      put 'WARNING for PAR macro run:  You did not specify FIXEDVAR.';
      put '   Only full PAR will be computed.';
      run;
      %end;
%end;

/******** make sure have a prevalence v-c matrix if using SURVEY ******/
%if &vprevopt eq SURVEY and %length(&vpdata) eq 0 %then %do;
  data _null_;
  put 'ERROR: in PAR macro call:  You specified VPREVOPT=SURVEY, but did not give a value for VPDATA.';
  put 'The PAR macro will stop.';
  file print;
  put 'ERROR: in PAR macro call:  You specified VPREVOPT=SURVEY, but did not give a value for VPDATA.';
  put 'The PAR macro will stop.';
  run;
%end;
%if &vprevopt eq SURVEY and %length(&vpdata) eq 0 %then %goto errend;;

   
   
   /*=====================================================================*/
data &bdata;
    set &bdata;
    %if &_nmissvar ne 0 %then %do; drop &missvarlist;%end;;
 run;
 
/**** read in the estimates/betas *****/
data betas;
   set &bdata;
   where _type_ eq: 'PAR';
   drop  intercept _link_ _status_ _lnlike_ _name_;
   run;

data betas;
    set betas;
    drop _TYPE_;
 run;

/****** check whether any of the modifiable variables has a negative coefficient ******/
data modbeta;
 set betas (keep=&modvar);
 modblt0=0;

 %do i=1 %to &_nmodvar;
    if &&_vvn_&i lt 0 then do;
    modblt0=1;
    put "ERROR: in PAR macro call:  &&_vvn_&i has negative coefficient:" &&_vvn_&i ;
    end; 
%end;
file print;
%do i=1 %to &_nmodvar;
    if &&_vvn_&i lt 0 then do;
    modblt0=1;
    put "ERROR: in PAR macro call:  &&_vvn_&i has negative coefficient, " &&_vvn_&i ;
    end; 
%end;

call symput ('_modblt0', modblt0);
run;

%if &_modblt0 %then %goto errend;


   /***  list of variable names in v-c matrix numbered by order in model 
    also counting number of variables in model ***/
/** RL: change bvlist to covs **/
*data bvlist;
 data covs;
   set &bdata;
   where _type_ eq: 'COV';
   _name_=left(upcase(_name_));
   if _name_ ^in ('INTERCEPT', 'INTERCPT', 'INTERCEP');
/* get rid of rows corresponding to missvarlist */
%if &_nmissvar ne 0 %then %do;

/*the following are ok, just feel a little bit odd, so change the position:
%do i=1 %to &_nmissvar;  if "&&_missvar&i" eq _name_ then delete;  %end; */
  %do i=1 %to &_nmissvar;  if _name_ eq "&&_missvar&i" then delete;  %end;
%end;
  * keep _name_;
 drop Intercept _link_ _status_ _lnlike_;
run;

proc sort; by _name_;
    
data bvlist;  
 set covs(keep=_name_)  end=_end_ ;
   length name $32 ;
   name=_name_;
   call symput('varname'||trim(left(_n_)), left(name));
   if _end_ then call symput('_nbet', trim(left(_n_)));
   run;
   
proc sort data=bvlist;  by name;  run;

data _maclist_;  set _fvnames_ _vvnames_;  run;
proc sort data=_maclist_;  by name;  run;

/**
%put _user_;

proc print data=_fvnames_;
    proc print data=_vvnames_;
proc print data=_maclist_;
proc print data=bvlist;
    **/

data _chk_;  merge bvlist (in=inb) _maclist_ (in=inm);  by name;
   inbb=inb;  inmm=inm;
   label inbb='in V-C matrix'  inmm='in macro call';
   if sum(inb, inm) ne 2 and name ne ' ' then problem=1;  else problem=0;
   run;
proc means noprint data=_chk_;  var problem;
   output  out=_chk1_  sum= ;
   run;
data _chk1_;  set _chk1_;
   call symput ('_varnameprb_', problem);
   run;
%if &_varnameprb_ ne 0 %then %do;
   data _null_;
   put 'ERROR: in PAR macro call: The variable names in the macro call do not';
   put '   match those in the V-C matrix of the coefficients.';
   put 'The PAR macro will stop.';
   put '   See the output, which shows the variables that do not match.';
   file print;
   put 'ERROR: in PAR macro call: The variable names in the macro call do not';
   put '   match those in the V-C matrix of the coefficients.';
   put 'The PAR macro will stop.';
   put '   See the output, which shows the variables that do not match.';
run;
   
   title5 'Variable names that are not the same in the macro call '; 
   title5 '  and the BDATA dataset';
   proc print data=_chk_ label;  var name inbb inmm;  where problem eq 1;
   run;
   title5;
%end;
 %if &_varnameprb_ ne 0 %then %goto errend;
/*
data _null_;
   nf=&_nfixedvar;  nv=&_nmodvar;  nb=&_nbet;
   nt=sum(nf, nv);
   if nb ne nt then incomp=1;  else incomp=0;
   if incomp eq 1 then do;
      put 'ERROR: in PAR macro run:  The number of variables specified in';
      put '  the macro call, ' nt ' is not the same as';
      put '  the number of coefficients in BDATA, ' nb '.';
      file print;
      put 'ERROR: in PAR macro run:  The number of variables specified in';
      put '  the macro call, ' nt ' is not the same as';
      put '  the number of coefficients in BDATA, ' nb '.';
      end;
   call symput ('_incomp', incomp);
   run;
%if &_incomp eq 1 %then %goto errend;
*/

    
   /*===========================================================================*/
   
%if "&vprevopt" eq "COHORT" and %length(&sampsz) eq 0 %then %do;
   %let vprevopt = FIXED ;
   data _null_;
   put 'WARNING for PAR macro run:  You specified COHORT variance-covariance';
   put '     matrix, but did not give a sample size.';
   put '  The PAR macro will use VPREVOPT=FIXED.';
   file print;
   put 'WARNING for PAR macro run:  You specified COHORT variance-covariance';
   put '     matrix, but did not give a sample size.';
   put '  The PAR macro will use VPREVOPT=FIXED.';
   run;
%end;



   /**********************now beginning real computations ******************/   
   /*=====================================================================*/
   
   /* make macro variables of betas -- macro name for variable X is b_X */
data _betas;  
   set betas;
   array betas{*} _betas1 - _betas&_nbet ;
%do i=1 %to &_nbet;
   *call symput('b_'||trim(left("&&varname&i")), trim(left(&&varname&i)));
   call symput('b_'||trim(left("&&varname&i")), &&varname&i);
   betas{&i} = &&varname&i;
%end;
   run;
 
   
   /*==============================================================================*/
/******* now get the v-c matrix ****/   
data covbs;  
   /*set &bdata;  where _type_ eq: 'COV'; */
   set covs;
   _name_=upcase(_name_);
/**
   if _name_ in ('INTERCEPT', 'INTERCPT', 'INTERCEP') then delete;
   if index("&missvarlist", _name_) gt 0 then delete;
**/
%if &_nmissvar ne 0 %then %do;
  %do i=1 %to &_nmissvar;
    if _name_ eq "&&_missvar&i" then delete;
  %end;
%end;
   keep _name_ &fixedvar &modvar  ;  /* keeps only columns from v-c matrix plus _name_ */
   run;

 
   
   /* make macro variables of variance-covariance matrix of betas */
   /* this is now set up so that the dynamic vbl _n_ counts the rows of the v-c matrix */
   /* the indexing of cvb is the same as for variable names */
data _covbs;  set covbs;
   array covs{*} _covbs1 - _covbs&_nbet;
%do i=1 %to &_nbet;
  
  call symput('cvb_'||trim(left(&i))||'_'||trim(left(_n_)), trim(left(&&varname&i)));

/* indexed as ith column, _n_th row */
   covs{&i}=&&varname&i ;
%end;
   run; 

   /*==========================================================================*/
   
   /*  fix up to get prevs  */
   /* regardless of value of n_or_p, end up with a vbl _prev_ that has the 
    prevalence in the data set &pdata */
/* if no missing indicators, pct of data used in par computation is 100 */
%if &_nmissvar eq 0 %then %do;
  %let _pctused = 100 ;
  %if &n_or_p eq N %then %do;
    proc means noprint data=&pdata;  var &n_or_pname;
    output  out=_totn_  sum=totn;
    run;
    data _pdat;  if _n_ eq 1 then set _totn_;
    set &pdata  end=_end_;
    _prev_=&n_or_pname/totn;

 /** RL added the following two lines to get _nlev, 6/1/2017; **/  
   if _end_ then call symput('_nlev', trim(left(_n_)));  
    run;

   %end;
  %else %do;  /* n_or_p=P */
    data _pdat;  set &pdata  end=_end_;  _prev_=&n_or_pname;

   ** RL added the following two lines to get _nlev, 6/1/2017;  
   if _end_ then call symput('_nlev', trim(left(_n_)));
 %end;
%end;  /* end of nmissvar=0 */

%else %do;  /* there are missing indicators */
%if "&n_or_p" eq "N" %then %do;
   /* get original total N */
   proc means noprint data=&pdata;  var &n_or_pname;
   output  out=_tn_  sum=_totn_;
   run;
  /*  14mar2016:  code to delete all lines with a missing indicator=1 */
   data newpdata;  set &pdata;
   %do i=1 %to &_nmissvar;  if &&_missvar&i eq 1 then delete;  %end;
   run;
   proc means noprint data=newpdata;  var &n_or_pname;
   output  out=_newtn_  sum=_newtotn_ ;
   run;
   data _null_;  if _n_ eq 1 then set _tn_;  set _newtn_;
   pctused=round(100*_newtotn_/_totn_, .1);
   call symput ('_pctused', trim(left(pctused)));
   run;
   data _pdat;
   if _n_ eq 1 then set _newtn_;
   set newpdata  end=_end_;
   _prev_= &n_or_pname / _newtotn_;
   call symput('_prev_'||trim(left(_n_)), trim(left(_prev_)));
   /* _prev_&n= prev of nth level in prevalence file */
   if _end_ then call symput('_nlev', trim(left(_n_)));
   run;
%end;
%else %do;  /*  prevalences */
   data _pdatnew;  set &pdata;
   %do i=1 %to &_nmissvar;  if &&_missvar&i eq 1 then delete;  %end; 
   run;
   proc means noprint data=_pdatnew;  var &n_or_pname;
   output  out=_newtp  sum=_newtotp_;
   run;
   
   data _newtp;  set _newtp;
       pctused=round(100*_newtotp_, .1);
       call symput ('_pctused', trim(left(pctused)));
   run;

   data _pdat;  if _n_ eq 1 then set _newtp;
        set _pdatnew  end=_end_;
  _prev_=&n_or_pname/_newtotp_;
   call symput('_prev_'||trim(left(_n_)), trim(left(_prev_)));
   /* _prev_&n= prev of nth level in prevalence file */
   if _end_ then call symput('_nlev', trim(left(_n_)));
   run;
%end;  /* end of n_or_p=P */
%end;  /*  end of else do--missvarlist not empty **/



options notes;
data _null_;
put 'After deleting observations with missing data,';
put "the percent of the data used in the PAR computation is &_pctused..";
file print;
put 'After deleting observations with missing data,';
put "the percent of the data used in the PAR computation is &_pctused..";
run;

options &notes;


data _pvec;  set _pdat (keep=_prev_);  run;
proc transpose data=_pvec  out=_pvect prefix=_pvect;  var _prev_;  run;
   
/*  may need to modify for missing indicators */
%if &vprevopt eq SURVEY %then %do;  data _pv_;  set &vpdata;
%do i=1 %to &_nlev;  _pv_&i =&vpvarname.&i;  %end;
   run;  %end;
  
   
   /*===========================================================================*/

   /*===========================================================================*/
   /******** this is the part where we get the sums and the partial deriviatives *******/
data _freqs_;
   retain _rrprob_ 0;
   if _n_ eq 1 then set _betas;
   set _pdat  end=_end_ ;
   _logrr_=0;  _fixlogrr_=0; 
   array _betas{*} _betas1 - _betas&_nbet;
%do q=1 %to &_nbet;  /* within a stratum, add up all the x(i)*beta(i) */
   _summand_=&&varname&q * _betas{&q};
   _logrr_=sum(_logrr_, _summand_);  /* logrr for all vbls */
   /* summing the fixed vbls and the modifiable vbls separately, if nec */
   /*  for partial par, but compute for everything */
   _fixlogrr_=sum(_fixlogrr_, &&&&_fix_&&varname&q * _summand_);
   /* logrr for the fixed (out) variables */
   %end;
   _rr_=exp(_logrr_);   label _rr_ =' rr for stratum' ;
   if _rr_ lt 1 then _rrprob_=1;
   _fixrr_=exp(_fixlogrr_);  label _fixrr_=' rr for fixed vbls for stratum';
   if _fixrr_ lt 1 then _rrprob_=1;
   _prr_=_prev_ * _rr_ ;  label _prr_='prevalence * rr';
   _pfrr_=_prev_ * _fixrr_ ;  label _pfrr_='prev*rr for fixed vbls';
   _sumpar_=_prev_ * (_rr_ -1);  label _sumpar_='summand--prev*(rr-1)';
* put _n_ _prev_= _rr_= _sumpar_= _fixrr_= _pfrr_= _rrprob_=;
   
%do q=1 %to &_nbet;
   sdpfbeta&q=_prr_ * &&varname&q ;  /* summand for numerator of partial deriv of 
				      full par by b&&varname&q */
   snum1dppx&q=_pfrr_*&&varname&q * &&&&_fix_&&varname&q ; /* partial first factor of first part for
							    partial deriv of partial par by b&&varname&q */			      
   snum2dppx&q = _prr_ * &&varname&q;  /* partial first factor of second part of
					numerator of partial deriv of partial
					par by b&&varname&q */
    /*  note:  this is the same as sdpfbeta&q */
%end;
   
   /* for the matrices needed to compute the variance of full par */
   call symput('_rr_'||trim(left(_n_)), _rr_);
   /*  rr&n = rr of nth combination in prevalence file */
   call symput('_fixrr_'||trim(left(_n_)), _pfrr_);
   /* _fixrr_&n = rr for fixed vbls in nth combination in prevalence file */
   call symput('_numdpfprev_'||trim(left(_n_)), trim(left(_rr_ - 1)));
   /* numerator for partial deriv of full par by prevalence _n_ */
   
   call symput ('_rrprob_', trim(left(_rrprob_)));
   run;
   
  /* some _rr_ or _fixrr_ is not ge 1  */

%if &_rrprob_ eq 1 %then %do;
   data _null_;
   put 'WARNING for PAR macro run:  Not all of the RRs for the combinations of';
   put '   all or fixed variables are at least 1.';
   put '   The PAR macro will continue, but please note that';
   put '   the computed full PAR is not the fraction of cases prevented';
   put '   if all the risk factors are eliminated, but the';
   put '   fraction of the risk averted if all';
   put '   risk factors in the model have their reference levels/values';
   file print;
   put 'WARNING for PAR macro run:  Not all of the RRs for the combinations of';
   put '   all or fixed variables are at least 1.';
   put '   The PAR macro will continue, but please note that';
   put '   the computed full PAR is not the fraction of cases prevented';
   put '   if all the risk factors are eliminated, but the';
   put '   fraction of the risk averted if all';
   put '   risk factors in the model have their reference levels/values';
   file print;
   run;
   
   proc print data=_freqs_ (obs=10);  var &modvar &fixedvar _rr_ _fixrr_;
   where _rr_ lt 1 or _fixrr_ lt 1;
   title5 'variable combinations with _rr_ or _fixrr_ < 1';
   run;
   title5;
%end;
/* %if &_rrprob_ eq 1 %then %goto errend;  */

   /*===========================================================================*/
   /* to get sum of pi*rri, etc., over levels of covariates */
proc means noprint data=_freqs_;
   var _prr_ _pfrr_ _sumpar_ 
   %do q=1 %to &_nbet ;  sdpfbeta&q   snum1dppx&q snum2dppx&q  %end;;
   output  out=_sumprr_  sum=_sumprr_ _sumpfrr_  _ssumpar_  
      %do q=1 %to &_nbet;  _numpfbeta&q   _num1ppbeta&q _num2ppbeta&q   %end;
   ;
   run;


  
   
data _sumprr_;  set _sumprr_;
   label _sumprr_='sum(prev*rr)'
      _sumpfrr_='sum(prev*rr for fixed vbls)'
      _ssumpar_='sum(prev*(rr-1))'
      ;
   _sumprr2_=_sumprr_ * _sumprr_;  /* denominator for partial derivatives for partial par */
   _oneplssumpar2_=(1 + _ssumpar_)**2;  /*  denominator for partial derivatives 
					 for full par  */

   
%if "&fullpar" eq "T" %then %do;
   /******  now compute all partial derivatives for full par  **********/
 /* by each prevalence level */
   array _pdfprev{*} _pdfprev1 - _pdfprev&_nlev;
%do i=1 %to &_nlev;
   _pdfprev&i=&&_numdpfprev_&i / _oneplssumpar2_ ;
* put "&i " _pdfprev&i;
   call symput("_pdfprev&i", _pdfprev&i);
%end;
    /* by each beta */
   array pdfbeta{*} _pdfbeta1 - _pdfbeta&_nbet;
%do q=1 %to &_nbet;
   pdfbeta{&q} = _numpfbeta&q / _oneplssumpar2_ ;
   call symput("_pdfbeta&q", _pdfbeta&q);
%end;
%end;


%if "&partialpar" eq "T" %then %do;
   /***********  now compute all partial derivatives for partial par  ************/
   /* by each prevalence level */
   array pdpprev{*} _pdpprev1 - _pdpprev&_nlev;
%do i=1 %to &_nlev;
   numpdprev&i = &&_fixrr_&i * _sumprr_ - &&_rr_&i * _sumpfrr_ ;
   _pdpprev&i= numpdprev&i / _sumprr2_;
* put "&i " _pdpprev&i;
   call symput("_pdpprev&i", trim(left(_pdpprev&i)));
   %end;
   /* by each beta */
   array pdpbeta{*} _pdpbeta1 - _pdpbeta&_nbet;
%do q=1 %to &_nbet;
   num1&q = _num1ppbeta&q  * _sumprr_;
   num2&q = _num2ppbeta&q * _sumpfrr_ ;
   pdpbeta{&q} = (num1&q - num2&q)/ _sumprr2_ ;
   call symput("_pdpbeta&q", trim(left(_pdpbeta&q)));
%end;
%end;
   
   _fullpar_=_ssumpar_/(1+_ssumpar_);  label _fullpar_='full PAR';
   _partialpar_=(_sumprr_ - _sumpfrr_)/_sumprr_ ;  label _partialpar_='partial PAR';
   call symput ('_fullpar_', trim(left(_fullpar_)));
   call symput ('_partlpar_', trim(left(_partialpar_)));  /* note:  macro vbl must
							   not have same name as parameter passed to macro */
   run;
   
  /* vector datasets for vprimematv */
%if &fullpar eq T %then %do;
data _pdfbeta;  set _sumprr_(keep=_pdfbeta1 - _pdfbeta&_nbet);
data _pdfprev;  set _sumprr_(keep=_pdfprev1 - _pdfprev&_nlev);
%end;
%if &partialpar eq T %then %do;
data _pdpbeta;  set _sumprr_(keep=_pdpbeta1 - _pdpbeta&_nbet);
data _pdpprev;  set _sumprr_(keep=_pdpprev1 - _pdpprev&_nlev);
%end;
   
  /*******  getting variances of full and partial par  *********/
%if "&fullpar" eq "T" %then %do;
  /*  full  */
%if &VPREVOPT eq SURVEY %then %do;
   %vprimematv(v=_pdfprev, mat=_pv_, dimv=&_nlev, xsum=_varfull1_);
%end;
%if &VPREVOPT eq COHORT %then %do;
   data _null_;
   _varfull1_=0;
   %do i=1 %to &_nlev;
      _varfull1_=sum(_varfull1_, (&&_pdfprev&i * &&_prev_&i * (1 - &&_prev_&i) * &&_pdfprev&i / &sampsz));
      %if &i gt 1 %then %do;
	 %do j = 1 %to %eval(&i -1);
	    _varfull1_ = sum(_varfull1_, &&_pdfprev&i * &&_prev_&i * &&_prev_&j * &&_pdfprev&j / &sampsz ) ;
			     %end;
			     %end;
      %end;
   call symput('_varfull1_', trim(left(_varfull1_)));
   run;
%end;
%if &VPREVOPT eq FIXED %then %do;
   %let _varfull1_ = 0 ;
%end;
   %vprimematv(v=_pdfbeta, mat=_covbs, dimv=&_nbet, xsum=_varfull2_);
   data _null_;  varfull=&_varfull1_ + &_varfull2_ ;
   call symput('_varfull_', trim(left(varfull)));
   run;
   %if &citransf eq LOG %then %do;
   data _null_;
      _uclfull_=1-(1-&_fullpar_) * exp(-1.96*sqrt(&_varfull_)/(1-&_fullpar));
      _lclfull_=1-(1-&_fullpar_) * exp(1.96*sqrt(&_varfull_)/(1-&_fullpar));
      call symput('_uclfull_'. trim(left(_uclfull_);
      call symput('_lclfull_'. trim(left(_lclfull_);
      run;
   %end;
   %else %do;
     /* fisher's z transforming and back */
      %fztransf(quant=&_fullpar_, varquant=&_varfull_, name=full);
   %end;
%end;
   
%if "&partialpar" eq "T" %then %do;
  /*  partial  */
%if &VPREVOPT eq FIXED %then %do;
   %let _varpart1_ = 0;
   %end;
%if &vprevopt eq COHORT %then %do;
   data _null_;
   _varpart1_=0;
   %do i=1 %to &_nlev;
      _varpart1_= sum(_varpart1_, (&&_pdpprev&i * &&_prev_&i * (1 - &&_prev_&i) * &&_pdpprev&i / &sampsz));
      %if &i gt 1 %then %do;
	 %do j = 1 %to %eval(&i -1);
	    _varpart1_ = sum(_varpart1_ , &&_pdpprev&i * &&_prev_&i * &&_prev_&j * &&_pdpprev&j / &sampsz ) ;
	    %end;
	 %end;
      %end;
   call symput('_varpart1_',trim(left(_varpart1_)));
   run;
%end;
%if &VPREVOPT eq SURVEY %then %do;
   %vprimematv(v=_pdpprev, mat=_pv_, dimv=&_nlev, xsum=_varpart1_);
%end;
   %vprimematv(v=_pdpbeta, mat=_covbs, dimv=&_nbet, xsum=_varpart2_);
   data _null_;  varpart=&_varpart1_ + &_varpart2_ ;
   call symput('_varpart_', trim(left(varpart)));
   run;
   %if &citransf eq LOG %then %do;
   data _null_;
      _uclpartial_=1-(1-&_partlpar_) * exp(-1.96*sqrt(&_varpart_)/(1-&_partlpar));
      _lclpartial_=1-(1-&_partlpar_) * exp(1.96*sqrt(&_varpart_)/(1-&_partlpar));
      call symput('_uclpartial_'. trim(left(_uclpartial_);
      call symput('_lclpartial_'. trim(left(_lclpartial_);
      run;
   %end;
   %else %do;
       /* fisher's z transforming and back */
      %fztransf(quant=&_partlpar_, varquant=&_varpart_, name=partial);
   %end;
%end;
   
  /*********  finally, printing out the results  ************/

title "&modvar AND &fixedvar";
data _null_;
   file print;
   put ' ';
   put "option for the variance-covariance matrix of the prevalences is &vprevopt .";
%if "&fullpar" eq "T" %then %do;
   fullpar=round(&_fullpar_, .001);
   uclfull=round(&_uclfull_, .001);
   lclfull=round(&_lclfull_, .001);
   put ' ';
   put "Full PAR (95% CI) for";
   put "  &modvar";
   put "  &_fixedvarp";
   put "           "  fullpar "(" lclfull " , " uclfull ")";
%end;
%if "&partialpar" eq "T" %then %do;
   partlpar=round(&_partlpar_, .001);
   uclpartial=round(&_uclpartial_, .001);
   lclpartial=round(&_lclpartial_, .001);
   put ' ';
   put "Partial PAR (95% CI) for";
   put "  modifiable vbls :  &modvar ";
   put "  fixed vbls :  &_fixedvarp ";
   put "           "  partlpar "(" lclpartial " , " uclpartial ")";
%end;
   run;
   
   /*===========================================================================*/
   
   options &notes;
proc datasets nolist;  delete betas covbs _tn_ bvlist _pdat _freqs_
   _betas _covbs _pvect _pv_ _pdpprev _pdpbeta _pdfpref _pdfbeta
   ;
   run;
   
%goto realend;
   
%errend:
   quit;
   
%realend:  
   quit;
%let xsum=.;
%let _lclfull_=.;  %let _lclpartial_=.;
%let _uclfull_=.;  %let _lclpartial_=.;
proc optload  data=_tmpoptsfromrun_;  run;

     
%mend;


  /***************************************************/
   
%macro vprimematv(v=, mat=, dimv=, xsum=);
/* does x` <matrix> x multiplication */
   options nonotes;
   /*  for pre- and post-multiplying a symmetric matrix by a vector and its transpose.
    assumes order of subscripts is ok and that everything is macro variables
    of the specified forms  */
  /* v is both the name of the data set and the name of the variable
   mat is both the name of the data set and the name of the variable */
%global &xsum;
data _mm1_;
   if _n_ eq 1 then set &v;
   set &mat;
   array vs{*} &v.1 - &v&dimv;
   array mats{*}  &mat.1 - &mat&dimv;
   rowsum=vs{_n_} * mats{_n_} * vs{_n_};  /* the diagonal term */
   if _n_ gt 1 then do j=1 to (_n_ - 1);
      rowsum=sum(rowsum, 2 * vs{j} * mats{_n_} * vs{_n_});
      end;
 /*  old version
  rowsum=0;
  do j=1 to &dimv;
  rowsum=sum(rowsum, vs{_n_} * mats{j} * vs{j} );
  */
   run;
   
proc means noprint data=_mm1_;  var rowsum;
   output  out=_mm2_  sum=mmsum;
   run;
   
data _mm2_;
   set _mm2_;
   call symput("&xsum", mmsum);
   run;
proc datasets nolist;  delete _mm1_ _mm2_;  run;
   options notes obs=max;
%mend;
/***************************************************/


   
%*                                                                            ;
%* This macro was developed by Carrie Wager,Programmer,Channing Laboratory 1990;
%* Modified AMcD 1993, e hertzmark 1994 and L Chen 1996;
%*                                                                            ;
   
%macro numargs(arg, delimit);
   %if %quote(&arg)= %then %do;
      0
	 %end;
   %else %do;
      %let n=1;
      %do %until (%qscan(%quote(&arg), %eval(&n), %str( ))=%str()); 
	 %let n=%eval(&n+1);
	 %end;
      %eval(&n-1)
	 %end;
%mend numargs;

%macro fztransf(quant=, varquant=, name=quant);
/*
%put variable is &name;
%put quantity is &quant;
%put variance of quantity is &varquant;
*/
%global _lcl&name._ _ucl&name._;
data _null_;
   varquant=&varquant;  quant=&quant;  /* making regular vbls */
   deriv=.5/(1+quant) + .5/(1-quant);
   varztq=varquant * deriv**2;
   seztq=sqrt(varztq);  /* std error */
   ztq=.5*log((1+quant)/(1-quant));  /* z-transform of quant */
   lclztq=ztq-1.96*seztq;  uclztq=ztq+1.96*seztq;  /* confidence bounds of transformed quant */
   lclquant=(exp(2*lclztq)-1) / (exp(2*lclztq)+1);  /* invert fisher's z transform to get orig scale*/
   uclquant=(exp(2*uclztq)-1) / (exp(2*uclztq)+1);  /* invert fisher's z transform to get orig scale*/
   
   call symput("_lcl&name._", trim(left(lclquant)));
   call symput("_ucl&name._", trim(left(uclquant)));
   run;
   options notes obs=max;
%mend;



/***
%put _user_;


%errend:
   quit;
   
%realend:  
   quit;

%mend;

%macro notyet();
***/    
