* vcbd.sas
* xref:
* input:
* output:
* does: A macro to calculate sums of squares for a nested
*       study of subjects within providers within practices, and to
*       calculate theta, lambda1 and lambda2, by equating expected
*       and observed sums of squares.
*
*      - final version of the VCBD macro (old name s010.sas)
*************************************************************************;
%macro vcbd(data=_last_,
               i=,
               mij=,
               yij=);
***************************************************;
%let err=0;
%if %length(&data)=0 %then %do;
   %put E R R O R : The input data must be specified in the macro call.;
   %let err=1;
%end;
%if %length(&i)=0 %then %do;
   %put E R R O R : i= must be specified in the macro call.;
   %let err=1;
%end;
%if %length(&mij)=0 %then %do;
   %put E R R O R : mij = must be specified in the macro call.;
   %let err=1;
%end;
%if %length(&yij)=0 %then %do;
   %put E R R O R : yij = must be specified in the macro call.;
   %let err=1;
%end;

%if &err=1 %then %goto STOPIT;
***************************************************;

Data _A_ (keep = &i &yij &mij mij2 aij);
  set &data;
  mij2 = &mij * &mij;
  aij  = &yij * &yij / &mij;
  run;
***************************************************;
proc summary data=_A_ nway;
  class &i;
  var              &yij &mij mij2 aij;
  output out=_B_ sum=yi  mi  qi   ai   n=ni;
  run;
***************************************************;
data _C_;
  set _B_;
  bi = yi * yi / mi;
  ri = qi / mi;
  hi = mi * mi;
  run;
***************************************************;
proc summary data=_C_;
  var                yi  mi qi ai bi ri ni hi;
  output out=_D_ sum=y   m  q  a  b  r  n  h   n=k;
  run;
***************************************************;
data _E_;
  set _D_;
  ss1 = y - a;
  ss2 = a - b;
  ss3 = b - y * y / m;
  x_ = ss1/(m-n);
  y_ = (ss2 - (n-m-k+r)*x_)/(m-r);
  z_ = (ss3 - (k-1-r+q/m)*x_ - (r-q/m-m+h/m)*y_)/(m-h/m);
  ybar = y/m;
  if ((4*z_) > 1) then th = ybar;
  else do;
    th1 = (1 + sqrt(1-4*z_)) / 2;
    th2 = (1 - sqrt(1-4*z_)) / 2;
    if (abs(ybar-th1) < abs(ybar-th2)) then th = th1; *theta;
    else th = th2;      * theta is the closest to ybar;
  end;
  la1 = 1 - y_ / z_; * lambda1;
  la2 = 1 - x_ / y_; * lambda2;
  if (la1 < 0) then la1 = 0; else if (la1 > 1) then la1 = 1;
  if (la2 < 0) then la2 = 0; else if (la2 > 1) then la2 = 1;
  comp1 = (1-la1)*(1-la2) * 100;
  comp2 = (1-la1)*la2 * 100;
  comp3 = la1 * 100;
  corr1 = (comp2 + comp3) / 100;
  corr2 = (comp3) / 100;

  u12 = th*th + corr1*th*(1-th);
  or1 = u12*(1-th-th+u12)/(th-u12)**2;
  u12 = th*th + corr2*th*(1-th);
  or2 = u12*(1-th-th+u12)/(th-u12)**2;

  put "Variance Components for Nested Binary Responses (Version 1.0)" ;
  put "(c) Bahjat Qaqish and Habib Moalem (1993)" ;
  put "Department of Biostatistics";
  put "CB 7400";
  put "Chapel Hill,  NC 27599-7420";
  put ;

  put "Data set: &data" ;
  put "Response: &yij / &mij"  ;
  put "Practice ID: &i" ;
  put ;

  put "Total number of practices:  "  k ;
  put "Total number of physicians: "  n ;
  put "Total number of patients:   "  m ;
  put "Sum of &yij:  "                 y ;
  put "Mean of &yij: "                 ybar;
  put ;

  put "Estimates:";
  put "    Theta:    " th;
  put "    Lambda_1: " la1;
  put "    Lambda_2: " la2;
  put ;

  put "Variance components:";
  put "    Between practices:                 " comp3 "%";
  put "    Between providers within practice: " comp2 "%";
  put "    Between patients Within provider:  " comp1 "%";
  put ;

  put "Correlation:";
  put "    Between providers within practice: " corr2;
  put "    Between patients within provider:  " corr1;
  put ;

  put "Odds ratio:";
  put "    Between providers within practice: " or2;
  put "    Between patients within provider:  " or1;
  put "  " / "  " / "Macro VCBD done.";
  run;
***************************************************;
  %STOPIT:
%mend vcbd;