*--------------------------------------------------------------;
* Ratio estimator of population mean or total from a two-stage ;
* cluster sample.                                              ;
*--------------------------------------------------------------;
%macro est_rat(sample=,setup=,npop=,n=,mpopi=,cluster=,ybar=,
  var=,phat=,counts=,mi=,mbar=,param=);

 %if %length(&cluster)= 0 %then %let cluster = %str(cluster);
 %if %length(&npop)   = 0 %then %let npop    = %str(npop);
 %if %length(&mpopi)  = 0 %then %let mpopi   = %str(mpopi);
 %if %length(&mi)     = 0 %then %let mi      = %str(mi);

 %if %length(&mbar)   = 0 %then %let mbar    = %str(s_mi_/n_);
 %if %length(&setup) > 0 %then %do;
     proc sort data = &setup;  by &cluster;
 %end;
 %if %length(&sample) > 0 %then %do;
     proc sort data = &sample;  by &cluster;
   %if %length(&setup) > 0 %then %do;
     data est_;  merge &setup(in=inset_) &sample(in=insamp_);  by &cluster;
        if insamp_ = 1;
   %end;
   %else %do;
     data est_;  set &sample;
   %end;
 %end;
 %else %do;
    data est_;  set &setup;
 %end;

 %if %length(&ybar) > 0 %then %do;
        ybar_ = &ybar;
         var_ = &var;
 %end;
 %if %length(&phat) > 0 %then %do;
        ybar_ = &phat;
         var_ = (&mi/(&mi-1))*ybar_*(1-ybar_);
 %end;
 %if %length(&counts) > 0 %then %do;
        ybar_ = &counts/&mi;
         var_ = (&mi/(&mi-1))*ybar_*(1-ybar_);
 %end;
      mi_ybi_ = &mpopi*ybar_;                        
      mi2_ybi_ = &mpopi**2*ybar_;                    
      fpci_ = (&mpopi-&mi)/&mpopi;                      
      term_ = &mpopi**2*fpci_*(var_)/&mi;                 

    proc means data = est_ noprint;                
      var mi_ybi_ &mpopi mi2_ybi_ term_;               
      output out = new_   uss = ssq_my_  ssq_mi_       
                          sum = sum_my_  s_mi_  sum_m2y_  s_term_
                          n = n_;

    data est_;  set new_;   
     mm = &mbar;
     mu_hat_ = sum_my_/s_mi_;                     
     sr2_ = (ssq_my_ - 2*mu_hat_*sum_m2y_ + mu_hat_**2*ssq_mi_)/
         (n_-1);                                           
     term1_ = (&npop-n_)/&npop*sr2_/n_/(&mbar)**2;                 
     term2_ = s_term_/n_/&npop/(&mbar)**2;                     
     var_mu_ = term1_ + term2_;                              
     std_mu_ = sqrt(var_mu_);                               
     bnd_mu_ = 2*std_mu_;                                    
     drop term1_ term2_;                                    
 %if %index(%upcase(&param),TOTAL) > 0 %then 
   %put %str( *** TOTAL is not an option with this macro *** );

 %if %index(%upcase(&param),MEAN) > 0 %then %do;
   proc print data = est_ noobs split='*';
     title1 'Ratio Estimate of Population Mean';
     title2 'Two-stage Cluster Design';      
     label mu_hat_ = 'Estimate';
     label std_mu_ = 'Standard*Error';
     label bnd_mu_ = 'Bound';
     label sr2_ = 's(r)^2';
     var mu_hat_ std_mu_ bnd_mu_ sr2_;
 %end;

 %if %index(%upcase(&param),PROP) > 0 %then %do;
   proc print data = est_ noobs split='*';
     title1 'Ratio Estimate of Population Proportion';
     title2 'Two-stage Cluster Design';      
     label mu_hat_ = 'Estimate';
     label std_mu_ = 'Standard*Error';
     label bnd_mu_ = 'Bound';
     var mu_hat_ std_mu_ bnd_mu_;
 %end;

run; title;
%mend est_rat;