/*----------------------------------------------------------------* | 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) x 2 | | 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 (see config) | | ='MAX' standardizes each table to a maximum | | cell frequency of 100. | | ='MAXALL' standardizes all tables to | | max(f[i,j,k])=100. | | config ={1 2}|{1}|{2} specifies the margins to standardize | | down number of panels down each page | | across number of panels across each page | | sangle angle for side labels (0|90) | | colors names of two colors to use for the smaller and | | larger diagonals of each 2x2 table. | | patterns names of two fill patterns. For grayscale, use | | {SOLID SOLID} and colors={GRAYC0 GRAY80}. | | alpha error rate for confidence rings on odds ratios; | | 0 to suppress | | conf type of confidence rings ='Individual'|'Joint' | | font font used for labels. | *----------------------------------------------------------------* * Author: Michael Friendly * * Citation: SAS/IML graphics for fourfold displays, Observations,* * 3rd Quarter 1994, 47-56. * * Created: 9 May 1991 19:12:12 Copyright (c) 1992 * * Revised: 8 Jun 1994 08:38:22 * * Version: 1.7 * * Changes: Added tests for homogeneity and cond. independence * * Added options for colors, patterns, font * * Added confidence rings for odds ratios * *----------------------------------------------------------------*/ goptions gunit=pct; start fourfold(dim, table, vnames, lnames) global (std, config, down, across, name, sangle, colors, patterns, alpha, conf, font ); if ncol(dim)<3 then k=1; * number of panels; else k = dim[3]; /*-- Set global defaults --*/ if type(std) ^='C' then std='MARG'; if type(config) ^='N' then config={1 2}; if type(down) ^='N' then down=min(k,2); if type(across) ^='N' then across=1; if type(name) ^='C' then name='FFOLD'; if type(sangle) ^='N' then sangle=90; if type(colors) ^='C' then colors= {RED BLUE}; if type(patterns)^='C' then patterns={L1 X1}; if type(alpha) ^='N' then alpha=.05; if type(conf) ^='C' then conf='Individual'; if type(font) ^='C' then font='DUPLEX'; /*-- 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, bounds); if k>1 then run tests(dim, table, vnames); 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 nport=1 | 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,2]); if alpha>0 then run conflim(dim, t, bounds[i,], table); if mod(i,nport)=0 | i=k then call gshow; end; call gclose; finish; start stdize(fit, t, table) global(std, config); /*-- standardize table to equal margins --*/ if std='MARG' then do; 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, colors, patterns, font); /*-- 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',font); 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 = (shape(patterns[,{1 2 2 1 2 1 1 2}],2))[1+(d>0),]; color = (shape( colors[,{1 2 2 1 2 1 1 2}],2))[1+(d>0),]; /*-- draw quarter circles, with color and shading */ do i = 1 to 4; call gpie(50,50, r[i], angle1[i], angle2[i], color[,i], 3, shade[,i]); 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,trim(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, bounds) global(alpha, conf); /*-- Calculate odds ratios for 2x2xK table --*/ if ncol(dim)<3 then do; k = 1; rl=''; end; else do; k = dim[3]; rl = lnames[3,]; end; 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]) ); se = se // sqrt( sum( 1 / (t + .01#(t=0)) ) ); end; z = log(odds)/se; p = 2#(1 - probnorm(abs(z))); odds = odds || log(odds) || se || z || p; title= 'Odds (' + trim(lnames[1,1]) + '|' + trim(lnames[2,1]) + ') / ('+ trim(lnames[1,1]) + '|' + trim(lnames[2,2]) + ')'; reset noname; print title; cl={'Odds Ratio' 'Log Odds' 'SE(Log)' 'Z' 'Pr>|Z|'}; print odds[r=rl c=cl format=10.4]; /* Find confidence intervals for log(odds) and odds */ if all(alpha=0) then return; if nrow(alpha)>1 then alpha=shape(alpha,1); cl={'Lower Odds' 'Upper Odds' 'Lower Log' 'Upper Log'}; if substr(conf,1,1)='I' then conf='Individual'; else conf='Joint'; do i=1 to ncol(alpha); if conf='Individual' then pr=alpha[i]; else pr= 1 - (1-alpha[i])##(1/k); z = probit(1 - pr/2); ci= (odds[,2] * j(1,2)) + (z # se * {-1 1}); co= exp(ci); bounds = bounds || co; ci= co || ci; print conf 'Confidence Intervals, alpha=' pr[f=6.4] ' z=' z[f=7.4 ], ci[r=rl c=cl f=10.4]; end; reset name; finish; start tests(dim,table, vnames); config = {1 1 2, /* Test homogeneity of odds ratios */ 2 3 3}; /* (no three-way association) */ call ipf(m, status, dim, table, config); chisq = sum ( (table - m)##2 / ( m + (m=0) ) ); chisq = 2 # sum ( table # log ( (table+(table=0)) / (m + (m=0)) ) ); df = dim[3] - 1; prob= 1- probchi(chisq, df); out = chisq || df || prob; test={'Homogeneity of Odds Ratios'}; print 'Test of Homogeneity of Odds Ratios (no 3-Way Association)', test chisq[f=8.3] df prob[f=8.4]; config = { 1 2, /* Conditional independence */ 3 3}; /* (Given no 3-way) */ call ipf(m, status, dim, table, config); chisq = 2 # sum ( table # log ( (table+(table=0)) / (m + (m=0)) ) ) - chisq; df = 1; prob= 1- probchi(chisq, df); *-- CMH test (assuming no 3-way); free m; k = dim[3]; do i=1 to k; r = 2#i; * row index, this table; t=table[((r-1):r),]; * current 2x2 table; n = n // t[1,1]; s = t[+,+]; m = m // t[+,1] # t[1,+] / s; v = v // (t[+,][#]) # (t[,+][#]) / (s#s#s-1); end; cmh = (abs(sum(n) - sum(m)) - .5)##2 / sum(v); chisq = chisq // cmh; df = df // df; prob = prob // (1 - probchi(cmh,1)); out = out // (chisq || df || prob); test = {'Likelihood-Ratio','Cochran-Mantel-Haenszel '}; head = 'Conditional Independence of '+trim(vnames[1]) +' and '+trim(vnames[2])+' | '+trim(vnames[3]); reset noname; print head, '(assuming Homogeneity)'; reset name; print test chisq[f=8.3] df prob[f=8.4]; finish; start conflim(dim, t, bounds, table) global(std); do l=1 to ncol(bounds); run limtab(dim, t, bounds[l], tbound); *-- standardize the fitted table the same way as data; run stdize(fit, tbound, table); s = shape(fit,1,4); * vector of scaled frequencies; r = 5 * sqrt(s); * radii of each quarter circle; angle1 = { 90 0 180 270 }; angle2 = {180 90 270 0 }; pat = 'EMPTY'; color='BLACK'; do i = 1 to 4; call gpie(50,50, r[i], angle1[i], angle2[i], color, 3, pat); end; end; finish; start limtab(dim, t, bounds, tbound) global(std); /* find 2x2 tables whose odds ratios are at the limits of the confidence interval for log(odds)=0 */ odds = bounds[1]; a = sqrt(odds)/(sqrt(odds)+1); b = 1 - a; ltab = (a || b) // (b || a); *-- construct a fitted table with given odds ratio, but same marginals as data; config = {1 2}; call ipf(tbound,status,{2 2},t,config,ltab); * print odds[f=7.3] ltab[f=6.3] t[f=5.0] tbound[f=8.2]; finish;