proc iml; /*-- Load graphics --*/ call gstart; /*--------------------*/ /*-- Define modules --*/ /*--------------------*/ /* Module : compute contours */ start contour(c,x,y,npoints,pvalues); /* This routine computes contours for a scatter plot */ /* c returns the contours as consecutive pairs of columns */ /* x and y are the x and y coordinates of the points */ /* npoints is the number of points in a contour */ /* pvalues is a column vector of contour probabilities */ /* the number of contours is controlled by the ncol(pvalue) */ xx=x||y; n=nrow(x); /* Correct for the mean */ mean=xx[+,]/n; xx=xx-mean@j(n,1,1); /* Find principal axes of ellipses */ xx=xx` *xx/n; call eigen(v,e,xx); /* Set contour levels */ c=-2*log(1-pvalues); a=sqrt(c*v[1]); b=sqrt(c*v[2]); /* Parameterize the ellipse by angle */ t=((1:npoints)-{1})#atan(1)#8/(npoints-1); s=sin(t); t=cos(t); s=s` *a; t=t` *b; /* Form contour points */ s=((e*(shape(s,1)//shape(t,1)))+mean`@j(1,npoints*ncol(c),1))`; c=shape(s,npoints); /* Returned as ncol pairs of columns for contours */ finish contour; /*-- Module : draw contour curves --*/ start gcontour(t1, t2); run contour(t12, t1, t2, 30, {.5 .8 .9}); window=(min(t12[,{1 3}],t1)||min(t12[,{2 4}],t2))// (max(t12[,{1 3}],t1)||max(t12[,{2 4}],t2)); call gwindow(window); call gdraw(t12[,1],t12[,2],,'blue'); call gdraw(t12[,3],t12[,4],,'blue'); call gdraw(t12[,5],t12[,6],,'blue'); call gpoint(t1,t2,,'red'); finish gcontour; /*-- Module : find median, quartiles for box and whisker plot --*/ start boxwhskr(x, u, q2, m, q1, l); rx=rank(x); s=x; s[rx,]=x; n=nrow(x); /*-- Median --*/ m=floor(((n+1)/2)||((n+2)/2)); m=(s[m,])[+,]/2; /*-- Compute quartiles --*/ q1=floor(((n+3)/4)||((n+6)/4)); q1=(s[q1,])[+,]/2; q2=ceil(((3*n+1)/4)||((3*n-2)/4)); q2=(s[q2,])[+,]/2; h=1.5*(q2-q1); /*-- step=1.5*(interquartile range) --*/ u=q2+h; l=q1-h; u=(u>s)[+,]; /*-- adjacent values -----------------*/ u=s[u,]; l=(l>s)[+,]; l=s[l+1,]; finish boxwhskr; /*-- Box and Whisker plot --*/ start gbxwhskr(t, ht); run boxwhskr(t, up, q2,med, q1, lo); /*---Adjust screen viewport and data window */ y=min(t)//max(t); call gwindow({0, 100} || y); mid = 50; wlen = 20; /*-- Add whiskers */ wstart=mid-(wlen/2); from=(wstart||up)//(wstart||lo); to=((wstart//wstart)+wlen)||from[,2]; /*-- Add box */ len=50; wstart=mid-(len/2); wstop=wstart+len; from=from//(wstart||q2)//(wstart||q1)// (wstart||q2)//(wstop||q2); to=to//(wstop||q2)//(wstop||q1)// (wstart||q1)//(wstop||q1); /*---Add median line */ from=from//(wstart||med); to=to//(wstop||med); /*---Attach whiskers to box */ from=from//(mid||up)//(mid||lo); to=to//(mid||q2)//(mid||q1); /*-- Draw box and whiskers */ call gdrawl(from, to,,'red'); /*---Add minimum and maximum data points */ call gpoint(mid, y ,3,'red'); /*---Label min, max, and mean */ y=med//y; s={'med' 'min' 'max'}; call gset("font","swiss"); call gset('height',13); call gscript(wstop+ht, y, char(y,5,2),,,,,'blue'); call gstrlen(len, s); call gscript(wstart-len-ht,y,s,,,,,'blue'); call gset('height'); finish gbxwhskr; /*-- Module : do scatter plot matrix --*/ start gscatmat(data, vname); call gopen('scatter'); nv=ncol(vname); if (nv=1) then nv=nrow(vname); cellwid=int(90/nv); dist=0.1*cellwid; width=cellwid-2*dist; xstart=int((90 -cellwid * nv)/2) + 5; xgrid=((0:nv)#cellwid + xstart)`; /*-- Delineate cells --*/ cell1=xgrid; cell1=cell1||(cell1[nv+1]//cell1[nv+1-(0:nv-1)]); cell2=j(nv+1, 1, xstart); cell2=cell1[,1]||cell2; call gdrawl(cell1, cell2); call gdrawl(cell1[,{2 1}], cell2[,{2 1}]); xstart = xstart + dist; ystart = xgrid[nv] + dist; /*-- Label variables ---*/ call gset("height", 5); call gset("font","swiss"); call gstrlen(len, vname); where=xgrid[1:nv] + (cellwid-len)/2; call gscript(where, 0, vname) ; len=len[nv-(0:nv-1)]; where=xgrid[1:nv] + (cellwid-len)/2; call gscript(4,where, vname[nv - (0:nv-1)],90); /*-- First viewport --*/ vp=(xstart || ystart)//((xstart || ystart) + width) ; /* Since the characters are scaled to the viewport */ /* (which is inversely porportional to the */ /* number of variables), */ /* enlarge it proportional to the number of variables */ ht=2*nv; call gset("height", ht); do i=1 to nv; do j=1 to i; call gportstk(vp); if (i=j) then run gbxwhskr(data[,i], ht); else run gcontour(data[,j], data[,i]); /*-- onto the next viewport --*/ vp[,1] = vp[,1] + cellwid; call gportpop; end; vp=(xstart // xstart + width) || (vp[,2] - cellwid); end; call gshow; finish gscatmat; /*-- Placement of text is based on the character height. */ /* The IML modules defined here assume percent as the unit of */ /* character height for device independent control. */ goptions gunit=pct; use pottery; vname={ Al, Fe, Mg, Ca, Na }; read all var vname into xyz; run gscatmat(xyz, vname); /*-- 3 x 3 scatter plot matrix --*/ quit; goptions gunit=cell; /*-- reset back to default --*/