%macro irreml(data=,
              fix=,
              rand=,
              class=,
              cont=,
              score=,
              lsmeans=,
              lsm_opt=,
              noparms=,
              iter=50,
              n_ij=,
              cutpoint=,
              fixres=yes,
              scoring=1,
              converge=1e-8,
              link=probit);

/***************************************************************************************/
/* This macro implements the IRREML algorithm of Keen A, Engel B 1996 Analysis of a    */
/* mixed model for ordinal data by iterative re-weighted REML.                         */
/* Statistica Neerlandica (in press)                                                   */
/*  data=      sas-dataset containing data; must be sorted by classes and ordered      */
/*             categories within classes; all categories must be included, even if     */
/*             cell count is 0                                                         */
/*  fix=       MODEL statement of MIXED                                                */
/*  rand=      RANDOM statement of MIXED                                               */
/*  class=     CLASS statement of MIXED                                                */
/*  cont=      non-classifier independent variables                                    */
/*  score=     name of variable holding score                                          */
/*  lsmeans=   LSMEANS statement of MIXED                                              */
/*  lsm_opt=   Options fo LSMEANS statement                                            */
/*  noparms=   How many variance components, including residual?                       */
/*  iter=      How many iterations?  Default is 50                                     */
/*  n_ij=      Variable that contains frequencies of scores within class               */
/*  cutpoint=  Number of cutpoints (=no of ordered categories -1)                      */
/*  fixres=    Do you want to fix the residual variance parameter (at value 1 for      */
/*             probit link, at value (pi**2)/3 for logit link)? Options:               */
/*             (fixres=yes) or not (fixres=no). In the latter case, an extra (over)    */
/*             dispersion parameter is fitted for the residual term                    */
/*  mix_iter   Number of iterations per call of MIXED to compute variance components.  */
/*             Fisher's scoring is used. Default is mix_iter=1.                        */
/*  converge   Convergence criterion; the iterations are terminated when               */
/*             sum of squares of successive cutpoint estimates < converge              */
/*             The default is converge=1e-8                                            */
/*  link       Link function for (conditional) predictor; options: probit and logit    */
/*             default is probit.                                                      */
/*                                                                                     */
/*  During iterations variance components < 0.0001 are set equal to 0.0001             */
/*  Macro requires SAS 6.12 modules BASE, STAT and IML                                 */
/*  Developed under MS Windows for Workgroups  3.11                                    */
/*  Macro written by Hans-Peter Piepho, University of Kassel, Germany                  */
/*  (piepho@wiz.uni-kassel.de)                                                         */
/*  Last update 11 Aug 1997                                                            */
/***************************************************************************************/



*proc printto print=NUL;
*proc printto log=NUL;


proc glmmod outdesign=des1 data=&data noprint;
class &class;
weight &n_ij;
model &score=&fix/noint;
run;

proc transpose data=des1 out=a;
data;
i=1;
set a NOBS=dim point=i;
dim2=dim-2;
call symput('dim_x',dim2);
output;
stop;

proc glmmod outdesign=des2 data=&data noprint;
class &class;
model &score=&rand/noint;
run;

proc transpose data=des2 out=a;
data;
%let dim_z=0;
i=1;
set a NOBS=dim point=i;
dim2=dim-1;
call symput('dim_z',dim2);
output;
stop;

data design;
set &data;             /*corrected May 28, 1997. Error detected by Dirk Seidel*/
keep &class &cont;

data solf;
do i=1 to (&cutpoint+&dim_x);
  estimate=0.1;
  output;
end;

data solr;
do i=1 to &dim_z;
  estimate=0;
  output;
end;

data dat_thet;
do i=1 to &cutpoint;
  theta=(i-1)/10;
  output;
end;

%if %UPCASE(&link)=PROBIT %then %let link=0;
%if %UPCASE(&link)=LOGIT  %then %let link=1;


data cov;
do i=1 to %eval(&noparms-1);
  estimate=0.1;
  output;
end;
estimate=1; output;

/*Daten mussen nach i, j sortiert sein!
  Fuer jedes i muessen alle j vorkommen!*/


%if %UPCASE(&fixres)=YES %then %let eqcons=&noparms;
%if %UPCASE(&fixres)=NO  %then %let eqcons=0;

%let x_last=%str("x&cutpoint");

%if &dim_z>0 %then %let _random_=%str(random &rand/solution;);
%if &dim_z>0 %then %let _make_R_=%str(make 'solutionR' out=solr;);

%if &dim_z=0 %then %let _random_=%str();
%if &dim_z=0 %then %let _make_R_=%str();

%let mix_iter=%eval(&scoring+1);



%DO it=1 %to &iter;

data old_thet;
set dat_thet;
old_thet=theta;
keep old_thet;

proc iml;

bool_z=0+&dim_z;

use des1;
read all var _num_ into x;
close des1;

if bool_z>0 then do;
  use des2;
  read all var _num_ into z;
  close des2;

  colz=ncol(z);
  z=z[,2:colz];

  use solr;
  read all var {estimate} into u;
  close solr;

end;

use solf;
read all var {estimate} into beta;
close solf;

use dat_thet;
read all var {theta} into theta;
close dat_thet;

start phi(c,d);
  if &link=0 then do;
    ppii=2*arsin(1);
    c=1/sqrt(2*ppii)*exp(-d*d/2);
  end;
  if &link=1 then c=exp(-d)/((1+exp(-d))**2);
finish phi;

start cdf(c,d);
  if &link=0 then c=probnorm(d);
  if &link=1 then c=1/(1+exp(-d));
finish cdf;

bon=x[,1];
n=x[,2];
x=x[,3:ncol(x)];

JJ=max(bon);
rows=nrow(x);
cols=ncol(x);
II=rows/JJ;

x_xi=j(rows,jj-1,0);
d_theta=beta[cols+1:cols+&cutpoint];
beta=beta[1:cols];
theta=theta+d_theta;


nn=(i(ii)@j(1,jj,1))*n;
zeta=j(rows,1,0);
r=j(rows,1,0);

eta=x*beta;
if bool_z>0 then eta=eta+z*u;

x_theta=j(rows,jj-1,0);
do i=1 to ii;
  do j=1 to jj;
    count=(i-1)*jj + j;
    if j<jj then x_theta[count,j]=1;
  end;
end;



  t_star= x_theta*theta- eta;



  do i=1 to II;
    do j=1 to JJ;
      count=(i-1)*jj + j;

      if j=jj then a=0;
      if j<JJ then call phi(a,t_star[count]);

      if j=1 then b=0;
      if j>1   then call phi(b,t_star[count-1]);

      if j=1  then d1=0;
      if j>1  then call phi(d1,t_star[count-1]);

      if j=jj then d2=0;
      if j<jj then call phi(d2,t_star[count]);

      c=d1-d2;
      if abs(c)<1e-6 then zeta[count]=.;           /*ignore offending observations*/
      else do;
        if j<jj then x_xi[count,j]=    a/c;
        if j>1  then x_xi[count,j-1]= -b/c;

        pi=0;
        if j<jj then do;
           call cdf(dum,t_star[count]);
           pi=pi+dum;
        end;
        if j=jj then pi=1;
        if j>1  then do;
          call cdf(dum, t_star[count-1]);
          pi=pi-dum;
        end;
        zeta[count]=eta[count]+(n[count]-nn[i]*pi)/c/nn[i];
        if pi<>0 then r[count]=nn[i]*c*c/pi;
      end; /*else*/
    end;
  end;

create dat_zeta var{zeta};
append;
close dat_zeta;


varnames='x1':&x_last;
create dat_x_xi from x_xi [colname=varnames];
append from x_xi;
close dat_x_xi;

create dat_r var{r};
append;
close dat_r;

create dat_thet var{theta};
append;
close dat_thet;




quit;

data temp;
merge design dat_zeta dat_x_xi dat_r;

proc mixed data=temp maxiter=&mix_iter scoring=&scoring;
class &class;
weight r;
model zeta=&fix x1-x&cutpoint/noint;
&_random_
repeated;
parms/parmsdata=cov hold=&eqcons;
make 'covparms' out=cov2;
run;

data cov2;
set cov2;
*if estimate<0.0001 then estimate=0.0001;   /*constraint removed, 21.5.2002*/


proc mixed data=temp;    /*noprofile added for 6.12*/
class &class;
weight r;
model zeta=&fix x1-x&cutpoint/solution noint;
&_random_
repeated;
make 'solutionF' out=solf;
&_make_R_;
parms/parmsdata=cov2 noiter;
make 'covparms' out=cov;
run;

data eval;
merge dat_thet old_thet;
sqr_diff=(theta-old_thet)**2;

proc univariate data=eval noprint;
var sqr_diff;
output out=eval2 sum=ssq;
run;


data eval2;
set eval2;
if ssq<&converge then dum=&iter;
    else dum=&it;
call symput('it',dum);
run;


%end;


proc printto print=print;
proc printto log=log;


proc mixed data=temp;     /*noprofile added for 6.12*/
class &class;
weight r;
model zeta=&fix x1-x&cutpoint/solution chisq noint;
&_random_
repeated;
parms/parmsdata=cov noiter;
lsmeans &lsmeans/&lsm_opt;
run;

proc print data=dat_thet;

data;
set eval2;
if ssq<&converge then put 'CONVERGENCE CRITERION MET';
else put '!!WARNING: DID NOT CONVERGE!!';
run;


%mend;