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

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

 %if %length(&mbar) = 0 %then %do;
    %let mbar = %str(s_mi_/n_);
 %end;
 %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_;                         
       fpci_ = (&mpopi-&mi)/&mpopi;                       
       term_ = &mpopi**2*fpci_*(var_)/&mi;
     proc means data = est_ noprint;                 
       var mi_ybi_ &mpopi term_;  
       output out = new_ sum = s_miyi_ s_mi_ s_term_  
            uss = sq_miyi_ n = n_; 

     data est_;  set new_;
       mtot_ = (&mbar)*&npop;                
       mu_unb_ = &npop*s_miyi_/(mtot_*n_); 
       sb2_ = (sq_miyi_ - 2*(&mbar)*mu_unb_*s_miyi_+n_*((&mbar)*mu_unb_)**2)  
              /(n_-1);                      
       part1_ = (&npop-n_)/&npop*sb2_/n_/(&mbar)**2;
       part2_ = s_term_/n_/&npop/(&mbar)**2; 
       var2_ = part1_ + part2_;        
       std_mu_ = sqrt(var2_);            
       bnd_mu_ = 2*std_mu_;   
       tau_hat_ = &npop*(&mbar)*mu_unb_; 
       std_tau_ = &npop*(&mbar)*std_mu_;         
       bnd_tau_ = &npop*(&mbar)*bnd_mu_;       

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

 %if %index(%upcase(&param),TOTAL) > 0 %then %do;
   proc print data = est_ noobs split='*';  
      title1 "Unbiased Estimate of Population Total";
      title2 "Two-stage Cluster Design";      
      label tau_hat_ = 'Estimate';
      label std_tau_ = 'Standard*Error';
      label bnd_tau_ = 'Bound';
      label sb2_ = 's(b)^2';
      var tau_hat_ std_tau_ bnd_tau_ sb2_;
 %end;

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

 run; title;
%mend est_unb;