R-Chart

Word Cloud from Blog RSS

2011-09-02T03:37:00.000-07:00

Crazy busy - no time to blog recently. Time enough for pretty pictures based upon previous words though...(thanks http://www.wordle.net).

Ten News Stories of 2010 - and the Statistics that Made Them

2011-01-01T12:43:00.000-08:00

Significance Magazine
According to Significance Magazine (jointly published by Royal Statistical Society and the American Statistical Association) the following are the top ten stories of 2010.

1. Progress in the prevention of HIV: Public health studies result in HIV treatment advancements.
2. Drug regulation: restrictions and retractions:   Related to breast cancer and type 2 diabetes.
3. Measuring a teacher's value: LA Times graded teachers based on standards tests results.
4. Political rhetoric finds a helpmeet in statistics:   "a statistical recovery and a human recession."
5. Census of Marine Life: The first census of the world's seas completed in 2010.
6. Death of Frederick Jelinek: a pioneer in speech recognition and statistical methods of NLP.
7. The genetic key to Shangri-La: Dr. Paola Sebastiani genetics advancements related to longevity.
8. Screening saves: CT Scanning definitively associated with a reduced risk of lung cancer mortality.
9. Fat kills: Quantitative reviews in various areas of health and nutrition.
10. Words, words words: Culturomics project produces the Google Ngram Viewer.

The details are parceled out in 5 articles: Part 1 | Part 2  | Part 3 | Part 4 | Part 5

Other Stories - or my $0.02.
The following are not exactly in the same category as the listings in significance magazine - but they involve personalities and events that affect many members of the R community and have some sort of analytical/statistical significance.

World Statistics Day
I mean, I missed picking up a greeting card - but the objective of the celebration is pretty worthwhile:

building support and better understanding for official statistics among the general public and the policy-makers worldwide.

R-Bloggers
For the R community, R-Bloggers has had a banner year and provided a great deal of visibility for the R community. They are looking for sponsorship - so please consider supporting them.

U.S. Economic News
News involved the use of additional zeros tacked on to end of numbers. The recovery.org web site has been somewhat underwhelming. Edward Tufte's nomination to serve on the Recovery Independent Advisory Panel was a fascinating development. His emphasis on clear and truthful presentation of information could be a Good Thing.

New Era of Data Journalism
The World Bank has continued to provide more data on economic and social topics.
A couple of blog posts covered this, and an R package is also available to access the World Bank Data API. There has been an increased refinement in data journalism as well as controversy surrounding WikiLeaks during 2010.

Data Marketplace
InfoChimps is pioneering an online marketplace for buying and selling data. Seems that they have a plausible idea - they recently landed 1.2 million dollars in funding.

Benoit Mandelbrot
Another noteworthy death this year that was not mentioned was the loss of the "Father of Fractals" - Benoit Mandelbrot.

R-Chart: Year End Wrap Up

2010-12-31T14:17:00.000-08:00

Thanks to everyone who visited and commented here at R-Chart over the last year! Blogging has forced me to crystallize my thoughts and I hope others have benefited a bit from these meanderings. It it great to interact with the knowledgeable, educated and friendly folks in the R community.

I make no claims to be an expert or authority on statistics, visualization, design or any of the myriad of other topics touched on over the past year. I appreciate all who have provided encouragement, suggestions and corrections. Unlike many of you more scientifically minded types who meticulously verify all conclusions before speaking, I tend to throw ideas out in the blog and make adjustments and corrections based upon feedback. This is really one of the great values of blogging - and so again, thank you for your responsiveness. It was unexpected and very helpful.

Lessons Learned
In case you blog or are thinking of blogging, I thought you might be interested in how things have worked here at R-Chart to this point.

Make Good Titles
It was interesting to find out which items were of most interest (based upon the number of hits per page). A great deal seems to be based upon the headline to the blog - never underestimate the value of a well-constructed-sound-byte of a title. This often dictates the future of a posting. Bad title = no response. I really never gave much thought to how important it is to construct a meaningful, attention grabbing title.

Blog Promotion
Promotion of each article also took more time than I expected. Tal over at R-Bloggers really does the R community a service - bloggers who sign up have content aggregated automatically. If you want to draw additional readers you have to do a certain amount of footwork yourself. I get about 15% of total traffic to the site from search engines - which is kind of low. Most of the generic sites that I submitted the blog to didn't send any traffic. Content that was of specific interest to a given community ended up resulting in the most traffic.

The top sites that have sent traffic this way are shown below.

www.reddit.com    15,218
www.google.com 8,932
news.ycombinator.com 7,211
www.r-bloggers.com 4,885
www.dzone.com 3,682
habrahabr.ru    1,167 (Hi to friends in Russia for this - the highest ranking non-English site)
twitter.com      689
www.google.co.in      565
www.google.co.uk      531
www.rubyflow.com      470

R is International
I was really amazed at the international response - folks from 164 countries around the world hit the blog since its inception. Germany was the top non-English site in total visits and France was also well represented.

This probably is of no surprise to many - R has been widely used in academic research and there are a relatively small number of highly specialized professionals around the world using R. It's obvious that the web reaches everywhere - it is not obvious who will end up visiting a given site.

Interest as Indicated by Traffic
A few other numbers of note:

96,928 R-Chart Pageviews all time history as of 01/31/2010.

620 Downloads of the free R-Chart iPhone application

237 Total days blogging at blogspot (as.Date('2010-12-31') - as.Date('2010-05-08'))

195 Days blog has lived at r-chart.com (as.Date('2010-12-31') - as.Date('2010-06-19'))

158 Comments on this blog

Advertising

Apologies to folks who are put off by the advertising. I had a goal to dip into this area a bit to come to offset costs and maybe buy a book or two. This may happen eventually...

$ 42.89 AdSense Revenue

$ 13.46 Advertising Revenue through Amazon affiliates

Again - thanks to all - and have a Happy New Year

Google AI Challenge: Scores/Rank by Language

2010-12-08T04:00:00.000-08:00

A quick follow up to the previous post: about the the scores in the 2010 Google AI competition relative to programming language. The chart above makes each language visible and discrete - and the scales are the same.

library(ggplot2)

df<- read.csv('googleAI2010.csv',sep=';',header=FALSE)

df$V7 <- NULL

names(df)<- c('rank', 'username','country','organization','language','elo_score')

ggplot(data=df, aes(x=rank, y=elo_score, color=language)) +

+ geom_point(size=1) +

+ facet_wrap(~ language) + opts(title='Google AI 2010: Score by Rank for each Language')

It is based upon a simple comparison of rank and score.

df<- read.csv('googleAI2010.csv',sep=';',header=FALSE)
df$V7 <- NULL
names(df)<- c('rank', 'username','country','organization','language','elo_score')

ggplot(data=df, aes(x=rank, y=elo_score)) + geom_point(size=1) + opts(title='Google AI Score by Rank')

Another approach to viewing this information is a histogram by score (which ignores rank). With a binwidth of 100 (and ignoring low scores of people who signed up but who dropped out relatively early) a (nearly) bimodal distribution appears.

qplot(data=df, x=elo_score, geom='histogram', binwidth=100)

Any ideas about why this is not normal? Is there some aspect of ELO scoring that leads to this shape? Or are there different types of programmers represented?

This can be broken down by language. To avoid difficulty distinguishing colors, the rainbow palette is used and a few languages are not reported (since they were not highly represented in the competition).

library(sqldf)

df2=sqldf("select * from df where language not in ('Groovy','Scala','Go','OCaml')")
df2$language=factor(df2$language)
qplot(data=df2, x=elo_score, fill=language, geom='histogram', binwidth=100) + scale_fill_manual(values=rainbow(12))

As mentioned in the previous post, the data is available at GitHub - feel free to post some of your own visualizations of this data.

Google AI Challenge: Languages Used by the Best Programmers

2010-12-02T16:18:00.000-08:00

The Google AI Challenge recently wrapped up with a Lisp developer from Hungary as the winner. The competition challenges contestants to create bots that push the limits of AI and game theory. These bots compete against one another, and a complete ranking of competitors is available. The big story today is that the winner (Gábor Melis) used Lisp to beat out over 4000 other contestants around the world using a host of different programming languages.

Paul Graham has stated that Java was designed for "average" programmers while other languages (like Lisp) are for good programmers. The fact that the winner of the competition wrote in Lisp seems to support this assertion. Or should we see Mr. Melis as an anomaly who happened to use Lisp for this task?

Programming Languages Usage

Java, C++, Python and C# were heavily used overall.

language count(*)

1 Java 1634
2 C++ 1232
3 Python 948
4 C# 485
5 PHP 80
6 Ruby 55
7 Haskell 51
8 Perl 42
9 Lisp 33
10 Javascript 19
11 C 18
12 OCaml 12
13 Go 6
14 Scala 4
15 Groovy 1

In the Top 200

language count(*)

1 Java 70

2 C++ 64

3 Python 34

4 C# 17

5 C 4

6 Haskell 3

7 PHP 3

8 Ruby 2

9 Javascript 1

10 Lisp 1

11 OCaml 1

Top 100

1 Java 33
2 C++ 32
3 Python 20
4 C# 9
5 C 3
6 Haskell 1
7 Lisp 1
8 OCaml 1

Top 10

language count(*)

1 Java 4

2 C++ 3

3 C# 2

4 Lisp 1

The plot above is a bit difficult to discern due to the number of languages represented (and similarity in colors). So here is a breakdown by language.

Lisp does appear to be skewed towards higher ranking. But even more striking are the C hippies:

The functional crowd represented with Haskell also ranked on the higher end:

How about Java? There is a trend towards the average - but a significantly larger number of entrants used Java. It also is a language taught in many colleges, and might reflect greater student participation in these languages (although MIT did focus on Lisp back in the day...).

How about representatives from the Microsoft? Einstein and Elvis showed up - Mort was not interested.

I can post charts of other languages if anyone asks - otherwise, download the files for yourself and draw your own conclusions. And congratulations to

Gábor Melis - I am again feeling the inspiration to delve into the mysteries of Lisp and meander among mountains of parenthesis...

Methodology Used
No need to proceed further unless you are interested in how the results listed above were derived.

Basically, I used Ruby to scrape the results from the Google AI Rankings site. The results were read into Ruby, and ggplot2 and sqldf libraries were used to analyze the results.

Get the Data into R
So to find out more...I whipped up a ruby script to create a delimited file from the 47 page listing online. (Feel free to get these from their GitHub location and do some additional validation/analysis of your own). Read this file into R:

df<- read.csv('googleAI2010.csv',sep=';',header=FALSE)
df$V7 <- NULL
names(df)<- c('rank', 'username','country','organization','language','elo_score')

Sanity Check
Most of this work can be done in idiomatic R (which has some significant Lisp influences) - which might be a better way to honor the winner. However, I find myself using sqlite more and more these days - particularly in mobile development. So I used the sqldf library which uses this database behind the scenes.

Country rankings are available online, and the following emulates these results. Specifically, the number of entrants in the top 200 ranked contestants from each country can be derived as follows:

library('sqldf')

top200=df[df$rank <= 200,]

sqldf('select country, count(*) from top200 group by country order by 2 desc')

Organization rankings are similar, representing the top organizations within the top 100. There are some anomalies here, the highest ranking "Other" is not shown in the online version for obvious reasons, and only the most of these have only one entrant in the top 100 an are listed in an arbitrary manner. However, the results are otherwise the same in R.

top100=df[df$rank <= 100,]
sqldf('select organization, count(*) from top100 group by organization order by 2 desc')

R Code
The following are additional snippets of R code used to generate the results above.

# Language Usage

sqldf('select language, count(*) from df group by language order by 2 desc')

sqldf('select language, count(*) from top200 group by language order by 2 desc')
sqldf('select language, count(*) from top100 group by language order by 2 desc')

top10=df[df$rank <= 10,]
sqldf('select language, count(*) from top10 group by language order by 2 desc')

If you fiddle enough with the bucket size for histograms, you might be able to draw some conclusions... but the density plot seemed like a nicer option.

library('ggplot2')

# Substitute your favorite language of those available for Lisp below
qplot(data=df[df$language=='Lisp',], x=rank, geom='histogram', binwidth=1000) + opts(title='Lisp')

# The density plot at the top of this posting:

ggplot(data=df, aes(rank, fill=language)) +
geom_density(alpha = 0.2) +

xlim(0,5000) +

opts(title='2010 Google AI Challenge Rankings')

ggsave('program_language_density_plot.png')

# Breakdown by language:

ggplot(data=df[df$language=='Scala',], aes(rank, fill=language)) + geom_density(alpha = 0.2) + xlim(0,5000) + opts(title='Scala')

Update: I have been keeping up with the comments - and sketched out some other ways of looking at the data in another post.

Mortgage Calculator (and Amortization Charts) with R

2010-11-10T16:24:00.000-08:00

Mortgage rates have been at historic lows recently. The rates are posted various places online along with simple mortgage calculators. Such calculators illustrate the payment schedule for a mortgage based upon selected terms. But with less than a dozen lines of R code, you can do a far more sophisticated analysis.

Mortgage Calculation Function
Rather than reinvent the wheel, you can work with this nice R function by Thomas Girke (Associate Professor of Bioinformatics over at UC Riverside). At the R prompt, you can grab it from its home online by calling source:

source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/mortgage.R")

This loads the function and outputs a helpful description of the function:

The monthly mortgage payments and amortization rates can be calculted with the mortgage() function like this:

   mortgage(P=500000, I=6, L=30, amort=T, plotData=T)
   P = principal (loan amount)
   I = annual interest rate
   L = length of the loan in years

So keep in mind that there is a huge amount of R code available online:

are just the beginning.

Instant R Graphical User Interfaces

Rather than simply calling the function directly, you can quickly construct a GUI input widget using the fgui library.

library(fgui)
gui(mortgage)

With this trivial invocation, a window pops up.

Not terribly fancy, but about the simplest way you will ever be able to construct a GUI! In this case a mortgage amount of $90,000 for 10 years at 3.75% is illustrated.

After entering these values, click OK to actually call the function. This results in a good deal of interesting output. Close the pop up window and look at the R Console:

Monthly payment: $900.5512 (stored in monthPay)

Total cost: $108066.1

As indicated in this message, an R object named monthPay contains the amount of the monthly payment and can be used in subsequent R commands and calculations. You also are greeted with a graph illustrating annual interest and payments as a stacked bar chart.

Plenty of useful information! But that's not all...

Beyond the Basics
You might have noticed a number of messages regarding data stored in R objects. This is where the power of R exceeds that of any standard mortgage calculator. These objects can serve as input to other calculations or charting operations.

The aDFmonth object contains amortization data for each month, while aDFyear contains annual information. In the following example, a new data frame is created from the monthly data that does not include the amortization information and plot it using ggplot2. (The amortization data is a significantly different scale and better viewed independently).

library(ggplot2)
DF=melt(aDFmonth[-1], id.vars='Year')

ggplot(DF, aes(x=Year,y=value, group=variable)) + geom_line() + facet_wrap(~ variable, ncol=1)

You can quickly manipulate the data frame to view amortization information instead. Use the exact same ggplot call (though the facet_wrap is removed below as unnecessary for a single variable) to create a chart scaled to fit the values relevant to the amortization.

DF=melt(aDFmonth[c(1,5)], id.vars='Year')

ggplot(DF, aes(x=Year,y=value, group=variable))+ geom_line()

The limits of calculations and visualizations available in a web calculator or Excel are reached pretty quickly. R provides the means to create relatively full featured solutions in only a few lines of code.

Don't be a Turkey

2010-11-09T16:50:00.000-08:00

'Indeed, I am moving on: my new project is about methods on how to domesticate the unknown, exploit randomness, figure out how to live in a world we don't understand very well. While most human thought (particularly since the enlightenment) has focused us on how to turn knowledge into decisions, my new mission is to build methods to turn lack of information, lack of understanding, and lack of "knowledge" into decisions—how, as we will see, not to be a "turkey".'

- Nassim Nicholas Taleb

With thanksgiving on the way, an economic lesson provided by a turkey's statistical department seems appropriate. Our turkey - let's call him auRthur - like most turkeys has a statistical department at his disposal. His department is in fact tracking an index - the Turkey Welfare Index which is a reflection of how much the human race cares about auRthur. Notice the relatively positive trend... until Thanksgiving Day...

Evidently, our auRthur's statistical department utilized a model that had some flaws - "past performance is not necessarily a predictor of future returns". This is because the harvesting of the turkey is a "rare event." Rare (unprecedented) events are difficult to predict. The story is not terribly amusing to turkeys to begin with - but becomes less amusing to humans when understood as a metaphor of the financial meltdown and statistical modeling in use by banking institutions. Essentially, banking institutions assumed a huge amount of risk because a catastrophic meltdown was simply outside the realm of consideration. It was not represented in most of the models in use.

A great and vivid illustration. See Nassim Nicholas Taleb's essay where this chart and illustration originally appeared at edge.org. This article discusses the limits of statistical thinking and is a good springboard to other writings by Taleb - who was a practitioner of risk as he ran a hedge fund for a number of years and saw many of the practices in the financial industry up close and personal.

The chart above was created using R and ggplot2. The data frame named DF was populated with data related to the Turkey Welfare Index.

> DF
   TWI Day color
1 14 1 black
2 15 2 black
3 16 3 black
4 17 4 black
5 18 5 black
6 19 6 black
7 20 7 black
8 -100 8 red

UPADTE: This can be entered in a few different ways. One is through a grid (which requires that you specify the Day as a factor).

  DF=edit(data.frame())
  DF$Day=factor(DF$Day)

Plotted using ggplot2:

library(ggplot2)

ggplot(data=DF, aes(x=Day, y=TWI, fill=color)) +
  geom_bar() +
  scale_fill_manual(value= c("black", "red")) +
  theme_bw() + scale_x_discrete(breaks = NA) +
  opts(legend.position='none', axis.title.x=theme_blank(),
   axis.title.y=theme_blank(),
   title='Turkey Welfare Index')

This included a couple of somewhat unusual settings to shut off labels and axes that results in the simple "plain" appearance you see above.

So - Happy Thanksgiving - understand statistics and don't be a turkey...

iPhone App Store Acceptance Time / Download Results

2010-11-03T17:16:00.000-07:00

Complaints about the iPhone App Store are not uncommon among developers. The submission process is frustrating at best, you can expect arbitrary rejections, and Apple's policies have not always been particularly open or welcoming. If you make it through the process and get an app accepted, it can be essentially buried where it will remain unused unless you dedicate significant energy to marketing it.

And so I figured I would do a small scale experiment to check App Store response in these areas (and App Store user behavior)....and post the results here in the hopes that developers planning to write iPhone apps would benefit from my experience.

Cheaper by the (Half) Dozen
Having a limited amount of time available, and appreciating the value of quick iterations in development or in any business process, I limited myself to submitting a few apps with similar functionality and different target audiences. What follows is my findings (with a bit of R code used for analysis).

Six free iPhone apps were submitted for publication during a 21 day period between July 23, 2010 and August 13, 2010. Four apps were accepted and two were rejected.

The four apps that were accepted:

R-Chart (R Programming Community News) has been available 96 days and averaged 4.44 downloads per day. (Download it if you want to keep up with this blog and/or R community news).
FRB (U.S. Federal Reserve Board News) has been available 82 days and averaged 4.15 downloads per day.
Duq News (Duquesne University News) has been available 92 days and averaged 2.52 downloads per day.
Visit the Lehigh Valley (Lehigh Valley Tourism Info) has been available 93 days and averaged1.31 downloads per day.

Two apps were rejected

DeSales U (DeSales University News)
Blender Buzz (Blender Software Blog/News )

App Functionality and Subject Area

All apps contained essentially equivalent functionality, but differed by subject area and graphic and styling qualities. The "R-Chart" and "Blender Buzz" apps reference resources of interest to software users and were intended to promote this blog and the Blender Buzz blog. They are topical and not limited by locale or institution. "Duq News" and "DeSales U" provide news from Duquesne and DeSales Universities. "Visit the Lehigh Valley" provides information about places and events for visitors to Eastern Pennsylvania. The "FRB" app provides latest publicly available news from the U.S. Federal Reserve Board.

Review Process and Time

It took between 7 and 9 days for the App Store to review an app and either accept or reject it. It does not appear that subject area of the app contributed to its acceptance. Both R-Chart and Blender Buzz were directed at programming communities - one was accepted and the other was rejected. Likewise, one of the University apps was accepted, the other was rejected. Each app was submitted once only. None were resubmitted after initial rejection.

Download Counts

A total of 1120 have been downloaded. As noted above, the result is between 1 and 4 downloads per day.

I have not expended much effort in promotion. I tweeted about the R-Chart App and mentioned it here on this blog but never promoted any of the others. So the downloads were the result of folks searching for an app that was of interest to them.

Conclusion...
So getting an app accepted by the app store is not an insurmountable process - but does require time and planning. It is not an activity that will just take care of itself. And because of inconsistency in the process, you would do better to allot a bit of extra time for the app store process. And unless your app fills a rather unique niche, you will need to do marketing in the same way that you would for a web site or any other resource.

R Code used in the Analysis is Below

> # Read in the data

> df = read.csv('app_stats.txt')

> df

App Downloads Submitted_Date Response_Date Response

1 R-Chart 426 2010-07-23 2010-07-30 Accepted

2 FRB 340 2010-08-04 2010-08-13 Accepted

3 Duq News 232 2010-07-26 2010-08-03 Accepted

4 Visit the Lehigh Valley 122 2010-07-25 2010-08-02 Accepted

5 DeSales U 0 2010-08-02 2010-08-09 Rejected

6 Blender Buzz 0 2010-07-26 2010-08-03 Rejected

> # Get the total downloads
> sum(df$Downloads)
>
> Do a plot of downloads

> ggplot(data=df, aes(x=App, y=Downloads, fill=Response))

+ geom_bar() + coord_flip() + theme_bw()
> ggsave('app_downloads.png')
>
> ggplot(data=df, aes(x=App, y=as.numeric(Processing_Time), fill=Response)) + geom_bar() + coord_flip() + theme_bw()
> ggsave('processing_time.png')
>

> # Cast the date columns as such

> df$Submitted_Date=as.Date(df$Submitted_Date)

> df$Response_Date=as.Date(df$Response_Date)

> # Find out total dates each application was on the market

> as.Date('2010-11-03') - min(df$Submitted_Date)

> # Determine total number of days apps were being submitted for review

> max(df$Response_Date) - min(df$Submitted_Date)

> # Processing time at the app store

> df$Processing_Time = df$Response_Date - df$Submitted_Date
>
> # Time each app has been available for download

> df$days_available = as.Date('2010-11-03') - df$Response_Date

> # Downloads per day

> df$downloads_per_day = df$Downloads/as.numeric(df$days_available)

> # Limit view to selected columns

> df[c(1,7,8)]

Hadley on a Postage Stamp?

2010-10-18T19:20:00.000-07:00

Yet another reason to check out the ggplot2 wiki!

Cramer's Stock Pick Recommendations Analyzed (Part II)

2010-10-18T15:12:00.000-07:00

This is the second post (previous one here) that provides an analysis of Cramer's stock recommendations based upon the Mad Money Stock Screener as of 10/15/2010.

Recommendations by Segment
As mentioned in the previous post, recommendations are referenced either by a number below or by name.
# Description
5 Buy
4 Positive
3 Hold
2 Negative
1 Sell

Referencing calls by number is used to provide average scores and to plot results. Also as noted previously, not all segments of the show are available through the Stock Screener.

Calls made during the Lightning Round and Mail Bag have an average of less than 4 (positive). Calls driven by questions by the audience have a lower average, while those in interviews are higher.

Segment Average Call
Discussed 4.387156
Featured 4.351351
Interview 4.895954
Lighting Round 3.764465
Mail Bag 3.605691

Calls Trail the Market

This is a bit more speculative a representation. The S & P and Down Jones Industrial Average are used in this comparison. The call (on a scale of 1 to 5) is multiplied by a factor so that a smoothed condition mean line is generated on the chart. The factor is arbitrary, it just makes the line fit on the chart in a reasonable location, so the directionality of the call line is relevant - not the degree.

The S & P appears as the red line and the average of the calls for each day (times a factor of 3000) appears in blue.

The DJI is multiplied by a factor of 3000 in the chart below.

Although it is not completely clear, it appears that stock picks tend to trail the market movement. For this reason they generally sound plausible. Calls tend to be more pessimistic at the time the market has been moving down.

5 of Cramer's Favorites
A few of the stocks that Cramer has recommended in the last year qualify in a special way as "favorites." They are the stocks that appear the most times in the data with only a buy recommendation.

Wynn Resorts (WYNN) had a range of $27.00 between the lowest and highest 19 buy recommendations.

Weatherford Int'l (WFT)

had a range of

$6.18 between its 16 buy recommendations.

NVIDIA (NVDA) had range of 9.05 between its 15 buy recommendations

Cypress Semiconductor (CY) had a range of 4.58 between its 13 buy recommendations.

Teva Pharmaceutical (TEVA) had a 17.03range in its 12 buy recommendations.

Conclusion
We are bombarded with predictions and promises regularly in the news and media. Over time, I have grown more and more suspicious of the ability of individuals to consistently predict stock prices by simply having a superficial knowledge of current market motions and a general awareness of current financial news.

I often think about verifying the claims, but often don't have the time... figured I would at least take a cursory look. I hope this sort of thing becomes more common to keep the media honest.
At least Jim Cramer's show is cast somewhat in the realm of entertainment. Again, I'll refer you to Bill Alpert of Barrons who has done more extensive analysis and reporting on Cramer's recommendations.

2010 ggplot2 Case Study Competition Winners

2010-10-16T15:56:00.000-07:00

The 2010 ggplot2 Case Study Competition Winners have been announced! Congratulations to the winners!

Grand Prize: David Kahle, Rice University
Finalist: Michael Lavine, UMass Amherst
Finalist: Claudia Beleites, TU Dresden & Uni. Trieste

Check out their entries to get a glimpse of what is possible with R and ggplot2.

Cramer's Stock Pick Recommendations Analyzed

2010-10-16T13:07:00.000-07:00

Bill Alpert of Barron’s has demonstrated the use of R in financial journalism as he criticized the performance of Jim Cramer’s stock picks. Patrick Burns was an advisor for the analysis done in the article. R was an important tool that allowed them to do their research as indicated by Burns in his article and by Alpert in his presentation and in an article published in R News. Their specific critique were centered around obtaining objective verification of claims by CNBC that by following Cramer's advice, one could beat the S & P 500 index. I recommend you to these sources if you are interested in a more comprehensive analysis.

The following is is based upon data available on the Mad Money Stock Screener as of 10/15/2010.

Mad Money Stock Screener

Since Cramer's Mad Money Stock Screener is available on line there is at least an "unofficial" group of Cramer's recommendations available to analyze. It is apparent that the data is not complete, as several program segments in the drop down are not represented. Selecting any of the following segments results in no records returned:

Caller's Stock
Game Plan
Sudden Death

In addition, the date provided includes month and day (but not year information). It appears that data from one calendar year is available through the web site. This was born out by plotting the price at the time of the recommendation on charts for individual equities.

Always A Bull Market Somewhere

According to his profile, "Jim Cramer believes that there is always a bull market somewhere, and he wants to help you find it". His optimistic, entertaining and confident approach that he exudes on screen is reflected in his history of stock picks. Recommendations are indicated either by a number below or by name:

# Description

5 Buy

4 Positive

3 Hold

2 Negative

1 Sell

Cramer is functioning as an entertainer with financial knowledge. The show is geared towards providing action oriented advice (buy/sell) and tends to be skewed towards positive actions. This fits with his profile description - if there is always a bull market, there is always something to buy.

Analysis of All Recommendations

The following jitter demonstrates that a 5 (Buy Recommendation) is most often given, and the smoothing indicates that Cramer is generally positive in his ratings.

The program segment in use might also shed some light on the recommendations given. As mentioned earlier, not every program segment is represented in the data.

The specific totals represented above:

Buy Hold Negative Positive Sell

2336 7 229 422 559

The vast majority of the time, Cramer gives a buy recommendation. The second most often provided recommendation is to sell. A clear "Buy" or "Sell" is certainly more entertaining to hear than a "Meh... hold." None of this in and of itself means that Cramer's ratings are bad or inferior to other sources. It simply serves to illustrate that the program is geared towards entertainment. At best, one might hope that only the clear winners are topics of conversation on the program. However, further analysis at least calls this into question.

Individual Stock Recommendations

Not every company has an equal number of recommendations. The top 5 (in terms of total recommendations made) are Apple, Citigroup, Intel, Bank of America and Ford Motor Company. In order to get a sense of the a given recommendation in the scope of wider market history, the data from the stock screener can be superimposed on a stock chart.

Apple

Apple was given a recommendation on 90 occasions (more than twice times as many than the next most popular companies recommended). The average recommendation was 4.933333 and Cramer recommended Buy 84 times and Positive 6 times. The clear upward trend in Apple's price justifies an optimistic view in recent history. None of the remaining stocks in the top 5 recommended fits this pattern though.

Citigroup

Citigroup was given a recommendation on 40 occasions and had an average recommendation of 4.850000, Cramer recommended Buy 36 times and Positive 3 times and Negative 1 time. The negative rating was on the 12/08/2009 Lighting Round.

Intel

Intel was tied for second place with Citigroup with 40 recommendations. It had an average recommendation of 4.900000. He recommended Buy on 36 occasions and Positive on 4 occasions.

Bank of America

Bank of America had an average of 4.820513 in the 39 times it was recommended. It was given a Buy recommendation 34 times, a Positive 4 times and a negative once (During the 11/24/2009 Lighting Round).

Ford Motor

Ford Motor Company was given an average of 4.777778 in the 36 times it was recommended (a Buy 28 times and a Positive 6 times).

A Negative Example

One additional stock that is of interest is British Petroleum (BP), which had an rough ride this year due to the Deepwater Horizon Oil Spill.

Cramer issued a Sell recommendation on 10 occasions, a Positive on 1 occasion and a Buy on 6 Occasions (a total of 17 recommendations).

Conclusion

I'll leave it to you to draw your own conclusions about how to interpret recommendations given by Jim Cramer. For myself, I find him entertaining and well informed about financial news for a wide range of stocks. He also has the track record as a successful hedge fund manager over the course of a number of years. However, I am skeptical about the ability of analysts to consistently predict the direction of the market.

Depending upon the reception of this post, I can provide additional information about the methods used to obtain the data and create the charts above and also show how other stocks recommended on Mad Money have performed. Let me know if you have any interest in the comments.

Benoit Mandelbrot (the Father of Fractals) dead at 85

2010-10-16T06:46:00.000-07:00

Nicholas Taleb and Recent tweets indicate that Benoit Mandelbrot has died at age 85. Mandelbrot was a French and American mathematician, best known as the father of fractal geometry. His official biography at Yale does not yet reflect this (it was last updated in March 2010) and his Wikipedia page is protected from editing for the moment.

UPDATE: The New York Times has confirmed this.

He is known for changing the way researchers perceive and characterize the phenomenon of natural growth and for the implications for his work for scientists and mathematicians. But fractals have become iconic to the masses as expressive of the idea that beauty and mathematics are inextricably related. As I mentioned elsewhere my brother and I were inspired at an early age by the wonderful designs that could be created by simple mathematical equations.

He will be missed.

Belgian Astronomers and Exercise Machines

2010-10-06T08:45:00.000-07:00

In the twisting paths of human discovery, you never quite know what intellectual enterprise is going to result in a world changing discovery. For instance, the mathematical notion of expected value did not grow up in a sterile, academic environment. In 1654 Blaise Pascal was approached by Chevalier de Méré who was interested in gambling problems. Pascal corresponded with Fermat and thus the mathematical theory of probabilities was born.

In recent days reports on economic upheaval have often cast financial industries as institutions based upon greed and power that contribute nothing of value to society. Defenders of the free market are quick to respond with the immediate economic benefits provided by such institutions as they mitigate risk and serve as "middle men" in modern financial markets. What is seldom considered is that discoveries in one area often find application in a separate area of life that was never considered during the initial investigation. And so one day, perhaps Wall Street calculations might be put to non-financial use that benefit mankind in other ways. There is historical precedent. For instance, interests in the insurance industry served to popularize and apply a 19th century Belgian calculation in a manner that is now used on modern exercise machines.

A Belgian Astronomer: Adolphe Quetelet
Adolphe Quetelet (1796–1874) was a Belgian mathematician, astronomer and statistician. While studying astronomical activities in Paris he interacted with Joseph Fourier (1768–1830), Siméon Poisson (1781–1840) and Pierre Laplace (1749–1827). He went on to put his new found appreciation of probability to practical use in the study of the human body (a subject he had initially approached as a painter and sculptor). One calculation he created, dubbed the Quetelet Index, is a number that expresses a relationship between a person's height and weight. Quetelet was not specifically interested in the use of his index for health purposes, but simply for defining the characteristics of "normal" or "average" man.

The Financial Industry
In the mid 20th century actuaries observed increased mortality in overweight policyholders. And so in an effort to construct more accurate mortality tables the relationship between weight and cardiovascular disease became the subject of epidemiological studies. Weight tables were first used to predict life expectancy as far back as 1913. But tables of ideal or desirable weight were developed by the Metropolitan Life Insurance Company in the 1940's. In the 1960s, a small group insurance industry experts began to use the Quetelet Index. But it remained for a an American scientist to perform a comparative study of available indices and rename the Quetelet Index to the form that we know it today where it has become a subject related to health and nutrition.

An American Oceanographer, Biologist, and Physiologist
Ancel Benjamin Keys was a scientist who wrote an article for the July 1972 "Journal of Chronic Diseases" that coined the phrase "body mass index" or BMI as a modern designation for the Quetelet Index. Interestingly enough, Keys early studies culminated in a B.A. in economics and political science in 1925. His first Ph.D. was in oceanography and biology but his later work was related to his second Ph.D. focused on physiology. He is best known for two dietary contributions - the K-Ration and the Mediterranean Diet.

Keys (and others today) considered the tendency in the insurance industry to equate relative body weight with excess risk of death to be somewhat simplistic. There is worldwide variation according to diet and physical activity habits. In most industrial countries people in the middle range of body weight are healthier than those at an extreme.

Because of these types of concerns, BMI is often considered along with other concerns that can indicate potential health risks. Specifically:

A BMI in the overweight category along with certain diseases
A BMI of less than 25 and a waist size above the standard (35 for women or 40 for men)

The actual BMI ranges considered healthy or at risk are still being debated. In 1998, the U.S. National Institutes of Health changed U.S. definition of normal from 27.8 to 25 to conform to World Health Organization Standards. In addition, other countries in the world are encouraging the upper limit for BMI to be even lower than 25.

BMI and R
One’s optimal weight can be derived using the BMI and height as follows:

optimal_weight = function (height, bmi){
round((height**2 * bmi) / 703)
}

A grid, similar to one found here and the chart above can be constructed a script found at GitHub.

For example:

> create_bmi_dataframe(bmi_end=30)
19 20 21 22 23 24 25 26 27 28 29 30
60 97 102 108 113 118 123 128 133 138 143 149 154
61 101 106 111 116 122 127 132 138 143 148 153 159
62 104 109 115 120 126 131 137 142 148 153 159 164
63 107 113 119 124 130 135 141 147 152 158 164 169
64 111 117 122 128 134 140 146 151 157 163 169 175
65 114 120 126 132 138 144 150 156 162 168 174 180
66 118 124 130 136 143 149 155 161 167 173 180 186
67 121 128 134 140 147 153 160 166 172 179 185 192
68 125 132 138 145 151 158 164 171 178 184 191 197
69 129 135 142 149 156 163 169 176 183 190 196 203
70 132 139 146 153 160 167 174 181 188 195 202 209
71 136 143 151 158 165 172 179 186 194 201 208 215
72 140 147 155 162 170 177 184 192 199 206 214 221
73 144 152 159 167 174 182 190 197 205 212 220 227
74 148 156 164 171 179 187 195 203 210 218 226 234
75 152 160 168 176 184 192 200 208 216 224 232 240
76 156 164 173 181 189 197 205 214 222 230 238 246
77 160 169 177 186 194 202 211 219 228 236 245 253
78 164 173 182 190 199 208 216 225 234 242 251 260
79 169 178 186 195 204 213 222 231 240 249 257 266
80 173 182 191 200 209 218 228 237 246 255 264 273

There is some variation with the government site – perhaps related to rounding.

Max Heart Rate Calculations Compared

2010-10-04T16:10:00.000-07:00

Physical fitness has become increasingly technical and data driven. I started running a bit in the last few months and have been delving into the prevailing wisdom related to assessing ones health as a baseline for pursuing various fitness goals. Some of the terms related to tracking a heart rate gave me visions of white lab coats, cardiac monitors and sophisticated formulas based upon years of scientific analysis. And while there may be truth to this, the practical reality is quite a bit simpler.

In many workout routines, a target heart rate is calculated which is supposed to identify a range (usually in beats per minute) during an exercise routine that will provide optimal cardiovascular value. The basic idea is that you want a work out that is rigorous enough to derive a benefit from the exercise without harming your body. It is appealing in that it provides an objective measure to evaluate your workout. And once your workout can be measured, it is possible to set goals and work to improve your heath.

What is implied in the idea of a target heart rate is that there is some upper limit that cannot safely be exceeded. You might think that you need to be hooked up to a bunch of cardiac sensors to find out this value - and although this might be optimal, it is not the technique used by most folks. Instead, there are relatively simple formulas that are used to calculate a maximum heart rate for an individual. They are usually based only upon age (although some calculations consider gender as well).

Maximum Heart Rate Calculations
Various formulas (most of them simple linear formulas at that) have been devised to estimate individual Maximum Heart Rates. However actual maximum heart rates vary significantly between individuals based upon physiology, physical fitness and other factors so the value of the metric is disputed. Nevertheless, I was interested in comparing the available formulas to get a sense of a range (based upon "ensembling" if you will) of what is being reported or suggested by health sites, software and machines that use this value.

One of the gizmos I have begun using is the Garmin GPS with heart rate monitor. I am impressed with its performance so far. It includes its own software that does most of the types of data aggregation and summary that you would like - but I look forward to geeking out and seeing what can be done with the data in R in later posts.

Method
A ruby script was used to create a semicolon delimited file with the maximum heart rate from ages 18 through 90 for various calculation methods described in the Wikipedia article. The resulting data can be read into an R script to produce the charts in this blog.

A summary that combines the calculations combined does not make a whole lot of sense since two of the calculations in use are for women only and one is for men only. However, all of the techniques fit within a relatively narrow range (since we human beings aren't quite that random). Besides, the two calculations for women are among the most divergent presented, and so cancel each other out in part (though they probably pull down the average for they younger and older ends of the spectrum).

This average is included in the chart below - which is easier to see if you generate it yourself and stretch it to a size suitable for your monitor.

The only input value considered in calculations is gender - two of which are specific to women and one for men. It seems that the most popular calculations don't bother with gender anyway.

There are a number of possibilities for using R with fitness devices that provide heart rate information, geographic data, time, distance, caloric intake and consumption, etc. I was was not able to find much in the way of open source fitness related calculation software APIs - so this could be an new area for R developers to address. (It also provides some balance to the relatively sedentary life of developing and blogging).

Find Duplicate Records in a File

2010-09-24T16:10:00.000-07:00

In the world of data preparation a common task is to identify duplicate records in a file or data set. A few years ago, I did most development work in Java, and shudder to think of the amount of code required to accomplish this sort of task. Since that time I been involved in many projects that did not require programming in a specific language, but simply "getting the job done." With that in mind, "removal of duplicate records in a file" can be construed as manipulating a data set rather than an exercise in file processing. The following shows how R compares with other technologies when performing this task.

A file named "file.txt" containing semicolon delimited records appears in the examples below:

1;abc;123;etc

1234;qwer;4321;etc

1235;asdf;4341;etc

1;abc;123;etc

1234;qwer;5555;etc

Identifying rows in which every field is duplicated is relatively straight forward using Unix utilities. The file can be outputted using the cat (concatenate) command, the results sorted, and a unique list of results (prefaced by the count of occurrences in the file) can be filtered by a regular expression that indicates any row that has a number of occurrences not equal to one.

cat file.txt | sort | uniq -c | grep "^ [^1 ]"

Often, it is more concise to speak Unix than English. This starts to break down a bit when considering each line in the file as a record with distinguishable fields.

General purpose programming languages can do the same thing, but are a bit more verbose. In ruby, the file can be read into an array of arrays.

a=File.open('file.txt').readlines.map{|l|l.split(';')}

A list of unique rows can be outputted using the following one liner:

a.uniq.each{|l|puts l.join(';')}

And with a bit more effort, you can write a program that will filter the results as needed. This type of processing can also be done declaratively in R where such a file is read in as a data frame.

df=read.csv('file.txt',sep=';',header=FALSE)

unique(df)

The duplicated function can also be used to identify the single row that is duplicated.

df[duplicated(df),]

That is to say, the all fields in all records in the following data frames are equal.

unique(df)==df[!duplicated(df),]

The situation gets a bit more complicated when you want to only use some of the delimited fields to identify duplicate records. In the data set above, rows 1 and 4 are identical. Consider the requirement to recognize lines 2 and 5 as identical (due to the first two fields matching). In Unix, you could use awk with the -F option to process the delimited fields. In ruby you could store key fields in variables and compare them with each row. If you come from the SQL world, you could use the R sqldf package to treat the data frame as a database table.

sqldf('select * from df group by V1, V2 having count(*) >1')

The R way of getting this information is to identify the indices of the duplicated rows.

df[duplicated(df[c(1,2)]),]

I enjoy looking at the overlapping aspects of programming languages and utilities (like this OTN Article from a few years ago). It is helpful for highlighting the right tool for the right job, and aids in communication with other technical professionals.

New World Bank Data Available

2010-09-22T16:19:00.000-07:00

Just announced: World Bank Data features and data are available. Previous posts have demonstrated how to access and plot this data using R (including the use of the R WDI package). The chart above can be created using the following program in R.

library(ggplot2)
library(WDI)

DF <- WDI(country=c("US","FR","DE","GB","CN","RU","IN"), indicator="SP.POP.TOTL", start=1990, end=2008)

ggplot(DF, aes(year, SP.POP.TOTL, ggplot(DF, aes(year, SP.POP.TOTL, color=country))+geom_line(stat="identity")+theme_bw()+xlab("Year")+opts(title="Total Population")+ylab("")

Elder Research Two Day Course

2010-09-18T11:29:00.000-07:00

... or what I did on my summer vacation...

Just got back from the Elder Research Two Day Course "Tools for Discovering Patterns in Data". It was a great course that (while not R specific) provides a great overview of Data Mining tools and techniques and insight into current applications in a wide variety of industries.

Dr. Elder is a coauthor of a book available online (and provided with the course) called "Handbook of Statistical Data Analysis and Data Mining Applications." This book contains a wealth of practical examples and tutorials (most using the Statsoft Statistica software). It has a decidedly practical emphasis that allows you to see how algorithms are used to discern patterns in the data and to evaluate and compare how effective they are with specific data sets. Functional areas covered in the tutorials include aviation safety, movie box office receipts, customer services, credit scoring, automobile brand review, quality control, business administration in a medical industry, psychological evaluation, dentistry and profit analysis. This is very helpful for those who prefer to work from the concrete to the general (rather than being provided mathematical abstractions that you then apply to specific situations). They might also be helpful for showing a business user why data mining matters and what value it brings to a business or organization.

The conference covered many of the same topics discussed in Introduction to Data Mining by Tan, Steinbach and Kumar. However, there were many more concrete examples and applications of techniques in specific areas of finance, industry, government and education. A section of the book on ensemble methods is included in a larger section simply titled "Classification: Alternative Techniques". Dr. Elder went into greater detail on these topics and demonstrated the effectiveness of combining multiple models into a single model that is usually more accurate than the best of the individual component classifiers. It seems that different classifiers "see" certain parts of data sets better than others, and that combining classifiers results in a final analysis in which the best (most accurate) elements of each classification are retained while the worst aspects are largely ignored. By combining classifiers and manipulating the training set and input features a more accurate final model can be obtained.

More detail about ensemble methods is available in another book coauthored by Dr. Elder entitled Ensemble Methods in Data Mining. This book goes into greater detail about how and when to use ensembling and includes some examples in R. The use of multiple classification techniques raises a number of interesting issues - on the one hand they seem to work in practice, but there use makes it more difficult to trace how a final combined model is constructed from the original data set. This has raised some interesting issues about the definition of complexity and the quest for simple accurate models.

Dr. Andrew Fast presented on Text Mining and Social Network Analysis - and provided some valuable insights into these rapidly developing fields. There were also a number of software demos and time to interact with other members of Elder Research staff and conference participants.

The conference took place in Charlottesville VA which is a great setting with many historical and recreational attractions nearby.

So that's what I did on my summer vacation...

Ah Bach...

2010-09-07T17:54:00.000-07:00

As announced by David Smith over at Revolution Analytics, a ggplot2 Case Study Competition is on...

Rather than blogging for the last few days, I cobbled together an entry. It is not a particularly mind bending use of ggplot2, but the subject matter is relatively original. It is an brief analysis and visualization of a J.S. Bach 2 Part Invention. And because Bach's music is so well structured, the visualization itself is nice looking and well balanced. Perhaps suitable for geeky tee shirts...

Check it out when you get a chance.

Bot Botany - K-Means and ggplot2

2010-09-02T17:34:00.000-07:00

So if you had a robot that was an expert at botany - would you have a bot botanist? Among other things, it would need to to distinguish flowers through vision and image processing, and be able to classify various kinds of plants based upon specific characteristics. What do both of these requirements have in common? They can be done using the k-means clustering. Image segmentation can be used to allow our robot to recognize objects. Based upon petal and sepal size, it could determine say - the species of an iris. The well-known iris data set has been featured in other posts.

K-Means in R
If you look up the k-means algorithm online or in a reference book, you will be met with a flurry a mathematical symbols and formal explanations. The basic principal (informally stated) is rather simple... given set of observations (picture a scatter plot of points), and a number of groups or clusters that you wish to group them in, the k-means algorithm finds the center of each group and associates observations with the groups with the "closest" center.

To use k-means in R, call the kmeans function with a matrix of values and the number of centers. The function seeks to partition the points into k groups (the number of centers) such that the sum of squares from points to the assigned cluster centers is minimized. Each observation (point) belongs to the cluster with the nearest mean.

How-To
To start, we will copy the iris data set to a separate data frame. Not strictly speaking necessary, but makes it easier me to reflexively enter df whenever the data frame is in view. Next we create a matrix object containing only the Petal Length and Width.

df=iris
m=as.matrix(cbind(df$Petal.Length, df$Petal.Width),ncol=2)

Now we will do the actual clustering.

cl=(kmeans(m,3))

Simple eh? The cl object contains a number of interesting attributes associated with the model.

cl$size
cl$withinss

Next we do a bit of data formatting and preparation for subsequent calls to graph the data. Notice that we add the cluster information back to our original data frame. This is a good organization of the data and also a requirement for working with ggplot2 which is designed to use data frames.

df$cluster=factor(cl$cluster)
centers=as.data.frame(cl$centers)

The following graph color codes the points by cluster. We also add the centers and a semi transparent halo around the center to emphasize the place of the center... and its role in classifying the observations into clusters.

library(ggplot2)

ggplot(data=df, aes(x=Petal.Length, y=Petal.Width, color=cluster )) +
geom_point() +
geom_point(data=centers, aes(x=V1,y=V2, color='Center')) +
geom_point(data=centers, aes(x=V1,y=V2, color='Center'), size=52, alpha=.3, legend=FALSE)

This plot is an interesting example of how several different sets of data (in this case the actual observations as well as the centers) in separate data frames can be included in a single ggplot2 chart.

Misclassified Observations
The models is pretty accurate, but not perfect. The following SQL statement highlights the few misclassified observations:

sqldf('select Species, cluster, count(*) from df group by Species, Cluster')

Species cluster count(*)
1 setosa 2 50
2 versicolor 1 48
3 versicolor 3 2
4 virginica 1 6
5 virginica 3 44

So we grab the outliers into their own data frame....

df2 = sqldf('select * from df where (Species || cluster) in

(select Species || cluster from df group by Species, Cluster having count(*) < 10)')

Now we can enhance the previous graph to put a diamond around misclassified points.

last_plot() + geom_point(data=df2, aes(x=Petal_Length, y=Petal_Width, shape=5, alpha=.7, size=4.5), legend=FALSE)

And so with a bit of Data Mining knowledge and the R programming language, even our machines can stop and smell the roses...

Better than Average

2010-08-31T17:51:00.000-07:00

The NIST's The Engineering Statistics Handbook includes an Introduction to Time Series Analysis which provides a great way of demonstrating how R can be used to make such calculations. This post replicates the analysis of the data set introduced under Averaging Methods using R.

As you might expect, Time Series Analysis is a broad subject that has been investigated in depth elsewhere. If you need more information, a book such as Time Series Analysis and Its Applications provides a much more in depth look at the mathematical theory involved as well as providing practical examples of the use of R for analysis and forecasting.

But back to the NIST handbook... the data set they used represents supplier deliveries to a warehouse. The calculations that follow demonstrate how to perform the calculations they do in this section of the handbook using R.

Supplier Amount(in 1000 of $)

1 9

2 8

3 9

4 12

5 9

6 12

7 11

8 7

9 13

10 9

11 11

12 10

Simple Average (Mean)
In R the series can be represented as a vector.

v=c(9,8,9,12,9,12,11,7,13,9,11,10)

The average of the series is 10.
mean(v)

The "error" amount that each entry in the vector differs from the mean can be calcuated as follows.
s - mean(s)

This value can serve as the basis for a measure to ascertain how well a model fits (Error Squared).
(v - mean(v))^2

Finally, the sum or mean of these results can be used to compute values that represent the overall fit (or amount of error) for the estimate.

sum((v - mean(v))^2) # SSE" is the sum of the squared errors.
mean((v - mean(v))^2) # MSE" is the mean of the squared errors.

Now that we have a simple values that indicate how good an estimate for a set is, we can test with other values.
Rather than writing out an entire calculation each time, we can create a function in R and apply the function to each value in a vector.

sse = function(x, series){sum((series - x)^2)}
mse = function(x, series){mean((series - x)^2)}

To compare the estimate (10) with 7, 9, and 12.

sapply(c(7,9,10,12),sse,v)
sapply(c(7,9,10,12),mse,v)

Analyzing Time Series Data
A time series is simply a sequence of data points in time. Time series data has unique characteristics which allow it to be processed in a similar manner regardless of the underlying data represented. Many disciplines deal with this type of data including statistics, signal processing, econometrics and mathematical finance. Such data appears in business in relationship to sales forecasting, budgetary analysis, yield projections, and in the process / quality control arena. In other blog entries, they are used in relation to stock market analysis and economic data. They are relevant to web sites and are available through tools like Google Analytics.

So time series data is widely applicable but has common features regardless of its application. It can be analyzed to identify its characteristics and patterns. This often leads to forecasting in which a model is used to
predict future events based upon past data.

All time series data has the following common qualities:

a natural temporal ordering
often events that are close together are generally more closely related than those further apart
in most cases, past values are assumed to influence future values (rather than the other way around)
usually spaced at uniform intervals

The data set we are working with is a bit odd to consider as a time series - a supplier is not a unit of time. However, it is useful for making the point that a "simple" average (or mean) of all past observations is only a useful estimate for when there are no trends. Not sure what to make of this. I emailed the government and asked for clarification. Will post the answer here if I receive a response.

In R, a vector can be cast to a time series object as follows:

s=as.ts(c(9,8,9,12,9,12,11,7,13,9,11,10))

Moving Average
A moving average is described in the NIST Handbook and is also referred to as "smoothing" - a term that comes up in ggplot2 (geom_smooth). There are a myriad of functions available in R that involves some sort of lagged calculation of a series of numbers. A simple example that almost does the trick involves rollapply:

rollapply(s, 3, mean)

This works, but it is not clear that the first two entries were skipped. Better to use a library that has additional checks coded in...

library(TTR)
SMA(s,3)

If you take a look at the code inside... you can get an idea of the additional verification and error checking (which accounts for missing values at the beginning of the list). To view the source, simply input the function name without any parenthesis:

SMA

You can drill down into the internally called methods in this case:

runMean
runSum

With this method available, we can calculate the Error and the Error Squared:

s - SMA(s,3) # Error
(s - SMA(s,3))^2 # Error Squared

Note that the calculated mean replaced missing entries as zeroes...
x=((s - SMA(s,3))^2)
x[ is.na(x) ] <- 0
mean(x)

Oh - in case you were interested in the plot:

library(ggplot2)

df = as.data.frame(as.ts(v))

df$idx = as.numeric(rownames(df))

df$x= as.numeric(df$x)

qplot(data=df, idx, x) + geom_line() + geom_smooth()

Fractals in R

2010-08-27T15:21:00.000-07:00

Atte Tenkanen had a blog on fractals using R for a time. Much of his source code is still available online. To produce his version of the Mandelbrot set:

source('http://users.utu.fi/attenka/mandelbrot_set.R')

Fractals (such as the Mandelbrot Set pictured above) are objects that display self-similarity on all scales. Fractal are mathematical concepts with practical applications. For example, fractal dimensions provide a solution to measuring a coastline - where you come up with different lengths depending upon the length of the ruler you use. (the Coastline Paradox).

There also pretty and fun to look at....

Back in the 80's my little brother and I would type in fractal equations into a Radio Shack Color Computer. We would spend half a day typing in a program and debugging it and the other half a day waiting for the image to render. He later improved the process by using a Commodore Amiga and more than two fingers for typing.

You can also do this type of plot using ggplot2 - and most of the effort is related to turning off axes and legends. The source is on github and can be run from there.

source('http://github.com/ezgraphs/R-Programs/raw/master/mandelbrot.R')

So now I can download source code from somewhere out on the interweb thingy and use open source software to render these images in a matter of seconds.

Just like we used to do when we were kids....

How Safe is Your Money?

2010-08-24T20:07:00.000-07:00

The FDIC regularly publishes a Failed Bank List and related statistics. This post uses data in the original XLS from the FDIC web site which is formatted for human consumption to produce the charts below using R. Note that 2010 data below is incomplete.

The chart above is sort of a theme for the analysis. It is interesting to note that there were a large number of "Problem" Banks in 1990-1992 - percentage wise comparable to the last few years. However the total number of banks has dropped every year since 1990 (indicated by the color above). This is true of both Commercial Banks and Savings Institutions (Value indicates the number of institutions).

Creation of New Charters has been decreasing and has not opposed this trend.

Consolidation been one reason in the decline in the total number of institutions (especially during the 1990s), but these have decreased during the last several years.

The number of Problem Institutions in the early 1990s (and their assets) was even higher than the last few years.

Dollar amounts are expressed in billions (on the following and subsequent charts).

Combined Deposit Insurance Fund
The Combined Deposit Insurance Fund is used to guarantee the safety on deposits in member banks. The amount of Insured Deposits has been growing.

However, the Fund Balance itself (along with the Reserve Ratio) been negative since 2009. The balance was also negative in 1991.

The Number of Failed Institutions (described in detail on the Failed Bank List) is actually less than the early 1990s, but the total assets involved are much higher.

The amount of assistance provided has also been greater during the last few years. There actually was assistance provided in 1990-1992, but the relative amount is so little that it does not appear in the charts below.

Definitions do not appear in the spreadsheet, but on other parts of the FDIC site, the "estimated loss is the difference between the amount disbursed from the Deposit Insurance Fund (DIF) to cover obligations to insured depositors and the amount estimated to be ultimately recovered from the liquidation of the receivership estate". Simply said, the amount paid out by the government that could not be recovered in some way from the failed bank.

The "Resolution Receivables" are also not defined in the spreadsheet. But elsewhere on the FDIC site appear to represent "payments made by the DIF to cover obligations to insured depositors, advances to receiverships and conservatorships for working capital, and administrative expenses paid on behalf of receiverships and conservatorships".

Number of Employees
The number of FDIC employees has been adjusted in some way - and so years prior to 2008 have been revised (according to the note in the spreadsheet) "according to a new full-time equivalent methodology".

Technical Notes
The original XLS from the FDIC web site was saved as a csv using Excel. This csv and the R Source used to create the plots above are available on github.

Data Preparation
Many texts on data analysis or data mining make some comment about the amount of time used to prepare data. It needs to be "cleansed" or "reshaped" or otherwise transformed into a form that can be manipulated in the intended manner. Because the spreadsheet from the FDIC is directed at a human audience (rather than a computer program) it has a number of irregularities. There is text information in the header and footer, subheadings over related sections, columns added for space, and spaces added to cells. In addition, naming the rows and columns and reshaping the data make it easier to plot. There are a number of mechanisms available using native R and available packages. I ended up using several of them - but because of overlapping functionality, there are many ways this work could be accomplished. A few examples are below, the rest can be seen in the source file.

To remove empty columns reference data frame column indices as a list of negative numbers.
df=df[c(-2,-5,-7,-10,-12,-14,-16,-18,-20,-22,-24,-26,-28,-30,-32,-34,-36,-38,-40,-42)]

Rows can be removed in a similar fashion.
df=df[c(-1,-5,-9,-10,-13,-14,-18,-21,-24, -27),]

To remove trailing rows reference a negative sequence as the first parameter in the data frame.
df=df[-(29:36),]

To change the column names, substitute out characters
colnames(df)=gsub('X','',colnames(df))

or explicitly set column names.
colnames(df)[1]='Statistic'

I am kind of fuzzy about the specific rules - but columns must be cast as factors in certain contexts and as characters in others to get the desired result. If there is a set rule about usage, please let me know. At this point, I am assuming that these types are not handled in a consistent way and you just have to know how a particular function behaves.

df$Statistic=as.character(df$Statistic)

To trim out leading space in all cells of a given column use an appropriate regular expression.

df$Statistic=gsub('^ *','',df$Statistic)

To include subheading information include it inline. This can be done by explicit naming of rows
df$Statistic[11]='Problem Institutions Number'

or in a somewhat generic fashion where the contents of one field are appended to another.
df$Statistic[3]=paste(df$Statistic[2],df$Statistic[3])

To pivot data, you can sometimes do a transpose and the result is a matrix.

t(df)

A better result is obtained using the reshape package.

library(reshape)
df.melted=melt(df, id="Statistic")

The data is to be treated as numeric. To strip out all commas and cast the values to numerics

df.melted$value=as.numeric(gsub(',','',as.character(df.melted$value)))

Again, casting factors as numerics produces funky results in any case.
df.melted$Year=as.numeric(as.character(df.melted$Year))

Subsets of each series contained in the melted data can be obtained by matching on the Statistic in view.

fdic_employees = df.melted[df$Statistic==
'Number of FDIC Employees***',c('Year','value')]

The chart at the top was obtained by combining results from several of these subsets. The sqldf package is used to combine the data frames and do some calculations

pin=df.melted[df.melted$Statistic=="Problem Institutions Number",]

cb= df.melted[df.melted$Statistic=="Commercial Banks",]

si=df.melted[df.melted$Statistic=="Savings Institutions",]

library('sqldf')

pct_prob=sqldf('SELECT

p.Year Year,

p.value prob, c.value + s.value "Total Institutions",

p.value / (c.value + s.value) * 100 "Percent Problem"

FROM pin p

JOIN cb c ON c.Year = p.Year

JOIN si s ON s.Year = p.Year')

Plotting the Data

The plotting could be done in individual chunks of code like the following.

library(ggplot2)

p = qplot(data=pct_prob, Year, `Percent Problem`, color=`Total Institutions`, group=1)

p = p + geom_line()

p + opts(axis.text.x=theme_text(angle=-90, hjust=0), title='Total FDIC Problem Institutions by Year')

ggsave('totalFDICProblemInstitutionsByYear.png')

However, many of the charts had the exact same format and simply applied to a different series. So I created a function that would produce the desired chart for any supplies Statistic.

fdicPlot=function(statistic){

f=df.melted[df.melted$Statistic==statistic,]

ggplot(data=f, aes(x=Year, y=value)) + geom_line() + geom_point() + opts(title=statistic)

img=paste(gsub(" |\$|\$|\\*",'',i),'.png',sep='')

ggsave(img)

}

The function could then be called for each available Statistic.

for (i in df$Statistic){fdicPlot(i)}

dev.off()

All in all, a useful exercise for reviewing FDIC data a bit more closely and seeing how R can be used to munge data into a format conducive to producing charts with ggplot2.

Map of Upcoming Ruby Conferences

2010-08-21T12:09:00.000-07:00

One of the top searches on rubyflow is “conference”. A recent post showed how to create a map with the location of the 2010 R User Conference. So why not expand on the subject and create a map with numerous conference locations throughout the world?

This post shows how to create a map of locations (upcoming Ruby conferences) them straight off the web using Ruby and R.

The Packages and APIs
Both R and Ruby have a ton of functionality baked in to accomplish this task. Ruby can scrape the web using Hpricot and geocode the information using Google. It can call R through Rserve, and the maps and ggplot2 libraries can be used to render the result.

The Process
Ruby code described here is available on Github. Rserve must be running when the program is run. The following block of code is used to create a connection – and if the connection is not available, starts Rserve.

begin
  puts "Creating a new Rserve Connection."
  $c = Connection.new
rescue
  puts "Could not create an Rserve Connection: #{$!}"
  puts "Trying to start one now..."
  File.open('tmp.R','w'){|f|f.puts "library(Rserve)\nRserve()"}
  system('"R.exe" --no-save < tmp.R')
  sleep 3
  $c = Connection.new
  puts "Rserve Started."
end

Getting a list of conferences involves parsing HTML. In R the XML packages includes useful functionality and can be handy when data is in an HTML table. In this case, the data was not in an HTML table. Instead, the Hpricot parser accepts an XPath expression and iterates over the relevant elements. The data is extracted and stored in an array of hashes.

def get_conference_list()
  u='http://blog.sphereinc.com/2010/08/13-upcoming-ruby-and-rails-conferences-you-dont-want-to-miss'
  doc=Hpricot(open(u))
  recs=[]
  (doc/"//div[@id='post-216']/div/p/strong").entries.each_with_index{|e,i|
   h={}
   e.inner_text.split("\n").each{|d|
   p=d.split(':')
   unless [nil,''].include?(p[0]) or [nil,''].include?(p[1])
   # puts ">>#{p[0].strip} = #{p[1].strip}<<"
   h[p[0].strip]= p[1].strip
   end
   }
   recs << h
  }
  recs
end

Hpricot can parse XML as well as HTML, and so is used to get the latitude and longitude for each location.

  def get_location(str)
   u=URI.encode(
   "http://maps.google.com/maps/api/geocode/xml?sensor=false&address=#{str}"
   )
   loc=(Hpricot.XML(open(u)))/'//location'
   h={}
   h['lat']=(loc/:lat).inner_text
   h['lng']=(loc/:lng).inner_text
   h
  end

A data file is created which contains semicolon delimited records. This will provide the input to the R program.

Finally, an R program is used to plot the data. Since the files in view will be in the current working directory (and Rserve has no reference to this) it is substituted in prior to executing the program.

Challenges
In many cases ggplot2 creates publication quality graphics with a simple call. In this case, data was grouped very unevenly, and attempts to automatically add text (in this case the city names) to the map can result in overlapping and large amounts of wasted space.

I thought of a variety of ways to address the problem.

The simplest was to resize the image various ways. This was not successful - but I did later crop the images in a manual post processing step.
Another approach was to modify the scale in use in some way (e.g. the use of a log scale with scatter plots as Tal pointed out in this post). The ggplot2 package includes a large number of map projections (ways of representing a three dimensional sphere in two dimensions). Available projections provided with ggplot2 are described in help:

mapproject

None of these solved the problem in and of themselves, but provided some interesting variations.

I decided instead to simply “zoom in” on the relevant part of the chart and split the chart (sort of ad hoc faceting). In some cases this worked perfectly.

In other cases the results were suboptimal. Paths were sometimes drawn between points on the border that had been cropped and so extraneous map lines appear in some zooms. I finally decided that a general purpose automatic solution was not readily available (or at least known to me). So I cleaned up the final images using image editing software (Gimp).

If anybody has ideas about how this kind of processing can be used to automatically generate graphs that are "cleaned up" please comment below. Scriptable solutions that create finished products with no manual intervention are preferable.

The code is available on Github.

Programming Language Popularity: StackOverflow and Ohloh

2010-08-17T14:35:00.000-07:00

In the following example, programming language popularity is measured based upon two data sets. The first is the number of contributors associated with a language on ohloh.net. The second is tag usage at stackoverflow.com.

SQL with no DDL
I admit it... in an age of NoSQL... I like SQL. I agree that fixed table schemas can be a real pain though... who wants the overhead of defining database tables for a quick comparison of two data sets?

Joining on the language name provides a simple, intuitive way to correlate the results sets. Of course there are limitations to this approach - after all SQL was designed for relational databases. Since a "join" is being done based upon language and tag name, some languages may be under represented. For instance

JQuery - a Javascript library - is a leading tag.
Objective C questions might appear under the iPhone tag.
C# Questions might appear under .NET, ASP.NET or other Microsoft tags.

Another problem that can occur in this type of comparison is that keyed data might not actually correlate with a matched key. For instance the C programming language might be compared with a record pertaining to the third letter of the alphabet. This is not a problem in the current example because the specific domain of both data sets is programming languages.

Top Languages
In this particular analysis, I am really interested in outliers - not the vast majority of the languages that appear in the data set. So the name of each point will be plotted beside it. For less popular languages, this chart is impossible to read and madly cluttered... but it is great for focusing on the most popular languages.   So rather than coming up with a publication-quality graphic, the chart above provides a "quick-and-dirty" perspective that can lead to helpful discussions for people familiar with the programming language domain.

In the previous post, Ruby ranked at the top. This demonstrates the Ruby centric nature of github, which was initially directed towards the ruby community. Similar trends affect the results in the current post (where Ruby ranks 12th in tag count and 16th in the number of contributors). R is 18th in tag count and 33rd in number of contributors.

The data was extracted over the last few days and is available on github in ohlo_2010-08-16.txt and stackoverflow.txt (warning 400MB file... all tags from stackoverflow are listed in it). The process to analyze the files involved the following R Code.

library(ggplot2)
library(sqldf)

SODF=read.csv('stackoverflow.txt',header=TRUE, sep=';')
OHLODF=read.csv('ohlo_2010-08-16.txt',header=TRUE, sep=';')

head(OHLODF)
head(SODF)

df = sqldf('select Name name, Count tag_count, Contributors contributors
from OHLODF o
join SODF s on LOWER(s.Tag) = LOWER(o.Name) order by 1')

ggplot(data=df,
   aes(x=tag_count, y=contributors, color=name)) +
  geom_point() +
  geom_text(aes(label = name))

The resulting chart is displayed above. To list the top 10 languages:

> head(df[order(df$contributors, decreasing=TRUE),],10)
   name tag_count contributors
57 XML 12374 133183
24 HTML 21936 106012
28 Java 62386 78098
9 C 17256 78023
13 CSS 16429 72060
11 C++ 38691 61831
29 JavaScript 46608 60677
33 Make 537 50328
44 Python 31852 38691
39 PHP 53884 36952

> head(df[order(df$tag_count, decreasing=TRUE),],10)
   name tag_count contributors
10 C# 101811 22198
28 Java 62386 78098
39 PHP 53884 36952
29 JavaScript 46608 60677
11 C++ 38691 61831
44 Python 31852 38691
48 SQL 25316 28069
24 HTML 21936 106012
9 C 17256 78023
37 Objective-C 17250 6555

All other things being equal, one might think that the relationship between contributors to projects and tag counts might be roughly linear. As it stands, that is not the case at all.

Web Oriented Languages
The languages represented show a significant representation of web applications related technologies. HTML, CSS, Java Script and PHP are used almost exclusively for such development, and Ruby, Python, Perl, Java, C#, SQL are also heavily used for web applications (though not exclusively). C, C++, Objective C and Make are related technologies that are geared less towards web development.

Microsoft
According to wikipedia StackOverflow is a Microsoft partner and stackoverflow itself was developed on the Microsoft platform. This might provide some explanation to the high representation of C#.

Simple Languages = Less Questions
XML and HTML are markup languages with relatively simple syntax, hence the relatively small tag count. CSS and Make are also relatively small languages with specific uses rather than general purpose programming languages. The fact that C++ was developed as an enhancement to the C programming language explains why there are more questions (and a larger tag count) for C++ than C. A more speculative suggestion is that Perl's relatively low tag count indicates that the "more than one way to do it" philosophy leads to less questions. An obvious alternative is that Perl users simply ask questions in other venues.

Conclusion
All measures of programming language popularity have their limitations. Correlating various sets of data can provide some additional insights into their prevalence and usage. R and sqld provide a convenient means of making such comparisons. And ggplot2 provides a great way of charting results.

Update

A log scale (as suggested by Tal in the comments) provides better insight into the majority of languages that appear clustered in the lower left hand corner of the chart. However, though this site might be considered R rated, the **** was added through later image editing to make it fit for all audiences.

ggplot(data=df,

aes(x=log(tag_count), y=log(contributors), color=name)) +

geom_point() +

geom_text(aes(label = name))

GitHub Stats on Programming Languages

2010-08-09T13:36:00.000-07:00

GitHub has become a popular site for Open Source Developers to stash code and collaborate on projects. The following are some stats and analysis related to programming languages in use based upon the number of users and repositories. The data was obtained from GitHub's searches. It and the R code are available in GitHub as well (a lovely recursive relationship I must say).

> df.top_ten_reps

Language Repositories Users
1 Ruby 104239 23123
2 JavaScript 44482 10895
3 Perl 34232 2178
4 Python 32150 8775
5 PHP 21685 8872
6 Java 17687 6618
7 C 16137 5558
8 C++ 12521 5595
9 Objective-C 8027 2520
10 C# 6061 2706

Ruby has a commanding lead in terms of the number of repositories with 32.17% - more than the next two (Javascript and Perl) combined. R is ranked 25th with 191 repositories or about 0.06% and only 6 projects behind the D programming language. The top 5 are scripting languages, Java ranks 6th and the C family rounds out the top ten. Relatively open languages lead the pack, followed by those with a proprietary focus (Objective-C for Apple and C# for Microsoft).

When ranked by number of users, the top two remain the same. There is a bit of shuffling with the remainder of the top ten.

> df.top_ten_users
Language Repositories Users
1 Ruby 104239 23123
2 JavaScript 44482 10895
3 PHP 21685 8872
4 Python 32150 8775
5 Java 17687 6618
6 C++ 12521 5595
7 C 16137 5558
8 C# 6061 2706
9 Objective-C 8027 2520
10 Perl 34232 2178

The most striking is that Perl drops to 9th place. There are significantly less users associated with Perl - particularly for the number of projects. I noticed that there was a migration of Perl language source code to GitHub - so perhaps modules were migrated as well...but I couldn't find any specific announcements that clarified this.

All other things being equal, you might expect there to be a relationship between of the number of users to repositories. There is to some degree -

> df.Ratio=df.Repositories / df.Users
> mean(df[df$Ratio > 0 & df$User > 0 & !is.na(df$Ratio), ])

5.064732

A linear model suggests a slightly lower value (between 4 and 5).

Here is a plot restricted to the Top 10.

Some of the lesser used languages have few users and more repositories like IO (19 per user) and CoffeeScript (17 per user). Perl has a remarkable 15 per user.

Full data Set
Language Repositories Users Rep.pct Ratio
1 Ruby 104239 23123 32.17 4.508022
2 JavaScript 44482 10895 13.73 4.082790
3 Perl 34232 2178 10.56 15.717172
4 Python 32150 8775 9.92 3.663818
5 PHP 21685 8872 6.69 2.444206
6 Java 17687 6618 5.46 2.672560
7 C 16137 5558 4.98 2.903383
8 C++ 12521 5595 3.86 2.237891
9 Objective-C 8027 2520 2.48 3.185317
10 C# 6061 2706 1.87 2.239837
11 Shell 4657 1011 1.44 4.606330
12 VimL 4248 1267 1.31 3.352802
13 ActionScript 2609 1104 0.81 2.363225
14 Erlang 2520 532 0.78 4.736842
15 Haskell 2290 641 0.71 3.572543
16 Scala 2154 539 0.66 3.996289
17 Clojure 2082 481 0.64 4.328482
18 Lua 1754 511 0.54 3.432485
19 Groovy 870 261 0.27 3.333333
20 Scheme 707 140 0.22 5.050000
21 Go 398 103 0.12 3.864078
22 OCaml 382 121 0.12 3.157025
23 Objective-J 355 109 0.11 3.256881
24 D 197 64 0.06 3.078125
25 R 191 69 0.06 2.768116
26 ColdFusion 180 56 0.06 3.214286
27 Tcl 125 39 0.04 3.205128
28 ooc 112 11 0.03 10.181818
29 FORTRAN 93 47 0.03 1.978723
30 ASP 88 35 0.03 2.514286
31 Smalltalk 80 14 0.02 5.714286
32 HaXe 75 14 0.02 5.357143
33 F# 74 5 0.02 14.800000
34 Verilog 74 26 0.02 2.846154
35 VHDL 64 14 0.02 4.571429
36 Io 57 3 0.02 19.000000
37 SuperCollider 53 11 0.02 4.818182
38 Arc 48 15 0.01 3.200000
39 Delphi 43 16 0.01 2.687500
40 Assembly 41 5 0.01 8.200000
41 Boo 41 6 0.01 6.833333
42 Nu 40 4 0.01 10.000000
43 Eiffel 39 15 0.01 2.600000
44 CoffeeScript 34 2 0.01 17.000000
45 Vala 27 3 0.01 9.000000
46 Racket 20 8 0.01 2.500000
47 Self 7 3 0.00 2.333333
48 Duby 4 0 0.00 Inf
49 Max/MSP 4 2 0.00 2.000000
50 sclang 2 0 0.00 Inf
51 Common Lisp 0 0 0.00 NaN
52 Emacs Lisp 0 0 0.00 NaN
53 Pure Data 0 0 0.00 NaN
54 Visual Basic 0 0 0.00 NaN

Iris Data Set Visualization Web App in < 100 LOC

2010-08-07T14:11:00.000-07:00

The iris data set pops up pretty regularly in statistical literature. It consists of 50 records from three species of Iris flowers (Iris setosa, Iris virginica and Iris versicolor). I came across it recently while reading Introduction to Data Mining. It comes up in several places in the book to demonstrate techniques for visualization and classification. There has been a number of articles, posts and videos about R and the web in recent times. This post presents a way of creating some plots for the data set using R and a Ruby Web Application.

R and the Web
There are a number of situations where it would make sense to expose a data set and wrap a certain amount of R functionality within a web application. Non-R users might need access to the data. You might want to provide a presentation of findings available through the web. You might even want to collaborate with other R developers by posting an HTML table that they can read in using XML. In time, I expect that some standard web application frameworks will emerge to fill in the gap. And so the current application is some “thinking out loud” on my part in that direction.

Prerequisites
In order to run this application on your machine, R and Ruby must be installed and functional. The R packages ggplot2, R2HTML, RServe are used as well as the iris data set. This web app was written and tested in Windows, but should run in *nix with small modifications.

Ruby Configuration
Install the package to allow communication with RServe

   gem install rserve-client

On *nix systems folks often sudo to install.

If for some reason you do not go the normal route of installing a gem (e.g. you downloaded from github at http://github.com/clbustos/Rserve-Ruby-client), make sure that your ruby $LOAD_PATH has the library available when you run the program.

You can even do this in line in the ruby program by including a line like the following at the beginning of the program:

$LOAD_PATH<<'C:\clbustos-Rserve-Ruby-client-v0.2.4-2-g47b0da7\lib'

The application itself is available on github.
Run the Web Application

Start Rserve.
  C:\Program Files\R\R-2.10.1\library\Rserve> Rserve

You should see output like this if it started successfully.
Rserve: Ok, ready to answer queries.

Start the Ruby web application - specify a port if you like with the -p option.
iris_data_set_webapp.rb -p 4445

With these steps complete, you should be able to hit the application at http://localhost:4445. Three links are available. The r_version link simply demonstrates that Ruby (through the sinatra framework and Rserve client) can communicate with R. Clicking this link causes the version of R to display in the browser.

The second link is to the iris data itself. This page displays a formatted HTML table rendered using the R2HTML package. Admittedly, this is a bit of a confusion of concerns (view information be provided by R) but it provided a convenient mechanism to convert a data frame to an HTML table.

The third link allows you to modify the aesthetics of the plot. Specifically, the x, y and color can be set to any of the available variables. The result is a "grammatically correct" chart.

Code Walkthrough and Commentary
The package declarations could be written out
   require 'rubygems'
   require 'sinatra'
   ...

Instead all of the packages are included in an array (surrounded by brackets). Then each require directive is issued as we iterate through each element in the array..

['rubygems', 'sinatra', 'rserve','fileutils','haml'].each{|r|require r}

The packages being used are

rubygems - the ruby packaging system itself
sinatra - a minimal web app DSL
rserve - to integrate with R
fileutils - some convenience methods for file system access
haml - well, this one requires some explanation...

HAML is one of the many Ruby mark up/templating languages that is in vogue today among Rubyists. It seems to save a few keystrokes from writing straight HTML, but it slows me down since I think in HTML and end up working backwards to writing the HAML. I kind of like the pythonesque interpretation of indentation being meaningful and that the HTML looks pretty.

Anyway, it is used here but I am still on the fence about it.

To experiment with haml using irb, just require haml, create an engine and output the results to HTML.

require 'haml'
Haml::Engine.new('%h3 hello world').to_html

Back to the web app. Create a global connection to Rserve.

include Rserve
$c = Connection.new

The following lines kinda-sorta reload Sinatra most if the time which allows you to change code and view the changes without starting and stopping the server. Only it does not always work :) …but it works enough for me that I included it and just restart if things are not updating the way I expect.

configure do
  Sinatra::Application.reset!
  use Rack::Reloader
end

This looks at a line that comes from the web app source file itself. Yep, kind of wild (echoes of camping). If the line matches the regexp and is one of the get functions below (other than the index itself), we pull out the url path and slap it in an HTML anchor. This is a convenient way to have a home screen during development where each get URL can be invoked.

def anchor(line)

  if line=~/get \'\/([a-z|A-Z])/
   l=line.split[1].gsub("'",'');
   haml "%a{:href => '#{l}'}> #{l} \n%br"
  end

end

Return a string of html with an heading that says “Links” and a link to each “get” URL available in the web application. The list of links is generated by reading the contents of this file, and creating a hyperlink (if possible) using the “anchor” method above.

get '/' do
  html=haml '%h3 Links'
  File.open(__FILE__).readlines.each{|l|html+=anchor(l).to_s}
  html
end

This is a simple example of how integration with R works. The connection to RServe named $c is sent a string of R code to evaluate. We expect a single result that we interpret as a string.

get '/r_version' do
  $c.eval("R.version.string").as_string
end

This method creates an R script, evaluates it and returns a link to an image that will appear in the public directory that is in the same directory with this file. The < < SCRIPT syntax is sometimes called a heredoc. It is just a convenient way to create multi line strings - you could use double quotes in this context as well. The variables x, y and color that are passed in are substituted where you see #{x}, #{y}, #{color}.

def irisplot(x,y,color)
   script= < < SCRIPT
library(ggplot2)
  ggplot(
   data=iris,
   aes(x=#{x}, y=#{y}, color=#{color})
   ) + geom_point()
   ggsave('#{FileUtils.pwd}/public/irisplot.png')
SCRIPT

  $c.eval(script)

  " < img src='irisplot.png' width='600', height='600' > "

end

Note: The spaces between the less than signs for the HEREDOC are artificial - they were required because blogger was not correctly interpreting them together. Similar problem with the image tag - I just added spaces to prevent rendering issues.

This example demonstrates how to open an independent R script and run it. See the iris.R script itself for more information about what is going on. In general, the R2HTML package is being used to create a file whos handle is returned. We then read the contents of the file in and these are returned as HTML. The contents of the file are an HTML table that represents the iris data frame.

get '/iris_data' do
  url=$c.eval(File.open('iris.R').readlines.join("\n")).as_string
  File.open(url).readlines()
end

This page can be scraped using R and two lines of code. You could read this data into R running on another computer on the network:

library(XML)
df=readHTMLTable('http://nameofmachine:4445/iris')

The following creates an iris plot using the parameters passed in.

post '/plot' do
  irisplot(params['x'],params['y'],params['color'])
end

Finally, returns the form that allows you to input which fields are used to create a plot for the iris data set.

get '/iris_plot_input' do

  # Retrieve the iris data set column names into a ruby class variable.
  # These will be used to populate dropdowns.
  @colnames=$c.eval('data(iris);colnames(iris);').as_strings

  # WARNING Hard Coded Defaults below. I used these so that
  # we would have reasonable values by default.

  # Put all of the HAML markup in a string
html=<
%form{ :action => "/plot", :method => "post"}
  %table
   %tr
   %td
   %label{:for => "name"} x:
   %td
   %select{:name=>'x'}
   = @colnames.each do |col|
   %option{:value=> col, :selected => (col == 'Sepal.Length')}
   =col
   %tr
   %td
   %label{:for => "name"} y:
   %td
   %select{:name=>'y'}
   = @colnames.each do |col|
   %option{:value=> col, :selected => (col == 'Sepal.Width')}
   =col
   %tr
   %td
   %label{:for => "name"} color:
   %td
   %select{:name=>'color'}
   = @colnames.each do |col|
   %option{:value=> col, :selected => (col == 'Species')}
   =col
   %tr
   %td
   %input{:type => "submit", :value => "Create Plot"}
HAML

# Render it with the HAML engine
haml html

end

Examples produced by the app:

This application is amazing in that it simply pulls together some of the best programming resources around. In well under 100 lines of code it is simple and easy to maintain. With that in mind, I have been thinking about other directions that could be used to generalize this approach. One is simply to bundle Sinatra with R (perhaps using JRuby). Sinatra web apps could then be dynamically based upon data sets (kind of like the current app) or around R functions (kind of like the fgui package). It seems like Hadley Wickham had a similar idea first and has a related project on Github. His approach is to port Sinatra to R so that web apps could be developed in R without the use of another language such as Ruby.

U.S. Unemployment Data: Animated Choropleth Maps

2010-08-04T14:59:00.000-07:00

I was curious about creating maps with ggplot2 - and was happy to find that great minds have long since tackled this topic. Folks who are interested in how to learn about R (and ggplot2) might be interested in the process. (There should be a "greatest hits" listing related to R topics like this - the best I can think of to do is revive them in a forum like this).

After looking through the ggplot2 book (which had some simple examples) and the online documentation, I did a few searches and came across the topic on stackoverflow. At the end of this post there is a reference to previous discussions in the R community. It seems that an article on Flowing Data (using Python) inspired a bit of discussion of implementations of choropleth maps with R. And so a challenge was issued to do this in R. Several solutions were created and reviewed - and I decided to replicate the first response. Fortunately Hadley Wickham included his solution in his github repository. A slightly modified version of his original chart is shown below.

My revised version included a title, legend name and call to coord_map. I followed the same process as the original challenge - downloaded the data from the Bureau of Labor Statistics and saved the relevant sheet as a csv.

m=ggplot(choropleth, aes(long, lat, group = group)) +
   geom_polygon(aes(fill = rate_d), colour = alpha("white", 1/2), size = 0.2) +
   geom_polygon(data = state_df, colour = "white", fill = NA) +
   scale_fill_brewer(pal = "PuRd", name="Rate")

m + coord_map() + opts(title="US Unemployment Data (2009)")

Worked like a charm. Note that there are some known problems with this solution. Some additional mapping/data cleanup is required to get a complete view of the data set (which explains why Louisiana is missing). I thought it would be neat to animate the results. So I subcontracted out animation using Gimp to a local expert. He created the following - a quick attempt to illustrate the changes over time. It is not ideal, you kind of have to pick a point stare at it and watch for changes over time.

Incidentally, blogger does not support animated gifs directly, ~~so I added the image in Google docs and linked to it from there~~. This does not seem to work due to authentication issues. Uploaded it to Picasa and you can view it there instead.

I know that there is an R animation package, but was not sure of its applicability in this scenario. How would you create an animation that illustrates changes on a map over time?

Suggestions for R-Chart iPhone App?

2010-08-03T14:52:00.000-07:00

Have any ideas for news sources to include in the new R-Chart iPhone App? I am hoping to make some revisions and fixes and figured I would get some feedback from anyone who happened download it. Specifically, I am planning on removing the message button and was thinking of including another RSS feed. What news sources on R would you like to see included? I am considering adding CRAN...although these updates already show up under the Twitter listing. Any other ideas?

RODM Article on the Oracle Technology Network

2010-08-03T14:38:00.000-07:00

A new article over at OTN demonstrates the new RODM package that allows you to control Oracle Data Mining within the database itself using R. This and the recent announcement over at Revolution Analytics both reveal more solutions for processing large data sets and suggest new applications for the use of R.

Thinking about Graphs

2010-07-30T16:10:00.000-07:00

A recent Wall Street Journal article ruminated about the degree that language shapes thought (rather than the other way around). This idea has rather profound implications in the more specific domain of programming languages. We initially learn a programming language but later “think” in terms of the language.

To some degree, we are constrained in our ability to solve problems if we only know a single language. This situation has been recognized different ways by the programming community. The Logo programming language was built based upon constructionist learning theory and was intended to provide a “mental model” for children to come to understand mathematical constructs. In recent times, many programmers have committed to being polyglots, learning new languages as a part of professional development. Their concern is not always to learn the latest language that they will need to work, but to find out new ways of conceptualizing problems and structuring solutions. This leads to a more subtle goal of ggplot2.

The ggplot2 package is appealing because it makes it possible to quickly create appealing graphs and charts. However, it is based upon an underlying “grammar of graphics”. This “Grammar” serves a number of purposes. It provides a structure for the API implementation. The API is designed so that you specify what you want rather than how to create it.

Another, perhaps more subtle effect is that it also can influence the way that an R programmer thinks about creating a graph. With this in mind, it is helpful to “think through” the process of creating a chart in the terms presented by ggplot2 in a more disciplined fashion.

Components of a Plot

According Hadley Wickham (the author of ggplot and the ggplot book), the following components make up a plot:

Data
Aesthetic Mappings
Geometric Objects
Statistical Transformations
Position Adjustment
Faceting
Coordinate System

The Reference Manual is also organized around these components:

Geoms (Geometric Objects)
Statistics (Statistical Transformations)
Scales
Coordinate System
Faceting
Position Adjustment

He has organized the material slightly differently in a presentation at Vanderbilt.

Default Data Set
Set of Aesthetic Mappings
Multiple Layers (points, jittered points, box plots, histogram
Scale for Each Aesthetic
Faceting Specification
Coordinate System

In this case a layer comprises several of the elements listed earlier.

Data set and Aesthetic Mapping
Geometric Object.
Statistics
Position Adjustment

Data is not included as a part of ggplot2. In addition, algebra (from a component identified by Wilkinson) is not included as it in the realm of data transformation rather than actual chart creation.

The individual components of the grammar are fairly well defined regardless of where they appear on a list. The possible interactions between the components are rather complex. The construction of traditional charts are defined by a distinct combination of components. For example, the combination of geom and a stat is significant. At other times, the coordinate system is a defining factor. (A pie chart is a one column stacked bar chart that is mapped to a polar coordinate system).

Chart Geom Stat Coordinate System

Scatterplot point identity cartesian

Histogram bar bin cartesian

Pie Chart bar identity polar

Iris Data Set

The iris data set is a well known set of multivariate data introduced by Ronald Fisher in the 1930s. The first few rows of the set are as follows:

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1 5.1 3.5 1.4 0.2 setosa

2 4.9 3.0 1.4 0.2 setosa

3 4.7 3.2 1.3 0.2 setosa

4 4.6 3.1 1.5 0.2 setosa

5 5.0 3.6 1.4 0.2 setosa

6 5.4 3.9 1.7 0.4 setosa

The charts below will use the following components:

Data - iris data set
Mapping to aesthetic

x - Petal.Length
y - Petal.Width
Color - Species

Geometric Object - point

A scatterpoint that includes these components can be created using qplot. However, see Harlan's comment below (and this his blog which I appear to be echoing) - this is probably not the best way to start if one wants to "think" in the grammar rather than simply produce a good looking graph in the smallest number of keystrokes. Better to use ggplot as demonstrated later in the post.

library(ggplot2)

qplot(Petal.Length, Petal.Width, data=iris, color=Species)

The main decision that needed to be made to construct this call was how to map the aesthetics. It is important to consider whether each variable being mapped is discrete or continuous to create a meaningful (and not just grammatically correct) result.

Discrete Continuous

Color Distinct color Gradient (red to blue)

Size Distinct steps Radius based on value

Shape Distinct Shape N/A

We can even include more information by mapping another attribute (sepal area – derived from sepal length times width) to size.

qplot(Petal.Length,

Petal.Width,

data=iris,

size=Sepal.Length * Sepal.Width,

color=Species)

A great deal of information can be encoded using the various data attributes. The plot gives an some indication regarding the petal length and width (based upon the position), species (based upon color) and sepal area (based upon size). However, not every value is clearly in view. There are a few changes that might provide an indication that values might overlap. A jitter is might be used. A better alternative is to set an alpha value provides a degree of transparency.

qplot(Petal.Length,

Petal.Width,

data=iris,

size=Sepal.Length * Sepal.Width,

color=Species,

alpha=0.3)

We can be more explicit about what is going on using ggplot rather than qplot. The basic scatterplot can be created and in this case will be stored in a variable.

p = ggplot(data=iris,

aes(Petal.Length,

Petal.Width,

color=Species)

) + geom_point()

With the original plot in a variable, we can add components and immediately see their effect as it is rendered. A line might help discern a trend in the original scatterplot. When applying a stat, you need to – well - think statistically. Consider the following.

p + stat_abline()

The line created doesn’t mean much – this is because it is simply a line with a slope of 1 and intercept of zero. A more meaningful line can be created by determining the line of best fit.

coef(lm(Petal.Width ~ Petal.Length, data=iris))

# this returns

# (Intercept) Petal.Length

# -0.3630755 0.4157554

p +

stat_abline(intercept=-0.363,

slope=0.416,

color='purple')

So calling a given stat did something in this case. To get it to do something meaningful required additional work.

Distinction between Grammar Components

The distinction between a statistic and geometric object is not always clear (at least in terms of the ggplot2 API). A line with a slope and intercept might be though of as a statistic or a geometric object.

p + geom_abline(intercept=-0.363,slope=0.416, color='purple')

Likewise a position adjustment (like a jitter) can be thought of as both a geometric and positional terms.

qplot(Petal.Length,

Petal.Width,

data=iris,

position='jitter') +

geom_abline(intercept=-0.363, slope=0.416)

qplot(Petal.Length, Petal.Width, data=iris) +

geom_jitter()

I point out these idiosyncrasies because – as with many formal abstractions of real world concepts – edge cases exist. For example, in Western music theory, one dutifully learns the rules of counterpoint only to find out that they are not always observed by composers in practice and that certain constructs are not easily classified. This doesn’t eliminate the usefulness to studying music theory. It simply highlights the difficulty in neatly categorizing every aspect of a specific creation in an accurate an meaningful way. And for what its worth, I think that Hadley Wickham as done a marvelous job – and appears to have taken an approach of providing an interface to underlying functionality when it appears in more than one category.

Order of Application

Note that the order in which geoms and stats are applied matters! For instance:

p+geom_boxplot()

The boxplot obscures the original points. These can be added back on after applying the boxplot.

scap+geom_boxplot()+geom_point()

This gives a glimpse of the flexibility and sophistication of the system. The fundamental elements of chart design that comprise the grammar can be combined in new and flexible ways. Not every grammatically correct possibility is aesthetically pleasing or accurate as interpreted by human perception. But ggplot2 is worth learning not only for its own sake, but for the insights it can provide into the creative activity of constructing charts and graphs.

Free R Chart iPhone App

2010-07-30T15:05:00.000-07:00

Download the free R-Chart iPhone Application for the latest news and updates from top R web sites... (and this blog too).

App Store Link

Leave a comment if you have any feedback or suggestions!