/*----------------------------------------------------------------*
 | FOURFOLD.SAS   IML modules for four-fold                       |
 |                display of 2x2xK tables                         |
 | Usage:                                                         |
 |   proc iml;                                                    |
 |   %include fourfold;                                           |
 |   run fourfold( dim, table, vnames, lnames );                  |
 | where:                                                         |
 |   dim      table dimensions: {2 2 k}                           |
 |   table    two- or three-way contingency table,                |
 |            of size (2k)x2                                      |
 |   vnames   variable names, 1x3 character matrix.               |
 |            vnames[,1]=column variable                          |
 |            vnames[,2]=row variable                             |
 |            vnames[,3]=panel variable                           |
 |   lnames   category names, 3 x k character matrix.             |
 |            vnames[1,]=col categories,                          |
 |            lnames[2,]=row categories,                          |
 |            lnames[3,]=panel categories.                        |
 |global input variables:                                         |
 |   std      ='MARG' standardizes each 2x2 table to equal        |
 |            margins, keeping the odds ratio fixed.              |
 |            ='MAX' standardizes each table to a maximum         |
 |            cell frequency of 100.                              |
 |            ='MAXALL' standardizes all tables to                |
 |            max(f[i,j,k])=100.                                  |
 |   down     number of panels down each page                     |
 |   across   number of panels across each page                   |
 |   sangle   angle for side labels (0|90)                        |
 *----------------------------------------------------------------*
 *  Author:  Michael Friendly         <Friendly@VM1.YorkU.CA>     *
 * Created:  9 May 1991 19:12:12      Copyright (c) 1992          *
 * Revised:  12 Nov 1993 14:59:39                                 *
 * Version:  1.5                                                  *
 *----------------------------------------------------------------*/
goptions gunit=pct;
start fourfold(dim, table, vnames, lnames)
  global (std, down, across, name, sangle );
  if type(std  )  ^='C' then std='MARG';
  if type(down)   ^='N' then down=2;
  if type(across) ^='N' then across=1;
  if type(name)   ^='C' then name='FFOLD';
  if type(sangle) ^='N' then sangle=0;

  /*-- Establish viewports --*/
  np = max(down,across);
  pd = np - (across||down);
  size = int(100 / np);
  do i = 1 to across;
     px = size # ((i-1) // i) + (pd[1] # size/2);
     do j = 1 to down;
        py = 100 - (size # (j//(j-1)) + (pd[2] # size/2));
        ports = ports // shape( (px||py), 1);
     end;
  end;
  nport=nrow(ports);

  run odds(dim, table, lnames, odds);
  if ncol(dim)<3 then k=1;     * number of panels;
                 else k = dim[3];
  page = 0;                    * number of pages;
  do i=1 to k;
     r = 2#i;                  * row index, this table;
     t=table[((r-1):r),];      * current 2x2 table;

     /* construct top label for this panel */
     title='';
     if k > 1 then do;
        if vnames[,3] = " " then title = lnames[3,i];
        else title=trim(vnames[,3])+': '+lnames[3,i];
        end;

     /* standardize table to fit 100x100 square */
     run stdize(fit, t, table);

     print title;
     print fit[c=(lnames[1,]) r=(lnames[2,]) f=8.2] ' '
             t[c=(lnames[1,]) r=(lnames[2,]) f=8.0] ;

     /*-- start new page if needed --*/
     if mod(i,nport)=1 then do;
        call gstart;
        page = page+1;                * count pages;
        gname =trim(name)+char(page,1,0);
        call gopen(gname);            * name uniquely;
        end;

     /*-- set viewport --*/
     ip = 1 + mod(i-1,nport);         * viewport number;
     port = ports[ip,];               * coordinates;
     call gport(port);

     /*-- draw this panel, display if end-of page --*/
     call gpie2x2(fit, t, lnames, vnames, title,
                  np, odds[i]);
     if mod(i,nport)=0 | i=k then call gshow;
   end;
   call gclose;
finish;

start stdize(fit, t, table) global(std);
  /*-- standardize table to equal margins --*/
  if std='MARG' then do;
     config = {1 2};
     newtab = {50 50 , 50 50 };
     call ipf(fit,status,{2 2},newtab,config,t);
  end;

  /*-- standardize to largest cell in EACH table --*/
  else if std='MAX' then do;
     n = t[+,+];
     fit = (t/n)#200 ;
     fit = fit# 100/max(fit);
  end;

  /*-- standardize to largest cell in ALL tables --*/
  else if std='MAXALL' then do;
  fit = t # 100 / max(table);
  end;
  else fit = t;   /* raw counts */
finish;

start gpie2x2(tab,freq,lnames,vnames,title,np,d)
      global(sangle);
  /*-- Draw one fourfold display --*/
  t = shape(tab,1,4);     * vector of scaled frequencies;
  r = 5 * sqrt(t);        * radii of each quarter circle;

  /*-- set graphic window, font, and text height */
  call gwindow({-16 -16 120 120});
  call gset('FONT','DUPLEX');
  ht = 2.0 # max(np,2);
  call gset('HEIGHT',ht);

  /*-- set shading patterns for each quadrant */
  /* cell:[1,1] [1,2] [2,1] [2,2] */
  angle1 = { 90    0   180   270 };
  angle2 = {180   90   270     0 };
  shade  = {'L1' 'X1'  'X1'  'L1',
            'X1' 'L1'  'L1'  'X1'}[1+(d>0),];

  /*-- draw quarter circles, with color and shading */
  do i = 1 to 4;
     pat = shade[,i];
     if pat='X1' then color='BLUE';
                 else color='RED';
     call gpie(50,50, r[i], angle1[i], angle2[i],
               color, 3, pat);
     end;

  /*-- draw frame and axes --*/
  call gxaxis({0 50},100,11,1,1);
  call gyaxis({50 0},100,11,1,1);
  call ggrid({0 100}, {0 100});

  /*-- set label coordinates, angles --*/
  lx = { 50, -.5, 50, 101};
  ly = { 99, 50, -1,  50};
  ang= {  0,  0,  0,   0};
  /*-- label justification position --*/
       /*  ab  lt  bl   rt  */
  posn = {  2,  4,  8,   6};
  if vnames[,1] = " " then vl1 = '';
     else vl1= trim(vnames[,1])+': ';
  vl2='';

  /*-- are side labels rotated? --*/
  if sangle=90 then do;
     ang[{2 4}] = sangle;
     posn[ {2 4}] = {2 8};
     if vnames[,2] ^= " "
        then vl2= trim(vnames[,2])+': ';
     end;
  labels = (vl2 + lnames[2,1])//
           (vl1 + lnames[1,1])//
           (vl2 + lnames[2,2])//
           (vl1 + lnames[1,2]);
  run justify(lx,ly,labels,ang,posn,ht,xnew,ynew,len);

  /*-- write actual frequencies in the corners --*/
  cells = char(shape(freq,4,1),4,0);
  lx = {  5, 95,  5, 95};
  ly = { 94, 94,  4,  4};
  ang= {  0,  0,  0,  0};
  posn={  6,  4,  6,  4};
  run justify(lx,ly,cells,ang,posn,ht,xnew,ynew,len);

  /*-- write panel title centered above */
  if length(title)>1 then do;
     ht=1.25#ht;
     call gstrlen(len,title,ht);
     call gscript((50-len/2),112,title,,,ht);
  end;
finish;

start justify(x, y, labels, ang, pos, ht, xnew, ynew, len);
  /* Justify strings a la Annotate POSITION variable.
     x, y, labels, ang and pos are equal-length vectors.
     Returns justified coordinates (xnew, ynew)
  */

  n = nrow(x);
  call gstrlen(len,labels);
  xnew = x; ynew =  y;

  /* x and y offset factors for each position   */
  /*       1   2   3   4   5   6   7   8    9   */
  off1 = {-1 -.5   0  -1 -.5   0  -1   -.5   0  };
  off2 = { 1   1   1   0   0   0  -1.6 -1.6 -1.6};

  do i=1 to n;
     if ang[i]=0 then do;
        xnew[i] = x[i] + (off1[pos[i]]#len[i]);
        ynew[i] = y[i] + (off2[pos[i]]#ht);
     end;
     else if ang[i]=90 then do;
        ynew[i] = y[i] + (off1[pos[i]]#len[i]);
        xnew[i] = x[i] - (off2[pos[i]]#ht);
     end;
  call gscript(xnew[i], ynew[i], labels[i], ang[i]);
  end;
  len = round(len,.01);
finish;

start odds(dim, table, lnames, odds);
  /*-- Calculate odds ratios for 2x2xK table --*/
  k = dim[3];
  do i=1 to k;
     r = 2#i;
     t=table[((r-1):r),];
     odds = odds // ( t[1,1]#t[2,2] / (t[1,2]#t[2,1]) );
     end;
  rl = lnames[3,];
  odds = odds || log(odds);

  title= 'Odds (' + trim(lnames[1,1]) + '|' + trim(lnames[2,1])
         + ') / ('+ trim(lnames[1,1]) + '|' + trim(lnames[2,2]) + ')';
  reset noname;
  print title;
  print odds[r=rl c={'Odds Ratio' 'Log Odds'} format=10.4];
  reset name;

  odds = odds[,2];      *-- return log odds ratios;
finish;