Example 8.33: Merging data sets one-to-many

It's often necessary to combine data from two data sets for further analysis. Such merging can be one-to-one, many-to-one, and many-to-many. The most common form is the one-to-one match, which we cover in section 1.5.7. Today we look at a one-to-many merge.

Since the Major League baseball season started last Thursday, we'll use baseball as an example. Sean Lahman has compiled a large set of baseball data. In fact, it's large enough that it's only hosted in zipped form. The zipped comma-delimited files, with data through the 2010 season, can be downloaded here.

One file you get is the batting.csv file. This contains yearly batting data for all players. However, if you want to identify players by name, you have to use the playerid variable to link to the master.csv data set, which has names and other information. We'll add the names to each row of the batting data set. Then we can pull up players' batting data directly by name instead of with the playerid

SAS

We'll start by reading the data from the csv files using proc import (section 1.1.4).

proc import datafile="c:\ken\baseball\master.csv" 
  out=bbmaster dbms=dlm;
  delimiter=',';
  getnames=yes;
run;

proc import datafile="c:\ken\baseball\batting.csv" 
  out=bbbatting dbms=dlm;
  delimiter=',';
  getnames=yes;
run;

Then we sort on playerid in both data sets and use the merge and by statements to link them (section 1.5.7). SAS replicates each row of the master data set for every row of the batting data set with the matching by values.

proc sort data=bbmaster; by playerid; run;
proc sort data=bbbatting; by playerid; run;

data bbboth; merge bbmaster bbbatting;
  by playerid;
run;

Then we can use the data! Here, we plot the annual (regular season) RBIs of Derek Jeter and David Ortiz. Note the use of the v= option to the symbol statement (section 5.2.2) to display the players' names at the plot locations, and the offset option in the axis definition to make extra space for those names to plot.

symbol1 i=none font=swissb v="JETER" h=2 c=red;
symbol2 i=none font=swissb v="ORTIZ" h=2 c=blue;
axis1 offset = (1.5cm,1.5cm);
axis2 offset = (.5cm,);
proc gplot data=bbboth;
where (namelast="Jeter" and namefirst="Derek") or 
      (namelast="Ortiz" and namefirst="David");
plot rbi * yearid = namelast / haxis=axis1 vaxis=axis2 nolegend;
run;

The result is shown above.

R
We start here by reading the data sets. Then we use the merge() function to generate the desired dataset. As with SAS, the default behavior of the merging facility does what we need in this case.

master = read.csv("bball/Master.csv")
batting = read.csv("bball/Batting.csv")
mergebb = merge(batting,master)

To make the plot, we first make a new data set with just the information about Jeter and Ortiz. This isn't really necessary, but it does make typing slightly less awkward in later steps. Then we make an empty plot. In order to make room for the names (which are much bigger than usual plot symbols) we have to set the x and y limits manually (section 5.3.7).

jo = mergebb[
  (mergebb$nameLast == "Jeter" & mergebb$nameFirst == "Derek") | 
  (mergebb$nameLast == "Ortiz" & mergebb$nameFirst == "David"),]
plot(jo$RBI~jo$yearID, type = "n",xlim = c(1993, 2012), 
  ylim = c(-10,160), xlab = "Year", ylab = "RBI")

Then we can add the text values, using the text() function (section 5.2.11). We do this by separately pulling rows from the new jo dataset. In the reduced data set, we can specify rows using last names only.

text(jo$yearID[jo$nameLast == "Jeter"], jo$RBI[jo$nameLast == "Jeter"],
  "JETER", col="red")
text(jo$yearID[jo$nameLast == "Ortiz"], jo$RBI[jo$nameLast == "Ortiz"],
  "ORTIZ", col = "blue")

The result is seen below. David Ortiz has driven in more runs than Derek Jeter since about 2002 (Go Sox!).

Note that if space requirements prevent making a single massive dataset with many replicated rows, you can generate a lookup vector using the match() function, and use this to make the short data set. This version won't have the names in it, though.

matchlist = match(batting$playerID, master$playerID)
lastname.batting = master$nameLast[matchlist]
firstname.batting = master$nameFirst[matchlist]

jo2 = batting[
  (lastname.batting == "Jeter" & firstname.batting == "Derek") | 
  (lastname.batting == "Ortiz" & firstname.batting == "David"),]

Example 8.31: Choropleth maps

In our book, we show a simple example of a map (section 6.4.2) where we read the boundary files as data sets and use SAS and R to plot them. But both SAS and R have complex functionality for using pre-compiled map data. To demonstrate them, we'll show how to make a simple choropleth map, using US Census data available here. The file also includes data on the whole country, and by region, which may complicate use.

One feature of maps of large areas is that in representing the surface of the earth in two dimensions, they necessarily involve some distortion. Cartographers use "projections" to specify how the globe should be represented.

SAS
We first read the data from the URL mentioned above (section 1.1.6). Note that the file is rather wide, and SAS fails to read all the variable names. We then create a new variable indicating whether the state had a population above or below 3,000,000 in 2000.

filename census url 
"http://www.census.gov/popest/national/files/NST_EST2009_ALLDATA.csv"; 

proc import datafile=census out=ds dbms=csv replace;
getnames=yes;
guessingrows = 58;
run;

data cpop;
set ds;
pop3mill = (census2000pop > 3000000);
run;

SAS/GRAPH comes with many map data sets; they are installed in the maps directory. The "states" data set has the boundaries of all 50 states and Washington, D.C. But plotting Alaska and Hawaii is a hassle that's beyond our scope today, so we drop them. We also don't need especially fine resolution. By selecting density < 4 we can reduce the size of the map data set by a factor of 5.

data us48;
   set maps.states;
   if (state ne 2) and (state ne 15) and (state ne 72);
   if density < 4;
run;

The gproject procedure takes a map stored as latitude and longitude, and projects it into x and y values for plotting. The list of available projections is somewhat overwhelming, and we arbitrarily picked one available in R. SAS will choose convenient projection parameters for you (the parallels, in this case) by default, but to ensure the same results as when using R, we set them to some nice round numbers.

proc gproject data=us48 out=us48proj
  parallel1 = 35 parallel2=50;
  id state;
run;

Now we're ready to run proc gmap. Both the SAS map data set and the Census data set included a state ID number that happened to mean the same thing. This is used in the id statement and links the regions in the map data to those in the data data set. Note here that the choro statement will accept a continuous variable and will divide it into groups for you using the /nlevels= option in the choro statement, but the finer control of creating your own cutpoints is useful.

title "48 States by population";
pattern1 c=red;
pattern2 c=blue;
proc gmap map=us48proj data=cpop all;
   id state;
   choro pop3mill;
run; quit;

The results are shown above.

R
There are several R packages that support mapping. Here we use the maps and mapproj packages. The maps package includes both functions and map data sets. First, we read in the data from the Census. Then we select out the state names and the population in 2000, changing the state names to a character, from a factor object. We also make a vector with colors corresponding to the size of the state. We used the RColorBrewer package in Example 8.13, and that would be a fine tool to use here for more categories of the population size.

statepopdata <- read.csv(
 "http://www.census.gov/popest/national/files/NST_EST2009_ALLDATA.csv")
popnames = as.character(statepopdata$NAME)
pop2000 = statepopdata$CENSUS2000POP
popdich = ifelse(pop2000 < 3000000, "red", "blue")

The data() function makes the state map file available. The names of the regions in the map can be extracted using the map() function; we'll need those to link to the population data.

library(maps)
library(mapproj)
data(state)
mapnames = map("state", plot=FALSE)$names

The regions in the map file include separate names for some islands and other disconnected regions. These have the format state:region and we need to drop the colon and remainder, if present. We can do that with a complex application of the regexpr() function (section 1.4.6). First we test for the presence of a colon. If it isn't there, we return the name of the state. If it is there, we use the substr() function (section 1.4.7) to clip the colon and remainder. Finally, since the names are recorded in lower case in the map, we'll convert the population data set's name to lower also (section 1.4.8).

mapnames.state = ifelse(regexpr(":",mapnames) < 0, 
  mapnames, substr(mapnames, 1, regexpr(":",mapnames)-1))
popnames.lower = tolower(popnames)

Now we can make a color vector that copies the colors for each state into the possibly multiple regions the map includes for each state. We do this using the match() function, which returns the positions of the first argument in its second argument. Called as the index of the dichotomous indicator we made above, it will replicate the color for Michigan in the population data to both the upper and lower peninsulas of Michigan. Pretty cool! Also, since Puerto Rico, Alaska, and Hawaii don't appear in the map names vector, their population plot colors are not used. Then the choropleth is easily made with the map() function. The (many, many) projections available through the mapproj library can be found in the help for the mappropject() function.

cols= popdich[match(mapnames.state,popnames.lower)]
map("state",fill=TRUE,col=cols,proj="albers",param=c(35,50))
title("48 States by population")

The results are found below. In general, the states of the western interior are not as populous.