Showing posts with label call symput. Show all posts
Showing posts with label call symput. Show all posts

Tuesday, November 29, 2011

Example 9.16: Small multiples



Small multiples are one of the great ideas of graphics visionary Edward Tufte (e.g., in Envisioning Information). Briefly, the idea is that if many variations on a theme are presented, differences quickly become apparent. Today we offer general guidance on creating figures with small multiples.

As an example, we'll show graphics for the popularity, salary, and unemployment rates for college majors. This data was discussed here where a scatterplot graphic was presented. We draw on data and code presented there as well. The scatterplot does not show the unemployment rate, and the width of the field names and arbitrarily sized text makes it difficult to determine the popularity ranking. In contrast, the small multiples plot, demonstrated above, makes each of these features easy to see. (Click on the picture for a clearer image.)

R

The graphics options in R, particularly par("mfrow") or par("mfcol"), are well-suited to small multiples. The main tip here is to minimize the space reserved for margins and titles. In the example below, we do this with the mar, mgp, and oma settings. We'll begin by setting up the data in a process that relies heavily on the code here. (Note that the zip file referred to includes data already in the R format-- since our point today is not data management, we don't replicate the process used to make this out of the raw data.)

clean = function(x) as.numeric(gsub("\\$|, ", "", x))
clean2 = function(x) as.numeric(gsub("%", "", x))
load("table.rda")
X[,2] = clean2(X[,2])
for (i in 3:5) X[,i] = clean(X[,i])


X$cols = NA
X$cols[grep("BUSI|ACC|FINAN",X[,1])] = 1
X$cols[grep("ENGINEERING",X[,1])] = 2
X$cols[grep("STAT|COMPU",X[,1])] = 3
X$cols[grep("BIOL|CHEMI|PHYSICS|MATHEM",X[,1])] = 4
X$cols[grep("ENGLISH|HISTORY|LANG|FINE|MUSIC|PHILOS|DRAMA|LIBERAL|ARCH|THEO",X[,1])] = 5
X$cols[grep("SOCIO|PSYCH|POLI|ECONO|JOURNAL",X[,1])] = 6
X$cols[grep("COMMUN|MARKET|COMMER|MULTI|MASS|ADVERT",X[,1])] = 7
X$cols[grep("NURS|CRIM|EDU|PHYSICAL|FAMI|SOCIAL|TREAT|HOSP|HUMAN",X[,1])] = 8

This removes some funny characters and groups the fields together in a coherent manner. Then we write a function to set the par() values we need to change, and plot a pie for each row of the data set. Here a for loop is used-- we're not aware of a vectorized version of the pie() function. Colors for each pie are assigned via the colors.

sortx = X[order(X$Popularity),]

smajors = function(mt) {
old.par = par(no.readonly=TRUE)
on.exit(par(old.par))
nrows = sqrt(nrow(mt)) + (ceiling(sqrt(nrow(mt))) != sqrt(nrow(mt)))
par(mfrow=c(nrows,nrows), mar=c(1,0,0,0), oma=c(0,0,0,0), mgp=c(0,0,0))
colors = c("blue", "purple", "purple", "gray", "orange", "gray", "red", "black")
for (i in 1:nrow(mt)) {
pie(c(mt[i,2], 100 - mt[i,2]), labels=NA, radius=mt[i,4]/max(mt[,4]),
col = c("white",colors[mt[i,7]]))
mtext(substr(mt[i,1],1,19), side=1, cex=.8)
}
}

smajors(sortx[1:49,])

The resulting plot is shown above. There may be too many colors, though statistics and computing were assigned the same color as engineering. We can easily read from the plot that computing, statistics, and engineering (purple) are the largest circles, and thus the best paying. Similarly, the humanities (orange) are the worst paying. They are also not terribly popular-- the first orange pie appears in the second row. The "professions" (nursing, teaching, policing, therapy) don't pay well but are fairly popular. Most pies have roughly equal unemployment, though nursing and the professions generally are notable for near full employment. All in all, the rank, salary, unemployment percentage, and field are all clearer than in the scatterplot.

SAS

In SAS, one can use the greplay procedure to reproduce images in miniature. One has to define the size and shape allotted for each replayed plot, in a stored "template." This allows enormous control, but at the cost of some complexity. For example, one can create a scatterplot matrix using proc gplot instead of proc sgplot, as in this implementation by Michael Friendly. If you can generate your multiple images with a by statement, the coding for this is not too painful. However, in this case, it was not obvious how to change the color and radius for each pie using a by statement in proc gchart which includes pie charts, and would thus have been the obvious choice. In such cases, it may be easier to plot the figures directly using an annotate data set.

However, having demonstrated the use of annotate previously (e.g., Example 8.13), we show here an application using greplay, though the coding is a little bit involved. In outline, we use a macro to call for a pie to be made from each observation of the original data set. Then we use a template-making macro found here to generate the 7X7 grid template. Finally, we replay the pies into the grid.

/* read the data-- note that text file edited to remove spaces in
variable names */
proc import datafile =
'c:\sas-r dictionary\p\book\web\blog\majors\majors\table.txt'
out = majors;
getnames = yes;
run;

/* set up the field categories and colors */
data m2;
set majors;
cols=" "; /* make a missing character variable to hold them */
mf = majorfield; /* just a copy to save keystrokes */
/* the find() function is discussed in section 1.4.6 */
if sum(find(mf,"BUSI"), find(mf,"ACC"), find(mf,"FINAN")) ne 0 then cols = "blue";
if find(mf,"ENGINEERING") ne 0 the cols = "purple";
if sum(find(mf,"STAT"), find(mf,"COMPU"), find(mf,"SYSTEMS")) ne 0 then cols = "purple";
if sum(find(mf,"BIOL"), find(mf,"CHEMI"), find(mf,"PHYSICS"), find(mf,"MATH")) ne 0 then cols = "gray";
if sum(find(mf,"ENGL"), find(mf,"HIST"), find(mf,"FRENCH"), find(mf,"FINE"),
find(mf,"MUSIC"), find(mf,"PHIL"), find(mf,"DRAMA"), find(mf,"LIBERAL"),
find(mf,"ARCH"), find(mf,"THEO")) ne 0 then cols = "orange";
if sum(find(mf,"SOCI"), find(mf,"PSYCH"), find(mf,"POLI"), find(mf,"ECON"),
find(mf,"JOURN"), find(mf,"LIBERAL")) ne 0 then cols = "gray";
if sum(find(mf,"COMMU"), find(mf,"MARKET"), find(mf,"COMMER"), find(mf,"MULTI"),
find(mf,"MASS"), find(mf,"ADVERT")) ne 0 then cols = "red";
if sum(find(mf,"NURS"), find(mf,"CRIM"), find(mf,"EDU"), find(mf,"PHYSICAL"),
find(mf,"FAMI"), find(mf,"SOCIAL"), find(mf,"TREAT"),
find(mf,"HOSP"), find(mf,"HUMAN")) ne 0 then cols = "black";
drop MF; /* get rid of that keystroke-saver */
run;

/* order by popularity */
proc sort data = m2; by popularity; run;

The next macro takes a line number as input. A data step then reads that line from the real data set and uses call symput (section A.8.2) to extract as macro variables the median earnings used for the radius, the color, and the major name. It then produces two rows of data-- one with the unemployed percent and the other with the employed percent. We need this for the proc gchart input, as shown in the second half of the macro.

%macro smpie(obs);
data t1;
set m2 (firstobs = &obs obs = &obs);
call symput('rm', medianearnings);
call symput('color', cols);
call symput('name', strip(substr(majorfield,1,19)));
employed = "No"; percent = unemploymentpercent; output;
employed = "Yes"; percent = 100 - unemploymentpercent; output;
run;

pattern1 color= white v=solid;
pattern2 color= &color v=solid; /* only pattern2 should be needed, I think, but */
pattern3 color= &color v=solid; /* sometimes SAS required pattern3, in my trials */
title h=7pct "&name";
proc gchart data = t1 gout = kkpies;
pie percent / group = majorfield sumvar = percent value = none
noheading nogroupheading radius = %sysevalf((&rm * 45)/ 105000)
name = "PIE&obs" ;
run;quit;
%mend smpie;

Of particular note in the forgoing are the gout and name options. The former specifies a location where output plots should be stored. The latter assigns a name to this particular plot. In addition, the %sysevalf macro function is needed to perform mathematical functions on the radius variable.

Next, we write and call another macro to call the first repeatedly. Making the image small to begin with (using the goptions statement [section 5.3]) reduces quality loss when replaying as small multiples.

%macro pies;
goptions hsize=1in vsize=1in;
%do i = 1 %to 49;
%smpie(&i);
%end;
%mend;

%pies;

Finally, we can make the template to replay the images, and replay them, both using proc greplay.

/* key elements of the %makegridtemplate macro:
how many rows and columns (down and across).
name for the template.
Note that the macro is called *inside* the proc greplay statement,
which allows the user access to pro greplay statment options */
proc greplay nofs tc=work.templates;
%makeGridTemplate (across=7, down=7, ordering=LRTB,
hgap=0, vgap=0, gapT=0, gapL=0, gapr=0, name=sevensq,bordercolor=white)
quit;

/* this macro just types out text for us the sequence 1:pie1 2:pie2 ... 49:pie49 */
/* we need that to replay the figures in proc greplay */
%macro pielist;
%do i = 1 %to 49; &i:pie&i %end;;
%mend;

filename pies "c:\temp\pies2.png";
goptions hsize=7in vsize=7in gsfname=pies device=png;

/* now ready to replay the plots
The proc greplay options say what template to use and where to find it,
and where to find the input and place the output */
/* The treplay statement plays the old plots to the locations specified
in the template
proc greplay template=sevensq tc=work.templates nofs gout=kkpies igout=kkpies;
treplay %pielist;
run;
quit;

The net outcome of this is shown below. The image is pretty disappointing-- the circles are not round,and the text is pretty blurry. However, the message is as clear as in the prettier R version.

Monday, July 12, 2010

Example 8.2: Digits of Pi, redux

In example 8.1, we considered some simple tests for the randomness of the digits of Pi. Here we develop a different test and implement it.

If each digit appears in each place with equal and independent probability, then the places between recurrences of a digit should be Pr(gap = x) = .9^x * .1-- the probability the digit recurs immediately is 0.1, the probability it recurs after one intervening digit is 0.09, and so forth. Our test will compare the observed distribution of gaps to the null hypothesis described by the equation above. We'll just look at the gaps between the occurrence of the numeral 1.

R

The piby1 object, created in example 8.1, is a vector of length 10,000,000 containing the digits of Pi. We use the grep() function (section 1.4.6) to identify the places in the vector where the string "1" appears. Then we subtract the adjacent locations to find the distance between them, subtracting 1 to count the spaces.

grep("1", piby1)
gap = pos1[2:length(pos1)] - pos1[1:length(pos1) -1] -1


> head(gap)
[1] 1 33 2 8 18 25

Recollecting that Pi begins 3.141, the first value, 1, in the gap vector reflects accurately that there is a single non-1 digit between the first and second occurrences of the numeral 1. There are then 33 non-1 digits before another 1 turns up.

We next need to calculate the probabilities of each gap length under the null. We'll do that with the : shorthand for the seq() function (section 1.11.3). We'll lump all gaps larger than 88 digits together into one bin.

probgap = 0.9^(0:88) * .1
probgap[89] = 1 - sum(probgap[1:88])

How unusual was that 33-digit gap between the second and third 1?

> probgap[34]
[1] 0.003090315

Under the null, that gap should happen only about once in every 300 appearances.

Next, we lump the observed gaps similarly, then run the one-way chi-squared test.

obsgap = c(gaptable[1:88], sum(gaptable[89:length(gaptable)]))
chisq.test(obsgap, p=probgap)


Chi-squared test for given probabilities
data: obsgap
X-squared = 82.3844, df = 88, p-value = 0.6488

This test, more stringent than the simple test for equal frequency of appearance, also reveals no violation of the null.


SAS

In SAS, this operation is much more complex, at least when limited to data steps. Possibly it would be easier with proc sql.

We begin by finding the spaces between appearances of the numeral 1, using a retain statement to "remember" the number of observations since the last time a 1 was observed.

data gapout;
set test;
retain gap 0;
target=1;
if digit eq target then do;
output;
gap=0;
end;
else gap = gap + 1;
run;


Next, we tabulate the number of times each gap appeared, using proc freq to save the result in a data set (section 2.3.1).

proc freq data=gapout (firstobs=2) noprint;
tables gap / out=c1gaps;
run;


Next, we calculate the expected probabilities, simultaneously lumping the larger gaps together. Here we manually summed the cases-- this would be a somewhat more involved procedure to automate. The final if statement (section 1.5.1) prevents outputting lines with gaps larger than 88.

set c1gaps;
retain sumprob 0;
if _n_ lt 89 then do;
prob = .9**(_n_ - 1) * .1;
sumprob = sumprob + prob;
end;
else if _n_ eq 89 then do;
prob = 1 - sumprob;
count=93;
end;
if _n_ le 89;
run;

As involved as the above may seem, the real problem in SAS is to get the null probabilities into proc freq for the test. SAS isn't set up to accept a data set or variable in that space. One option would be to insert it manually, but it may be worthwhile showing how to do it via a global macro variable.

First, we use proc transpose (section 1.5.4) to make a single observation with all of the probabilities in it.

proc transpose data=c1gaps2 out=c1gaps3;
var prob;
run;


Next we make a global macro variable to store the printed values of the probabilities, and use the call symput function (section A.8.2) to put the values into the variable; the final trick is to use the of syntax (section 1.11.4) to print all of the probabilities at once.

%global explist;

data probs;
set c1gaps3;
call symput ("explist", catx(' ', of col1-col89));
run;


We're finally ready to run the test. The global macro is called in the tables statement to define the expected probabilities. The weight statement tells SAS that there are count observations in each row.

proc freq data=count1gaps2;
tables gap / chisq testp= (&explist);
weight count;
run;
The FREQ Procedure

Chi-Square 82.3844
DF 88
Pr > ChiSq 0.6488

Sample Size = 999332