Calculation of level probabilities


Abstract

The following SAS/IML-program calculates level probabilities, which are required in manyfold applications. The calculation of level probabilities is based on a recursive relationship of orthant normal probabilities with a tridiagonal correlation matrix (see Robertson et al., 1988, p. 77). A SAS/IML module for calculating the orthant probabilities up to dimension 11 can be found here. A required SAS/IML module for calculating the maximum likelihood estimates under total order restriction can be found here. In the following the recursive algorithm decribed by Bretz and Seidel (1999) is presented. The use of the algorithm is demonstrated by an data set example for calculating the p-value of the likelihood ratio test under total order restriction. This example is easily extended to the general unbalanced designs. The invokation of the modules below by
 
 
 
n={6 6 6 6};
xbar={-76.2 -73.5 -73.4 -74.4};
s=60.078;

p_value=LRT(xbar,s,n);
print p_value [format=6.4];

 /* sample size vector */
/* arithmetic means    */
/* mean square error  */

yields the following output:
 
 
 P_VALUE
       0.4634

 

SAS/IML-Source

start plkwunba(gewicht,p_lkw);
  k=ncol(gewicht);
  p_1_k=j(k-1,k-1,1);
  p_1_k[2,]=j(1,k-1,0.5);
  p_lkw=j(k,1,0);

  do t=3 to k;
    do i=2 to t;
      mengen=0;
      anz=erzall(t,i, mengen);
      nrowmen=nrow(mengen);

      if t<k then index=j(nrowmen,t,1)||j(nrowmen,k-t,0);
      else index=j(nrowmen,k,1);

      if i>1 then do;
        do s=1 to nrowmen;
          gew=j(1,i,0); top=0;
          prod=1;
          do u=1 to i;
            gew[u]=sum(gewicht[top+1:top+mengen[s,u]]);
            if mengen[s,u]=2 then prod=prod*0.5;
            else if mengen[s,u]>2 then prod=prod*p_1_k[mengen[s,u],top+1];
            top=top+mengen[s,u];
          end;
          if i=2 then z=0.5;
          else z=orthant(i,gew);
          if s=1 then y=z;
          else y=y//z;

          if t<k then p_1_k[t,1]=p_1_k[t,1]-z*prod;
          else p_lkw[i]=p_lkw[i]+z*prod;
        end;
      end;

      gewi=gewicht[,1:k];
      do j=2 to k-t+1;
        index=j(nrowmen,1,0)||index[,1:k-1];
        gewi=gewi[,2:ncol(gewi)];
        mengen1=mengen||index;
        if i>1 then do;
          do s=1 to nrowmen;
            gew=j(1,i,0);
            top=0;
            prod=1;
            do u=1 to i;
              gew[u]=sum(gewi[top+1:top+mengen[s,u]]);
              if mengen[s,u]=2 then prod=prod*0.5;
              else if mengen[s,u]>2 then
                prod=prod*p_1_k[mengen[s,u],top+j-1+1];
              top=top+mengen[s,u];
            end;
            if i=2 then z=0.5;
            else z=orthant(i,gew);
            if s=1 then y=z;
            else y=y//z;
            p_1_k[t,j]=p_1_k[t,j]-z*prod;
          end;
        end;
      end;
    end;
  end;
finish;


start erzall(n,k,mengen);
  feld1=j(1,k,0);
  first=0;
  mtc1=0;
  t=0;
  h=0;
  do while(mtc1<2);
    i=nexcom2(n, k, feld1, mtc1,t,h);
    if all(feld1) then
      if first=0 then do;
        mengen=feld1;
        first=1;
      end;
    else mengen=mengen//feld1;
  end;
  return(nrow(mengen));
finish;


start nexcom2(n_ges, c, feld, mtc, t, h);
  if (mtc=0) then do;
    feld[1]=n_ges;
    if(c>1)then do i=2 to c;
      feld[i]=0;
    end;
    t=n_ges;
    h=0;
  end;
  else goto marke2;

  marke1: if (feld[c]=n_ges) then mtc=2;
    else mtc=1;
    return(1);

  marke2: if (t>1)then h=0;
    h=h+1;
    t=feld[h];
    feld[h]=0;
    feld[1]=t-1;
    feld[h+1]=feld[h+1]+1;
    goto marke1;
finish;



start LRT(x,s,n);
  k=ncol(n);
  nu=sum(n)-k;
  m=MLE(x,n);
  xquer=sum(x#n/n[+]);

  x1=sum(n#(m-xquer)##2);
  x2=sum(n#(m-x)##2);
  t=x1/(x1+x2+nu*s);
  run plkwunba(n,p_lkw);

  p=0;
  do i=2 to k;
    p=p+p_lkw[i]*(1-probf(t*(sum(n)-i)/((i-1)*(1-t)),i-1,sum(n)-i));
  end;
  return(p);
finish;