/* Categorical Linear Models                    */
   /* by Least Squares and Maximum Likelihood      */
   /*  CATLIN                                      */
   /*  Input:                                      */
   /*     n the s by p matrix of response counts   */
   /*     x the s by r design matrix               */

proc iml ;
start catlin;

   /*---find dimensions---*/
   s = nrow(n);                     /* number of populations */
   r = ncol(n);                       /* number of responses */
   q = r-1;                     /* number of function values */
   d = ncol(x);               /* number of design parameters */
   qd = q*d;                   /* total number of parameters */

   /*---get probability estimates---*/
   rown = n[,+];                               /* row totals */
   pr = n/(rown*repeat(1,1,r));     /* probability estimates */
   p = shape(pr[,1:q] ,0,1);     /* cut and shaped to vector */
   print "INITIAL PROBABILITY ESTIMATES" ,pr; 

      /*   estimate by the GSK method   */

      /* function of probabilities */
   f = log(p)-log(pr[,r])@repeat(1,q,1);

      /* inverse covariance of f */
   si = (diag(p)-p*p`)#(diag(rown)@repeat(1,q,q));
   z = x@i(q);                     /* expanded design matrix */
   h = z`*si*z;                      /* crossproducts matrix */
   g = z`*si*f;                              /* cross with f */
   beta = solve(h,g);              /* least squares solution */
   stderr = sqrt(vecdiag(inv(h)));        /* standard errors */
   run prob;
   print ,"GSK ESTIMATES" , beta stderr ,pi;

      /*   iterations for ML solution   */
   crit = 1;
   do it = 1 to 8 while(crit>.0005);/* iterate until converge*/

      /* block diagonal weighting  */
      si = (diag(pi)-pi*pi`)#(diag(rown)@repeat(1,q,q));
      g = z`*(rown@repeat(1,q,1)#(p-pi));        /* gradient */
      h = z`*si*z;                                /* hessian */
      delta = solve(h,g);            /* solve for correction */
      beta = beta+delta;             /* apply the correction */
      run prob;                    /* compute prob estimates */
      crit = max(abs(delta));       /* convergence criterion */
   end;
   stderr = sqrt(vecdiag(inv(h)));        /* standard errors */
   print , "ML Estimates", beta stderr, pi;
   print , "Iterations" it "Criterion" crit;
finish catlin;

   /*   subroutine to compute new prob estimates @ parameters  */
start prob;
   la = exp(x*shape(beta,0,q));
   pi = la/((1+la[,+] )*repeat(1,1,q));
   pi = shape(pi,0,1);
finish prob;

   /*---prepare frequency data and design matrix---*/
n= { 58 11 05,
     75 19 07,
     49 14 10,
     58 17 08,
     33 18 15,
     45 22 10,
     15 13 15,
     39 22 18,
     04 12 17,
     05 15 08};     /* frequency counts*/

x= { 1 1 1 0 0 0,
     1 -1 1 0 0 0,
     1 1 0 1 0 0,
     1 -1 0 1 0 0,
     1 1 0 0 1 0,
     1 -1 0 0 1 0,
     1 1 0 0 0 1,
     1 -1 0 0 0 1,
     1 1 -1 -1 -1 -1,
     1 -1 -1 -1 -1 -1};    /* design matrix*/

run catlin;