# redoing the cdc obesity maps with ggplot2

 rebuild cdc obesity maps with ggplot

### Slideshow

• first obese rates
• second overweight rates
• and if the slide show does not work - here is the link to the
pictures

### 1 get the data

• for the years 1984 - 1997 I use read.xport() (foreign package) on the sas xpt files
• then the data sets became to large, so I used the ascii files read.fortran() and choose just a few columns
• here is a resulting example data set (2006 - I computed the bmi2 column for checking)
• 2012-09: I added the maps for 2011 since the new data were out

head(x2006)


State month day year age weight height sex htm   wkg  bmi bmigr bmirisk
1     1     5   2 2006  66    263    503   2 160 11955 4669     3       2
2     1     9  19 2006  56    290    603   1 191 13182 3632     3       2
3     1    12  12 2006  40    230    511   1 180 10455 3215     3       2
4     1     4  29 2006  38    320    603   1 191 14545 4008     3       2
5     1     4  29 2006  52    120    504   2 163  5455 2064     1       1
6     1     8   2 2006  32    165    510   2 178  7500 2372     1       1
heightcm  weightkg     bmi2
1   160.02 119.29417 46.58764
2   190.50 131.54110 36.24695
3   180.34 104.32570 32.07799
4   190.50 145.14880 39.99664
5   162.56  54.43080 20.59763
6   177.80  74.84235 23.67467


### 2 compute rates and plot the graphs

library(ggplot2)
library(scales)
library(plyr)
library(maps)

## map of the states (part of the map package)
states_map <- map_data("state")
states_map$region <- factor(states_map$region)

## got fips form here and saved it as txt; http://www.epa.gov/enviro/html/codes/state.html
fips$State.Name <- tolower(as.character(fips$State.Name))

## build the graphs

filenames <- paste("bmi",1984:2010,".rdata",sep="")
for(file in filenames){
year <- substr(file,4,7)
x <- get(paste("x",year,sep=""))

## for adding the year to the plot
testdf <- data.frame(x2=-70,y2=49,year=year)
## for the first 4 years was no bmi in the data set
## I named my computed one "bmi" so I need another "bmi2" for the loop, not very sophisticated,
## but it works

if(!("bmi2" %in% names(x))){
print(file)
x$bmi2 <- x$bmi
}

## bmi groups
x$bmi2gr <- cut(x$bmi2,breaks=c(0,25,30,300),include.lowest=T,labels=c("1","2","3"))

## count
x <- ddply(x,.(State),transform,perstate=sum(!is.na(bmi2)))
x <- ddply(x,.(State,bmi2gr),transform,perstate.gr=sum(!is.na(bmi2)))

dats <- unique(x[,c("State","bmi2gr","perstate","perstate.gr")])
dats <- dats[!is.na(dats$bmi2gr),] ## percents dats$perc <- dats$perstate.gr/dats$perstate
dats$ow <- as.numeric(dats$bmi2gr) > 1

## I just want the obese and overweight
## >= 25
dats2 <- dats[dats$ow==T,] dats2 <- ddply(dats2,.(State),summarize,perc=sum(perc)) dats2$gr <- "ow"

## >= 30
dats3 <- dats[dats$bmi2gr=="3",c("State","perc")] dats3$gr <- "obese"

dats <- rbind(dats2,dats3)

## identify the states in the data set using the region names in the map (fips coded)
dats <- merge(dats,fips[,2:3],by.x="State",by.y="FIPS.Code",all=T)
dats$gr[is.na(dats$gr)] <- "obese"
dats$State.Name <- factor(dats$State.Name)

## graph
ggplot(dats[dats$gr=="obese",],aes(map_id = State.Name)) + geom_map(aes(fill=perc),colour="black",map = states_map) + expand_limits(x = states_map$long, y = states_map\$lat) +
geom_text(data=testdf,aes(x=x2,y=y2,label=year),inherit.aes=F)

## save image
ggsave(file=paste("obese",substr(file,4,7),".png",sep=""))
}


output
[1] "bmi1984.rdata"
Saving 12.7 x 7.01 in image
[1] "bmi1985.rdata"
Saving 12.7 x 7.01 in image
[1] "bmi1986.rdata"
Saving 12.7 x 7.01 in image
[1] "bmi1987.rdata"
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image
Saving 12.7 x 7.01 in image


Date: 2012-08-04 21:24:07 CEST

Org version 7.8.02 with Emacs version 23

Validate XHTML 1.0