/**************************************************
 * Program for two-sample-LR-test of equality of  *
 * means with monotone missing data.              *
 **************************************************/

proc iml;

/*data from example 16.6.1*/
X={137 245,
   132 260,
   141 276,
   142 299,
   128 239,
   147 262,
   136 278,
   128 255,
   128 244,
   146 276,
   128 242,
   147 263,
   121   .,
   138   .,
   138   .,
   150   .};

Y={184 305,
   133 237,
   166 300,
   162 273,
   163 297,
   160 308,
   166 301,
   141 208,
   146   .,
   171   .};

/*Sample size and number of variables:*/
nx=nrow(X);
ny=nrow(Y);
p=ncol(X);
print "no. of obs. and variables:" nx ny p;

 xindex=((X^=.)#(1:p))[,<>];
 yindex=((Y^=.)#(1:p))[,<>];
do i=1 to p;
 nxi=(xindex>=i)[+]; 
/*number of sample units that have an obs. on the 
  ith variable (and, because of monotony, on every 
   lower variable)*/
 nyi=(yindex>=i)[+];
 ai=(nxi#nyi)/(nxi+nyi);
 ti=nxi+nyi;
 Xi=X[1:nxi,1:i];
 Yi=Y[1:nyi,1:i];
    muvecx=(Xi[:,])`;
    muvecy=(Yi[:,])`;
 Smat=Xi`*(I(nxi)-j(nxi,nxi,1/nxi))*Xi+
      Yi`*(I(nyi)-j(nyi,nyi,1/nyi))*Yi;
 Numat=Smat+ai#(muvecx-muvecy)*(muvecx-muvecy)`;
 if i>1 then do; Smatd=Smat[1:(i-1),1:(i-1)];
                 Numatd=Numat[1:(i-1),1:(i-1)]; end;
 if i=1 then do; lambda=(det(Smat)/det(Numat))##(ti/2);
                 helprho=(2#i+3)/ti; end;
 if i>1 then do; lambda=lambda#
(det(Smat)#det(numatd)/(det(Numat)#det(Smatd)))##(ti/2);
                 helprho=helprho+(2#i+3)/ti; end;
end;
rho=1-helprho/(2#p);
ll=-2#log(lambda);
lladj=rho#ll;
pvalue=1-probchi(lladj,p);
print lambda ll;
print "Test of H0:mu=0:" lladj p pvalue;
quit;
run;