/*----------------------------------------------------------------*/
/*                                                                */
/* This program is based on the article "Response Transformations */
/* in Repeated Measures and Growth Curve Models," by R. J. Boik,  */
/* C. Todd, and S. K. Hyde (Communications in Statistics---Theory */
/* and Methods, 2000, 29, 669-773).                               */
/*                                                                */
/*----------------------------------------------------------------*/

DM 'LOG;CLEAR;OUT;CLEAR;';
options ls=79 pagesize=60 nodate;

data raw;
 infile 'wei.dat';
 input grp y1-y5;

/*------------------------------------------------------------*/
/*--------- PROC GLM: Analysis of Repeated Measures  ---------*/
/*------------------------------------------------------------*/

/* Executation of proc glm is optional. */

/*
proc glm data = raw;
class grp;
model y1-y5 =grp /nouni;
repeated time 5;
*/

/*------------------------------------------------------------*/
/*----------- PROC GLMMOD: Create Design Matrix --------------*/
/*------------------------------------------------------------*/

/* Executation of proc glmmod is required. */

proc glmmod data = raw outparm=outp outdesign=outX noprint;
class grp;
model y1-y5=grp;
run;


/*-----------------------------------------------------------------*/
/*-- PROC IML: Compute MLE of Lambda & Perform Lack of Fit Tests --*/
/*-----------------------------------------------------------------*/

proc iml;

/*------------------------------------------------------------*/
/*------------------ User Defined Variables ------------------*/
/*------------------------------------------------------------*/

/* This iml program computes the MLE of the Box-Cox or folded-power   */
/* transformation parameter, lambda, in a repeated measures or in a   */
/* growth curve model.  The user must customize the program for her/  */
/* his particular application.  A maximum of three within-subject     */
/* factors are allowed.  The following variables can be assigned      */
/* customized values:                                                 */
/*                                                                    */
/* levels: This is a row vector that specifies the number of levels   */
/*    of the within subject factors.  The levels should be specified  */
/*    in the same order as on the REPEATED statement in PROC GLM. For */
/*    example, if there are two within-subject factors with 5 and 4   */
/*    levels, then levels would be specified as "levels = {4 5};" or  */
/*    as "levels = {5 4};" depending on which index varies fastest.   */
/*    If the linear model has no repeated measures, then levels       */
/*    should be specified as "levels = {1};."                         */
/*                                                                    */
/* Polydeg: This is a matrix that specifies the polynomial degree of  */
/*    the growth curve that is desired for each within-subject factor.*/
/*    To obtain a repeated measures model, a saturated polynomial     */
/*    curve should be specified.  That is, the polynomial degree for  */
/*    each factor should be one fewer that the number of levels of    */
/*    factor.  For example, if "levels = {5 4};" then a repeated      */
/*    measures model is obtained by setting "Polydeg = {4 3};."  If   */
/*    linear model has no repeated measures, then Polydeg should be   */
/*    specified as "Polydeg = {0};."                                  */
/*                                                                    */
/*    To obtain a growth curve model, a lower order polynomial curve  */
/*    should be specified for one or more within-subject factor.      */
/*    Each row of Polydeg specifies the set of polynomial degrees for */
/*    a single growth curve model.  For example, if "levels = {4 5};,"*/
/*    then "Polydeg = {3 4, 3 3, 2 3};" specifies the polynomial      */
/*    degrees for three growth curve models.  In the first model,     */
/*    Factor 1 is modeled as a third degree polynomial and Factor 2   */
/*    is modeled as a fourth degree polynomial.  Note that this is a  */
/*    saturated model and is equivalent to a repeated measures model. */
/*    In the second model, both Factors are modeled as third degree   */
/*    polynomials.  A polynomial curve with degree zero contains an   */
/*    intercept term only.                                            */
/*                                                                    */
/* time1: This is a row vector that specifies the time points at      */
/*    which observations were made on Factor 1.  The variables        */
/*    "time2" and "time3" are defined analogously.  If the levels of  */
/*    Factor 1 represent equal time intervals, then time1 can be      */
/*    specified as "time1 = 1:levels[1];."  If the levels of Factor 1 */
/*    represent time points 1, 3, 4, 6, then time1 can be specified   */
/*    as "time1 = {1 3 4 6};." If there is only one within-subject    */
/*    factor, then specify "time2 = {0};" and "time3 = {0};."         */
/*                                                                    */
/* confid: This is a scalar that represents the desired confidence    */
/*    coefficient.  Typically, confid is set to 0.95 or 0.99.  the    */
/*    value of confid must be in the interval (0,1).                  */
/*                                                                    */
/* boxcox: This is a character variable.  Specify "boxcox = 'yes';"   */
/*    if the Box-Cox family of transformations is to be employed.     */
/*    Otherwise, specify "boxcox = 'no';." Either boxcox or foldpow,  */
/*    but not both, must be specified as 'yes.'                       */
/*                                                                    */
/* foldpow: This is a character variable.  Specify "foldpow = 'yes';" */
/*    if the folded-power family of transformations is to be employed.*/
/*    Otherwise, specify "foldpow = 'no';."                           */
/*                                                                    */
/* guesslam: This is a scalar that represents an initial guess for    */
/*    lambda.  The value 1 appears to work well, but for some data    */
/*    sets, the initial guess might need to be modified.              */
/*                                                                    */
/* converge: This is a scalar convergence criterion used when         */
/*    computing the MLE of the covariance matrix under local          */
/*    sphericity.  Iteration stops when improvement in the evaluated  */
/*    log likelihood function is smaller than "converge."  The value  */
/*    1e-6 appears to work well and might not need to be modified.    */
/*                                                                    */
/* small:  This is a scalar that is used to determine the rank of a   */
/*    matrix.  The rank of a matrix is computed as the number of      */
/*    singular values that are larger than the maximum singular value */
/*    multiplied by small.  The value 1.e-10 appears to work well and */
/*    might not need to be modified.                                  */

levels = {5};
Polydeg={4,3,2,1,0};
time1 = {0 6 12 20 24}; time2 = {0}; time3 = {0};
confid = .95;
boxcox='yes';
foldpow='no';
guesslam = 1;
converge = 1e-6;
small = 1e-10;

/*------------------------------------------------------------*/
/*------------End of User Defined Variables ------------------*/
/*------------------------------------------------------------*/

/*-------------------------------------------------*/
/*-------------- Define Sub-Modules ---------------*/
/*-------------------------------------------------*/

/*-------------------------------------------------*/
/* Module to stack columns of matrix into a vector */
/*-------------------------------------------------*/    
start vec(mat);
  n = nrow(mat);
  p = ncol(mat);
  tot = n*p;
  Y = shape(mat,tot,1);
  return(Y);
finish vec;

/* ----------------------------------------------------------*/
/* Module to unstack vector into a matrix with the specified */
/* number of rows.  No checking is made if elements are lost */
/* or duplicated at the end                                  */
/* ----------------------------------------------------------*/    
start dvec(mat,nr);
  tot = nrow(mat);
  nc = tot/nr;
  Y = shape(mat,nr,nc);
  return(Y);
finish dvec;

/*-----------------------------------------*/
/* Module to compute the rank of a matrix  */
/*-----------------------------------------*/
start rankM(A) global(small);
  m = nrow(A); n = ncol(A);
  call svd(U,S,V,A);
  if m = 1
    then S = S[1];
    else if m = 0 then S = {0}; 
  tol = max(S) * small;
  r = max(loc(S>tol)); /*rank of A = number of nonzero singular values */
  return(r);
finish rankM;

/*-------------------------------------------*/ 
/* Module to compute the commutation matrix. */
/*-------------------------------------------*/
start commute(p,q);
  Ip = I(p);
  Iq = I(q);
  k = shape({0},p*q,p*q);
  if p <= q
    then
      do i = 1 to p;
        k = k + (Ip[,i]`@Iq)@Ip[,i];
      end;
    else
      do i = 1 to q;
        k = k + (Iq[,i]@Ip)@Iq[,i]`;
      end;
  return(k);
finish commute;

/*-------------------------------------------*/
/* Module to compute elementary matrices     */
/*-------------------------------------------*/
start makeE(m,i);
  nr = max(ncol(m),nrow(m));
  if nr = 1 then
    do;
      E = I(m[1]);
      return(E);
    end;
  if i = 1 then do; 
    s2 = m[2:nr][+,+];
    E = I(m[1]) // J(s2,m[1],0);
    return(E);
  end;
  if i = nr then do;
    s1 = m[1:nr-1][+,+];
    E = J(s1,m[nr],0) // I(m[nr]);
    return(E);
  end;
  s1 = m[1:i-1][+,+];
  s2 = m[i+1:nr][+,+];
  E = J(s1,m[i],0) // I(m[i]) // J(s2,m[i],0);
  return(E);
finish makeE;

/*-------------------------------------------*/
/* Module to compute the C and G matrices    */
/*-------------------------------------------*/
start makeCG(C,G,q,d,levels,w) global(time1,time2,time3); 
  nfact = max(nrow(levels),ncol(levels));
  
  if nfact = 1 then do;
    t1 = J(levels[1],1,1)/sqrt(levels[1]);
    if w[1] = 0 then do;
       G=t1;
       d=0;
       q=0;
       return;
    end;
    G = orpol(time1,w[1]);
    C = G[,2:w[1]+1];
    q = w[1];
    d = 1;
    return;
  end;
  else if nfact = 2 then do;
    t1 = J(levels[1],1,1)/sqrt(levels[1]);
    t2 = J(levels[2],1,1)/sqrt(levels[2]);
    if w[+] = 0 then goto AB00;
    if w[1]*w[2] > 0 then goto AB11;
    if w[1] > 0 then goto AB10;
    if w[2] > 0 then goto AB01;
    AB00:
       G=t1@t2;
       d=0;
       q=0;
       return;
    AB11:
       G1 = orpol(time1,w[1]);
       G2 = orpol(time2,w[2]);
       C1 = G1[,2:w[1]+1];
       C2 = G2[,2:w[2]+1];
       G = G1@G2;
       C=  (C1@t2) || (t1@C2) || (C1@C2);
       q = w[1] || w[2] || w[1]*w[2];
       d = 3;
       return;
    AB10:
       G1 = orpol(time1,w[1]);
       C1 = G1[,2:w[1]+1];
       G=G1@t2;
       C=C1@t2;
       q = w[1];
       d = 1;
       return;
    AB01:
       G2 = orpol(time2,w[2]);
       C2 = G2[,2:w[2]+1];
       G=t1@G2;
       C=t1@C2;
       q = w[2];
       d = 1;
       return;
  end;
  else if nfact = 3 then do;
    t1 = J(levels[1],1,1)/sqrt(levels[1]);
    t2 = J(levels[2],1,1)/sqrt(levels[2]);
    t3 = J(levels[3],1,1)/sqrt(levels[3]);
    if w[+] = 0 then goto ABC000;
    if w[1]*w[2]*w[3] > 0 then goto ABC111;
    if w[1]*w[2] > 0 then goto ABC110;
    if w[1]*w[3] > 0 then goto ABC101;
    if w[2]*w[3] > 0 then goto ABC011;
    if w[1] > 0 then goto ABC100;
    if w[2] > 0 then goto ABC010;
    if w[3] > 0 then goto ABC001;
    ABC000:
       G=(t1@t2)@t3;
       q=0;
       d=0;
       return;
    ABC111:
       G1 = orpol(time1,w[1]);
       G2 = orpol(time2,w[2]);
       G3 = orpol(time3,w[3]);
       C1 = G1[,2:w[1]+1];
       C2 = G2[,2:w[2]+1];
       C3 = G3[,2:w[3]+1];
       G = (G1@G2)@G3;
       C = (C1@t2)@t3 || (t1@C2)@t3 || (t1@t2)@C3 || (C1@C2)@t3 ||
           (C1@t2)@C3 || (t1@C2)@C3 || (C1@C2)@C3;
       q = w[1] || w[2] || w[3] || w[1]*w[2] || w[1]*w[3] || w[2]*w[3]
            || w[1]*w[2]*w[3];
       d = 7;
       return;
    ABC110:
       G1 = orpol(time1,w[1]);
       G2 = orpol(time2,w[2]);
       C1 = G1[,2:w[1]+1];
       C2 = G2[,2:w[2]+1];
       G = (G1@G2)@t3;
       C = (C1@t2)@t3 || (t1@C2)@t3 || (C1@C2)@t3;
       q = w[1] || w[2] || w[1]*w[2];
       d = 3;
       return;
    ABC101:
       G1 = orpol(time1,w[1]);
       G3 = orpol(time3,w[3]);
       C1 = G1[,2:w[1]+1];
       C3 = G3[,2:w[3]+1];
       G = (G1@t2)@G3;
       C = (C1@t2)@t3 || (t1@t2)@C3 || (C1@t2)@C3;
       q = w[1] || w[3] || w[1]*w[3];
       d = 3;
       return;
    ABC011:
       G2 = orpol(time2,w[2]);
       G3 = orpol(time3,w[3]);
       C2 = G2[,2:w[2]+1];
       C3 = G3[,2:w[3]+1];
       G = (t1@G2)@G3;
       C = (t1@C2)@t3 || (t1@t2)@C3 || (t1@C2)@C3;
       q = w[2] || w[3] || w[2]*w[3];
       d = 3;
       return;
    ABC100:
       G1 = orpol(time1,w[1]);
       C1 = G1[,2:w[1]+1];
       G = (G1@t2)@t3;
       C = (C1@t2)@t3;
       q = w[1];
       d = 1;
       return;
    ABC010:
       G2 = orpol(time2,w[2]);
       C2 = G2[,2:w[2]+1];
       G = (t1@G2)@t3;
       C = (t1@C2)@t3;
       q = w[2];
       d = 1;
       return;
    ABC001:
       G3 = orpol(time3,w[3]);
       C3 = G3[,2:w[3]+1];
       G = (t1@t2)@G3;
       C = (t1@t2)@C3;
       q = w[3];
       d = 1;
       return;
  end;
finish makeCG;

/*-----------------------------------------------------------------*/
/* Module to compute the mle of Omega under local sphericity       */
/* The algorithm is based on Fisher scoring                        */
/*-----------------------------------------------------------------*/

start mlesigma(A,n,m,converge) global(sig_iter,GG_save, G2_save);
   mt = m[+,+];
   Vhat = A/n;
   d = max(ncol(m),nrow(m));
   if d = 1 then do;
     v = trace(Vhat)/mt;
     Sigmahat = v*I(mt);
     return(Sigmahat);
   end;
   if m[##] = d then  return(Vhat);
   if sig_iter = 0 then do;
      GG = J(1,mt*mt,0); /* Initialize with a row of zeros */
      do i =1 to d;
         Ei = makeE(m,i);
         do j = i+1 to d;
            Ej = makeE(m,j);
            GG = GG // (Ej` @ Ei`);
         end;
      end;
      nr = nrow(GG);
      GG = GG[2:nr,];     /* Remove row of zeros */
      mmt=mt;
      G2 = (I(mt*mt)+commute(mt,mmt))*GG`;
      GG_save=GG;
      G2_save=G2;
   end;     
   old = {0};
   do i = 1 to d;
      if i = 1 then s = 0;
      else s = m[1:i-1][+,+];
      mi = m[i];
      Temp=A[s+1:s+mi,s+1:s+mi];
      Temp=(Temp+Temp`)/2;
      call eigen(eigs,dumb,Temp);
      vv  = eigs[<>]/n;
      old = old // vv;
      Vhat[s+1:s+mi,s+1:s+mi] = I(mi)*vv;
   end;
   old = old[2:d+1];
   if sig_iter > 0 then do;
      GG=GG_save;
      G2=G2_save;
   end;
   sig_iter=sig_iter+1;
   old2 = GG*vec(Vhat);
   old = old // old2;
   del = 1;
   guess = Vhat;
   Vi = inv(Vhat);
   logl = -trace(A*Vi)/2 - n*log(det(Vhat))/2;
   d1 = (mt*mt-m[##])/2;

   iter = 0;

   do while(del > converge);
     H11 = J(d,d,0);
     TT = Vi*A*Vi/n;
     mi1 = m[1];
     Vi1 = Vi[,1:mi1];
     Vii = Vi[1:mi1,1:mi1];
     Tii = TT[1:mi1,1:mi1];
     h1 = trace(Tii-Vii)/2;
     U = Vi1*Vi1`;
     H21 = GG*vec(U);
     do j1 = 1 to d;
       if j1 = 1 then sj1 = 0;
       else sj1 = m[1:j1-1][+,+];
       mj1 = m[j1];
       Vi1j1 = Vi[1:mi1,sj1+1:sj1+mj1];
       v = vec(Vi1j1);
       H11[1,j1] = v`*v/2;
       H11[j1,1] = H11[1,j1];
     end;

     do i1 = 2 to d;
       si1 = m[1:i1-1][+,+];
       mi1 = m[i1];
       Vi1 = Vi[,si1+1:si1+mi1];
       Vii = Vi[si1+1:si1+mi1,si1+1:si1+mi1];
       Tii = TT[si1+1:si1+mi1,si1+1:si1+mi1];
       h1 =  h1 // trace(Tii-Vii)/2;
       U = Vi1*Vi1`;
       H21 = H21 || GG*vec(U);
       do j1 = i1 to d;
         sj1 = m[1:j1-1][+,+];
         mj1 = m[j1];
         Vi1j1 = Vi[si1+1:si1+mi1,sj1+1:sj1+mj1];
         v = vec(Vi1j1);
         H11[i1,j1] = v`*v/2;
         H11[j1,i1] = H11[i1,j1];
       end;
     end;

     H22 = GG*(Vi@Vi)*G2;
     H = - ((H11 || H21`) // (H21 || H22));
     h2 = GG*vec(TT-Vi);
     hh = h1 // h2;
     change = -inv(H)*hh;
     c = 2;
     dd = -1;
     do while(dd < -1.e-10);
       c = c/2;
       new = old+c*change;
       new1 = new[1:d];
       new2 = new[d+1:d+d1];
       Vhat = dvec(GG`*new2,mt);
       Vhat = Vhat + Vhat`;
       do i = 1 to d;
         if i = 1 then s = 0;
         else s = m[1:i-1][+,+];
         mi = m[i];
         Vhat[s+1:s+mi,s+1:s+mi] = I(mi)*new1[i];
       end;
       call eigen(eigs,dumb,Vhat);
       dd = eigs[><];
       if dd > 0 then do;
          dd1=dd;
          vi=inv(Vhat);
          newl = -trace(A*Vi)/2-log(det(Vhat))*n/2;
          dd = newl - logl;
       end;
     end;
     iter = iter + 1;
     old = new;
     del = abs(logl-newl);
     logl = newl;
   end;
   return(Vhat);
finish mlesigma;

/*---------------------------------------------------------------*/
/* Module for computing the log likelihood for the following     */
/* four covariance structures.  In the growth curve model, these */
/* are conditional covariance structures.                        */
/*    1.  Unconstrained.                                         */
/*    2.  Local Sphericity                                       */
/*    3.  Local Sphericity with Independence.                    */
/*    4.  Omnibus Sphericity.                                    */
/*---------------------------------------------------------------*/
start lnlike(lam) global(Y,X,XPXiX,C,G,logg,n,q,t,wd,ci,
                  boxcox,foldpow,converge,del_max,maxlnlik,lnlik);
  if boxcox = 'yes'  then do;
     lnJ=(lam-1)*logg;
     if abs(lam) > .0001 then Z = (Y##lam-1)/lam;
     else Z = log(Y);
  end;
  if foldpow = 'yes' then do;
     lnJ=log(Y##(lam-1) + (J(n,t,1)-Y)##(lam-1));
     lnJ=lnJ[+,+];
     if abs(lam) > .0001 then Z = (Y##lam -(J(n,t,1)-Y)##lam)/lam;
     else Z= log(Y/(J(n,t,1)-Y));
  end;
  B = XPXiX*Z;
  S = (Z`*Z - Z`*X*B)/n;
  grow = 'no';
  if ncol(G) < t then grow = 'yes';
  SigGQ = S;
  if grow = 'yes' then do; 
     Sinv = inv(S);
     temp1 = G*ginv(G`*Sinv*G);
     Bg = B*Sinv*temp1;
     SigGQ = S + (B - Bg*G`)`*X`*X*(B - Bg*G`)/n;
  end;
  if q[+] = 0 then do;
      lnlik = n*t/2 + n/2*log(det(SigGQ)) - lnJ;
           /*Confidence interval calculation*/
       if ci=1 then lnlik = 10000*(lnlik - maxlnlik)**2;
       else lnlik=lnlik-del_max;
  return(lnlik);
  end;
  /* Compute Omegahat and Sigmahat for various covariance structures */
  if wd = 1 then Sigmag = SigGQ; /* Unconstrained */
  else do;
    Sig2=S;
    if grow='yes' then Sig2 = temp1*G`;
    Omegat = C`*Sig2*C;
    Omegati = inv(Omegat);
    /* Local Sphericity*/
    if wd = 2 then Omegah = mlesigma(n*Omegat,n,q,converge);
    /* Local with independence*/
    if wd = 3 then do;
      mt = q[+,+];
      Omegah = J(mt,mt,0);
      s2 = 0;
      do i = 1 to max(nrow(q),ncol(q));
        s1 = s2 + 1; 
        s2 = s2 + q[i];
        thetai = trace(Omegat[s1:s2,s1:s2])/q[i];
        Omegah[s1:s2,s1:s2] = thetai*I(q[i]);
      end;
    end;
    /* Omnibus */
    if wd = 4 then do;
      sizeO = ncol(C);
      Omegah = trace(Omegat)*I(sizeO)/sizeO;
    end;
    Sigmag = SigGQ + Sig2*C*Omegati*(Omegah-Omegat)*Omegati*C`*Sig2;
  end;

  lnlik = n*t/2 + n/2*log(det(Sigmag)) - lnJ;
  
  /*Confidence interval calculation*/
  if ci=1 then lnlik = 10000*(lnlik - maxlnlik)**2;
  else lnlik=lnlik-del_max;
  return(lnlik);
finish lnlike;

/*-------------------------------------------------*/
/*------------------ Main Module ----------------- */
/*-------------------------------------------------*/
do;
use outx;
read all var _num_ into X;
t = levels[#,#];          /*Product of levels of within-subject factors*/
n = nrow(X);
p2 = ncol(X);
Y = X[,1:t];
X = X[,t+1:p2];
r = rankM(X);
maxY=1;
error='no';
if boxcox = 'yes' then do;
   maxY=max(Y);   /* Normalization of Data for numerical stability*/
   Y = Y/maxY;    /* for Box-Cox family only                      */
   if foldpow = 'yes' then do;
      print 'Input error: Both transformation families have been selected.';
      print 'Select either Box-Cox or folded-power,';
      error = 'yes';
   end;
end;
minY=min(Y);
if minY <= 0 then do;
  print 'Data Error: One or more data values are non-positive';
  error='yes';
end;
if foldpow='yes' then do;
  maxYY = max(Y);
  if maxYY > 1 then do;
      print 'Data Error: One or more data values are greater than 1';
      print 'For the folded power family, all data values must be  ';
      print 'in (0, 1).';
      error = 'yes';
  end;
end;
checkdeg=J(nrow(Polydeg),1,1)*levels-Polydeg;
if min(checkdeg) < 1 then do;
   print 'Input Error: Polynomial degree cannot exceed';
   print 'the number of levels minus one.';
   print levels;
   print Polydeg;
   error = 'yes';
end;
if error = 'yes' then goto exit;
logg = sum(log(Y));
XPXiX = ginv(X`*X)*X`;

/*----------------------------------------------------------------*/    
/* Main module to compute mle of lambda (transformation parameter)*/
/* and a confidence interval for lambda under various covariance  */
/* structures.  In the growth curve model, these are conditional  */
/* covariance structures.                                         */
/*    wd = 1 -->  Unconstrained                                   */
/*    wd = 2 -->  Local Sphericity                                */
/*    wd = 3 -->  Local Sphericity with Independence              */
/*    wd = 4 -->  Omnibus Sphericity                              */
/*----------------------------------------------------------------*/

sizew = nrow(Polydeg);
crit = 2*gaminv(confid,.5); /*percentile of Chi-square 1 distribution*/
critn = sqrt(crit);         /*percentile of standard normal dist*/
do wj = 0 to sizew;

  if wj = 0 then do;
    /*Find the mle for unconstrained full model */
    wd = 1;                   /*Unconstrained*/
    w = levels - J(1,ncol(levels),1);
    run makeCG(C,G,q,d,levels,w);
    f = ncol(G);
    guess = guesslam;         /*Initial guess for lambda   */
    ci = 0;                   /* Flag for finding the mle  */
    maxlnlik = 0;
    sig_iter=0;
    del_max=0;
    call nlpnrr(rc,lambda,"lnlike",guess);
    lam_full=lambda;
    guess=lambda;
    bdim = r*t + t*(t+1)/2 + 1;
    blnlik = -lnlik -n*t*log(maxY);
    lnlik_m=-lnlik;
  end;
  else do;
    lofv = { 0 0 0 0 0 };
    estm = { 0 0 0 0 0 };
    w = Polydeg[wj,];
    run makeCG(C,G,q,d,levels,w);
    f = ncol(G);
    nwd=4;
    if max(nrow(levels),ncol(levels)) = 1 then nwd=2;
    if levels[+] = 1 then nwd=1;
    do wdtemp = 1 to nwd;
      wd=wdtemp;
      if nwd=2 & wdtemp=2 then wd=4;
      /* Find the mle estimate of the transformation parameter */
      ci = 0;                   /* Flag for finding the mle */
      maxlnlik = 0;
      sig_iter=0;
      full='no';
      if wd=1 & f=t then do;
         lambda=lam_full;
         logl=lnlik_m;
         full = 'yes';
      end;
      if full = 'no' then do;
         del_max=-lnlik_m;
         call nlpnrr(rc,lambda,"lnlike",guess);
         guess=lambda;
         logl = -lnlik-del_max;
      end;
      /* Estimate the standard error of lambdahat */
      sig_iter=0;
      del_max=0;
      call nlpfdd(maxlik,gradient,se,"lnlike",lambda);
      se = 1/sqrt(se);
      
      /* Find lower confidence interval limit for lambda */
      ci = 1;                   /* Flag for finding an endpoint of the ci */
      maxlnlik = maxlik + crit/2;
      guess1 = lambda - critn*se;
      sig_iter=0;
      low=guess1;
      tc=J(13,1,0); tc[1]=200; tc[2]=500; tc[3]=1.e-6; tc[4]=1.e-8;
      tc[6]=1.e-5; tc[7]=2.2204e-16; tc[10]=1;
      call nlpnrr(rc,low,"lnlike",guess1) tc=tc;
      /* Find upper confidence interval limit for lambda */
      guess2 = lambda + critn*se;
      sig_iter=0;
      high=guess2;
      call nlpnrr(rc,high,"lnlike",guess2) tc=tc;
      /* Calculate Lack of Fit Statistic and corresponding p-value*/
      /* plus the dimension of the model and AIC   */

      /* Calculate dimensions*/
      qq = q[+,+];
      if wd = 1 then dim = r*f + t*(t+1)/2 + 1;
      if wd = 2 then dim = r*f + (t*(t+1) - qq - ssq(q))/2 + d + 1;
      if wd = 3 then dim = r*f + (t-qq)*(t+qq+1)/2 + d + 1;
      if wd = 4 then dim = r*f + (t-qq)*(t+qq+1)/2 + 1;
      logl=logl-n*t*log(maxY);
      lof = -2*(logl - blnlik);
      df = bdim - dim;
      if lof > 0
        then pv = 1 - probchi(lof,df);
        else do;
          df=.;
          pv = .;
          lof = .;
        end;
      aic = logl - dim;
      lofv = lofv // ( lambda || logl || lof || df || pv );
      estm = estm // ( lambda || low || high || dim || aic );

    end;

    lofv = lofv[2:nwd+1,];
    estm = estm[2:nwd+1,];
    rows1 = {'Unconstrained','Local Sphericity','Local w Indep','Omnibus'};
    rows2 = {'Unconstrained','Omnibus Spher.'};
    rows=rows1;
    if nwd=2 then rows=rows2;
    if w = levels - J(1,ncol(levels),1) then do;
         if levels[+] > 1 then print 'Repeated Measures Model';
         else print 'Linear Model with No Repeated Measures';
         mattrib lofv
            rowname = rows
            colname = {'Estimate','LogLik','LoF','df','p-value'}
            label = 'Cov. Structure'
            format = F9.4;
    end;
    else do;
       print 'Growth Curve Model';
       mattrib lofv
         rowname = rows
         colname = {'Estimate','LogLik','LoF','df','p-value'}
         label = 'Cond Cov Struct'
         format = F9.4;
       mattrib w
          label = ''
          colname = ''
          rowname={'Polynomial Degree ='}
          format=F2.0;
       print w;
    end;

    print lofv;
    clevel=confid*100;
    dum2 = char(clevel,3);
    lowlim = concat(dum2,'% Low');
    highlim = concat(dum2,'% High');
    cols = {'Estimate'} // lowlim // highlim // {'Dimen', 'AIC'};
    mattrib estm
        rowname = rows
        colname = cols
        label = ''
        format = F9.4;
    print estm;
    if wj < sizew then do;
    print /;
  end;
end;

end;
/*------------------------------------------------------------*/

exit:
if error = 'yes' then abort;
end;
quit;