mages' blog

Show me the data! Or how to digitize plots

2012-02-27T07:30:00.000Z

I had mentioned the Guardian's data blog and the need for more data journalism earlier here. What I really like about the Guardian's approach in particular is that they share the data of their articles and encourage readers to use it.

Of course there are perfectly valuable reasons for only displaying a chart and not making the underlying data available, e.g. to generate leads, as potential customers may get in touch with you asking for the underlying data, or technology issues that don't allow you to upload data, etc.

I personally believe that when I show a chart I should also make the underlying data available. Pretty pictures give you the attention, but the underlying data will offer you an opportunity to engage with your reader on a different level. This might be similar to open source software. In most cases users don't want to see and read the code, but having the knowledge that they could provides more credibility.

Screen shot of plot digitizer using Guy Carpenter's
global property catastrophe rate on line index

Here is another reason why I should make the data available: Because it is easy to extract the data from a chart anyhow, thanks to digitizing software like the Java application plot digitizer. While in the past I may have used graph paper and a ruler, nowadays it only takes a few minutes to extract the information.

Big data seminar in London on 1 March 2012

2012-02-20T19:00:00.000Z

Removable disk packs in 1975
Source: Wikipedia, via Deutsche Fotothek

David Chan from City University is organising an interdisciplinary symposium on tackling the ‘Big Data’ challenge on 1 March 2012.

It is an open seminar trying to bring together academics and practitioners from across industry to tackle the challenges posed by "big data" - the growing amount of information that needs to be stored, searched, analysed and visualised in the digital age.

The event will take place in the Oliver Thompson Lecture Theatre, Northampton Square, London EC1V 0HB. Booking is required if you would like to attend. For more details check out the event page.

See you there.

Reshaping the IT world

2012-02-16T06:56:00.000Z

During my university time I worked on the IT help desk for a while. One day I received a call from a professor, who said that his printer had stopped working. So I asked him, if there was a message on the display and if he could read it to me. "Oh yes", he said, "it says: 'Load A4 paper.'"

Rachel King quotes a study by Cisco on ZDnet, which believes to have found out that college students and young employees under the age of 30 would rather take a lower salary than having no social media freedom, device flexibility and work mobility.

It feels like the 1960's in a lot of offices and IT departments of today. A younger generation is demanding more freedom and fun. It just not called rock music, mini skirts or of course the anti baby pill, which the generation of my professor was fighting for. That's all established now. It is the digital equivalent of those rights and I can understand that IT departments are concerned about this.

It is fascinating to observe how consumer IT goods and software enter offices, e.g. employees wanting to switch from their corporate QWERTY-smartphones to their own touch-screen devices. Of course there have always been geeks and youngsters who tried to push the boundaries of the corporate IT policy, but today they are joined by the most senior guys: the C-suite, who want to show of their new gadgets and tablets to their peers.

Risk mitigation is the key. The growing demand for a more flexible IT infrastructure has seen new companies sprung up to solve those business problems, be it to create a secure connection from private devices into corporate networks or to support open source software like Linux and R. The Economist published an insightful special report on this subject in October 2011, so did Accenture, and the German news magazine Der Spiegel titled an article "Ferrari at home and horse-drawn carriage in the office".

There is one company that has reshaped the IT world more than others over the last few years, and that is Apple. Their shareholders have been rewarded for its innovations and bravery of entering new markets, see Apple's stock market performance below. Maybe worthwhile mentioning when you have challenging conversations with your IT folks.

Loading ...
The chart above shows Apple's weekly closing share prices and volumes over time since 1984, sourced from Yahoo! Finance and displayed as a time line using the Google Visualisation API in R with googleVis. You find a very similar example in the googleVis help page of gvisAnnotatedTimeLine.

The reshape function

2012-02-09T15:00:00.000Z

The other day I wrote about the R functions by, apply and friends, which allow me to operate on subsets of data. All those functions work nicely, if the data is given in the right format. More often than not it isn't and I have to reshape the data beforehand. Thus, time to discuss the reshape function. I will focus on the reshape function in base R, and not the package of the same name.

I use Fischer's iris data set again, as it is readily available after starting R. The iris data set has 150 observation and the first 6 rows look like this:

head(iris)
  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

I would like to create a box whisker plot, showing the measurements of the observations for each of the species, as in the chart below.

I know, that if I had all measurements in one column and the dimension in another column, I could produce a graph like this in one line with lattice.

library(lattice)
bwplot(Measurement ~ Species | Dimension, data=reshaped.iris)

Hence the reshape function is what I need. From the help file I learn that I want to transform my data from a wide format into a long format (direction="long"). In the long format I would like a varibale with the measurements (v.names="Measurement"), which I get by running through the first four columns (varying=1:4). I know which measurement I am reading by looking at the column names (times=names(iris)[1:4]), and I capture the dimension names in a new variable (timevar="Dimension"). This gives me the following statement:

reshaped.iris <- reshape(iris, varying=1:4, v.names="Measurement", 
                         timevar="Dimension", times=names(iris)[1:4], 
                         idvar="Measure ID", direction="long") 

head(reshaped.iris)
               Species    Dimension Measurement Measure ID
1.Sepal.Length  setosa Sepal.Length         5.1          1
2.Sepal.Length  setosa Sepal.Length         4.9          2
3.Sepal.Length  setosa Sepal.Length         4.7          3
4.Sepal.Length  setosa Sepal.Length         4.6          4
5.Sepal.Length  setosa Sepal.Length         5.0          5
6.Sepal.Length  setosa Sepal.Length         5.4          6

That's it, I can create the lattice box-whisker plot.

In my next example I would like the measurements of length and width in separate columns and capture the flower part in a new variable, so I can create scatterplots of length against width. Tweaking the reshape statement slightly gives me:

reshaped.iris.sp <- reshape(iris, varying=list(c(1,3),c(2,4)),
                            v.names=c("Length", "Width"), 
                            timevar="Part", times=c("Sepal", "Petal"),
                            idvar="Measure ID", direction="long")

head(reshaped.iris.sp)
        Species  Part Length Width Measure ID
1.Sepal  setosa Sepal    5.1   3.5          1
2.Sepal  setosa Sepal    4.9   3.0          2
3.Sepal  setosa Sepal    4.7   3.2          3
4.Sepal  setosa Sepal    4.6   3.1          4
5.Sepal  setosa Sepal    5.0   3.6          5
6.Sepal  setosa Sepal    5.4   3.9          6

xyplot(Length ~  Width | Species, groups=Part, 
       data=reshaped.iris.sp, auto.key=list(space="right"))

Let's swap Part against Species.

xyplot(Length ~  Width | Part, groups=Species, 
       data=reshaped.iris.sp, auto.key=list(space="right"))

I think, the charts illustrate quite nicely why the iris data set has become a typical test case for many classification techniques in machine learning.

googleVis 0.2.14 is released

2012-02-05T20:31:00.000Z

Version 0.2.14 of the googleVis package was released on CRAN today.

Changes

The help files have been checked against changes of the Google Visualisation API, typos in the vignette have been ironed out (thanks to Pat Burns for pointing them out), a new section on dealing with apostrophes in column names has been added and the example in the section "Setting options" has been reviewed. For more details and demos check out the project site.

New Feature

Additionally a new visualisation function has been added: gvisBubbleChart, which provides an interface to the bubble chart of the Google Visualisation API. You could think of it as a static version of the motion chart. Here are some examples, followed by a motion chart.

Loading ...

library(googleVis)
data(Fruits)

P <- lapply(2008:2010, function(x)
      gvisBubbleChart(subset(Fruits, Year %in% x), idvar="Fruit",
       xvar="Sales", yvar="Expenses",
       colorvar="Location", sizevar="Profit",
       options=list(width=400, height=300, 
        colors='["#B2EE2C", "#3F4FFF"]',
        title=paste("Fruit data ",x,", bubble size reflects profit", sep=""),
        sizeAxis="{minValue: 0,  maxSize: 12}",
        vAxis=paste("{title: 'Expenses', viewWindow:{min:65, max:95},",
             "baselineColor:'#EEEEEE', gridlines:{color:'#EEEEEE'}}"),
        hAxis=paste("{title: 'Sales', viewWindow:{min:70, max:115},",
             "baselineColor:'#EEEEEE', gridlines:{color:'#EEEEEE'}}")
             ))
            )
bubbleCharts <- gvisMerge(P[[1]], gvisMerge(P[[2]], P[[3]]))

M <- gvisMotionChart(Fruits, "Fruit", "Year", 
                     options=list(width=430, height=360))
plot(gvisMerge(bubbleCharts, M))

R is the easiest language to speak badly

2012-02-01T07:53:00.000Z

I am amazed by the number of comments I received on my recent blog entry about "by", "apply" and friends. I had started my post by pointing out that R is a language. Well indeed, I have come to the conclusion, that it is a language with lots of irregular expressions and dialects. It feels a bit like German or French where you have to learn and memorise the different articles. The Germans have three singular definite articles: der (male), die (female) and das (neutral), the French have two: le (male) and la (female). Of course there is no mapping between them, and how do you explain that a girl in German is neutral (das Mädchen), while manhood is female (die Männlichkeit)?

Back to R. As I found out, there are lots of different ways to calculate the means on subsets of data. I begin to wonder, why so many different interfaces and functions have been developed over the years, and also why I didn't use the aggregate function more often in the past?

Can we blame internet search engines? Why should I learn a programming language properly, when I can find approximate answers to my problem online. I may not end up with the best answer, but with something which will work after all: Don't know why, but it works.

And sometimes the help files can be more difficult to understand than the code in the examples. Hence, I end up playing around with the example code until it works, and only then I try to figure out how it works. That was my experience with reshape.

Maybe this is a bit harsh. It is always up to the individual to improve his language skills, but you can get drunk in a pub as well, by only being able to order beer. I think it was George Bernard Shaw, who said: "R is the easiest language to speak badly." No, actually he said: "English is the easiest language to speak badly." Maybe that explains the success of English and R?

Reading helps. More and more books have been published on R over the last years, and not only in English. But which should you pick? Xi'an's review on the Art of R Programming suggests that it might be a good start.

Back to aggregate. Has anyone noticed, that the formula interface of aggregate is different to summaryBy?

aggregate(cbind(Sepal.Width, Petal.Width) ~ Species, data=iris, FUN=mean)
     Species Sepal.Width Petal.Width
1     setosa       3.428       0.246
2 versicolor       2.770       1.326
3  virginica       2.974       2.026

versus


library(doBy)
summaryBy(Sepal.Width + Petal.Width ~ Species, data=iris, FUN=mean)
     Species Sepal.Width.mean Petal.Width.mean
1     setosa            3.428            0.246
2 versicolor            2.770            1.326
3  virginica            2.974            2.026

And another slightly more complex example:

aggregate(cbind(ncases, ncontrols) ~ alcgp + tobgp, data = esoph, FUN=sum)
summaryBy(ncases + ncontrols ~ alcgp + tobgp, data = esoph, FUN=sum)

Say it in R with "by", "apply" and friends

2012-01-28T14:43:00.000Z

Iris versicolor (Source: Wikipedia)

R is a language, as Luis Apiolaza pointed out in his recent post. This is absolutely true, and learning a programming language is not much different from learning a foreign language. It takes time and a lot of practice to be proficient in it. I started using R when I moved to the UK and I wonder, if I have a better understanding of English or R by now.

Languages are full of surprises, in particular for non-native speakers. The other day I learned that there is courtesy and curtsey. Both words sounded very similar to me, but of course created some laughter when I mixed them up in an email.

With languages you can get into habits of using certain words and phrases, but sometimes you see or hear something, which shakes you up again. So did the following two lines in R with me:


f <- function(x) x^2
sapply(1:10, f)
[1]   1   4   9  16  25  36  49  64  81 100

It reminded me of the phrase that everything is a list in R. It showed me again how easily a for loop can be turned into a statement using the apply family of functions and how little I know about all the subtleties of R.

I remember how happy I felt, when I finally understood the by function in R. I started to use it all the time, closing my eyes on aggregate and the apply functions family. Here is an example where we calculate the means of the various measurements by species of the famous iris data set using by.

by


do.call("rbind", as.list(
  by(iris, list(Species=iris$Species), function(x){
    y <- subset(x, select= -Species)
    apply(y, 2, mean)
  }
)))

           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa            5.006       3.428        1.462       0.246
versicolor        5.936       2.770        4.260       1.326
virginica         6.588       2.974        5.552       2.026

Now let's find alternative ways of expressing ourselves, using other words/functions of the R language, such as aggregate, apply, sapply, tapply, data.table, ddply, sqldf, and summaryBy.

aggregate

The aggregate function splits the data into subsets and computes summary statistics for each of them. The output of aggregate is a data.frame, including a column for species.


iris.x <- subset(iris, select= -Species)
iris.s <- subset(iris, select= Species)
aggregate(iris.x, iris.s, mean)

     Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1     setosa        5.006       3.428        1.462       0.246
2 versicolor        5.936       2.770        4.260       1.326
3  virginica        6.588       2.974        5.552       2.026

Addition: As John Christie points out in the comments, aggregate has also a formula interface, which simplifies the call to:

aggregate( . ~ Species, iris, mean)

apply and tapply

The combination of tapply and apply achieves a similar result, but this time the output is a matrix and hence we lose the column with species. The species are now the row names.

apply(iris.x, 2, function(x) tapply(x, iris.s, mean))

           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa            5.006       3.428        1.462       0.246
versicolor        5.936       2.770        4.260       1.326
virginica         6.588       2.974        5.552       2.026

split and apply

Here we split the data first into subsets for each of the species and calculate then the mean for each column in the subset. The output is a matrix again, but transposed.

sapply(split(iris.x, iris.s), function(x) apply(x, 2, mean))

             setosa versicolor virginica
Sepal.Length  5.006      5.936     6.588
Sepal.Width   3.428      2.770     2.974
Petal.Length  1.462      4.260     5.552
Petal.Width   0.246      1.326     2.026

ddply

Hadley Wickham's plyr package provides tools for splitting, applying and combining data. The function ddply is similar to the by function, but it returns a data.frame instead of a by list and maintains the column for the species.

library(plyr)
ddply(iris, "Species", function(x){
    y <- subset(x, select= -Species)
    apply(y, 2, mean)
  })

     Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1     setosa        5.006       3.428        1.462       0.246
2 versicolor        5.936       2.770        4.260       1.326
3  virginica        6.588       2.974        5.552       2.026

Addition: Sean mentions in the comments an alternative, using the colMeans function, while Andrew reminds us of the reshape package with its functions melt and cast.

ddply(iris, "Species", function(x) colMeans(subset(x, select= -Species)))
## or
ddply(iris, "Species", colwise(mean)) 
## same output as above
library(reshape)
cast(melt(iris, id.vars='Species'),formula=Species ~ variable,mean)
## same output as above

summaryBy

The summaryBy function of the doBy package by Søren Højsgaard and Ulrich Halekoh has a very intuitive interface, using formulas.

library(doBy)
summaryBy(Sepal.Length + Sepal.Width + Petal.Length + Petal.Width ~ Species, data=iris, FUN=mean)

     Species Sepal.Length.mean Sepal.Width.mean Petal.Length.mean Petal.Width.mean
1     setosa             5.006            3.428             1.462            0.246
2 versicolor             5.936            2.770             4.260            1.326
3  virginica             6.588            2.974             5.552            2.026

sqldf

If you are fluent in SQL, then the sqldf library by Gabor Grothendieck might be the one for you.

library(sqldf)
sqldf("select Species, avg(Sepal_Length), avg(Sepal_Width), 
    avg(Petal_Length), avg(Petal_Width) from iris 
    group by Species")

     Species avg(Sepal_Length) avg(Sepal_Width) avg(Petal_Length) avg(Petal_Width)
1     setosa             5.006            3.428             1.462            0.246
2 versicolor             5.936            2.770             4.260            1.326
3  virginica             6.588            2.974             5.552            2.026

data.table

The data.table package by M Dowle, T Short and S Lianoglou is the real rock star to me. It provides an elegant and fast way to complete our task. The statement reads in plain English from right to left: take columns 1 to 4, split them by the factor in column "Species" and calculate on the sub data (.SD) the means.

library(data.table)
iris.dt <- data.table(iris)
iris.dt[,lapply(.SD,mean),by="Species",.SDcols=1:4]

        Species Sepal.Length Sepal.Width Petal.Length Petal.Width
[1,]     setosa        5.006       3.428        1.462       0.246
[2,] versicolor        5.936       2.770        4.260       1.326
[3,]  virginica        6.588       2.974        5.552       2.026

apply

I should mention that R provides the iris data set also in an array form. The third dimension of the iris3 array holds the species information. Therefore we can use the apply function again, we go down the third and then the second dimension to calculate the means.

apply(iris3, c(3,2), mean)

           Sepal L. Sepal W. Petal L. Petal W.
Setosa        5.006    3.428    1.462    0.246
Versicolor    5.936    2.770    4.260    1.326
Virginica     6.588    2.974    5.552    2.026

Conclusion

Many roads lead to Rome, and there are endless ways of explaining how to get there. I only showed a few I know of, and I am curious to hear yours.

As a matter of courtesy I should mention the unknownR package by Matthew Dowle. It helps you to discover what you don't know that you don't know in R. Thus, it can help to build your R vocabulary.

Of course there is a key difference between R and English. R tells me right away when I make a mistake. Human readers are far more forgiving, but please do point out to me where I made mistakes. I am still hopeful that I can improve, but I need your help.

R code

The R code of the examples is available on github. For more examples on the apply family see also Neil Saunders' post.

Credit rating by country

2012-01-17T22:16:00.000Z

The financial crisis has put a lot of pressure on countries' long-term foreign currency credit ratings, with France recently being downgraded by S&P. Wikipedia provides a list of countries by credit ratings as report by US rating agencies S&P, Fitch, Moody's and Dagong, a Chinese rating agency.

So, was does the world look like today through the eyes of those rating agencies?

We use the R packages XML and googleVis to read and display the data from Wikipedia with just a few lines.

Loading ...
The geo maps show the data as sourced from Wikipedia on 17 January 2012.

Is it a surprise that the views of the US and Chinese rating agencies differ on their home countries?

Here is the R code to replicate the geo maps:

Managing change

2012-01-14T22:23:00.000Z

Why the old and the new need to share time together

It takes time to appreciate the new. Even if the new is much better than the old. It is easy to forget when you yourself created the exciting new.

At the end of August 2011 Google announced a new Blogger interface. The new interface offered about the same functionality, but had a different look and feel. At first I was reluctant to use it. I just wanted to get my job done. I knew the old interface, so why disturbe me with something else which does just the same?

Old Blogger interface

I was relieved to see that there was a button, which allowed me to switch back to the old interface. So I continued to use the old interface, although occasionally Google would offer me the new interface after login. Eventually I started to play around with the new interface. At times when I was bored and looking for a good excuse to waste my time. After a while I started to like the new interface and used it more often and then I stopped switching back to the old when Google offered me the new on login. It appeared cleaner, more contemporary and organised than the old. Now I have switched completely to the new interface. After a four months transition period!

New Blogger interface

There is something to be learned from this experience

Suppose you write analysis reports for others, e.g. colleagues, clients, friends, students, etc. and you have a new idea, which will make your report so much better. Maybe a better structure, a new R package, a new plot, something really sophisticated. You are excited about this, you implement it and send your report out to your readers.

Then: No feedback.

Appreciating change takes time.

Here is a suggestion: Implement your next great idea to improve your analysis report. Just because it is so much fun and your are really excited about the new. However, continue to produce the old style report and send it out with the new. Hopefully your excitement will sparkle over to the reader. But acknowledge that not all of them will have time to embrace the new immediately. They read the report for a specific purpose, e.g. table X on page Y. Hence give your reader the space to find his own time to review the new report, eventually he will listen and hopefully switch. But do listen to the reader as well and value his feedback. Easier said than done.

There is another aspect to it. Sending out the new with the old allows you to fail. You have a safety net. Some of your ideas might actually be not that great after all. I believe, we can only learn and succeed by making lots of small failures.

PS:

Producing the old and the new at the same time is of course only possible, if you have a slick automated reporting system. You find plenty of ideas on R-Bloggers.
I wonder, what would have happened if the transition period between the old national coins and notes and the Euro would have been longer than two months in 2002?

Feedback from vignette survey

2012-01-08T21:45:00.000Z

Many thanks to all who participated in the survey about writing R package vignettes.

Following my post last Thursday the responses came in quickly in the evening and all day on Friday. Since Saturday the response rate has been decreasing constantly and I think it is time for a summary based on the 56 responses received.

Summary - How to write a good vignette

Length: Trust yourself, but aim for about 20 pages.
Language: Don't use language which assumes that the reader is an R and/or subject expert.
Structure: Include at least the following sections:
- Examples
- Introduction
- Case studies
- References
It would be nice to include also sections on:
- Support
- Motivation
- Road map
Examples: Use lots of examples and don't repeat just the examples from the help pages.
Get inspiration from: Rcpp, reshape, plyr, vegan, and see below for more.
Secrets of good vignettes:
- Provide an introduction with a clear purpose of the package.
- Work with case studies, walk the reader through a task from start to finish.
- Demonstrate the non-default arguments of the package functions, highlight why and when you want to change them.
- Write briefly and concisely, but provide reference/footnotes to relevant literature and further help.
- Provide dummy data to play with.
- Discuss limitations.
What else: Potentially split the vignette into several documents, see Rcpp for an example.

About the survey participants

Expertise: Most participants seem to be medium to advanced R users. This is not a surprise - how else would they know about my blog or R-bloggers? Over 50% write their own functions and a further 45% create R packages.
Usage of vignettes: Over 75% use them often (based on ratings with 4/5 or 5/5). Thus, vignettes are worth wile writing. They are often the starting point to learn more about a new package.
Writing experience: About 1/5 of the participants have written vignettes themselves and about another third have published papers, but not written vignettes or help pages.

Conclusions

This was an interesting experiment and I am really pleased with the turn out. The survey was completed by a good mixture of people, with varied experience and it showed the value of vignettes.

My big learning point is to focus more on the non-default arguments of the key package functions. I can see, that this really helps the reader to get a better understanding of the package and will help them to apply it better to their own problems and needs.

I was surprised by the relative high request for a support section in the vignette. Over 60% say that a support section would be at least nice to have. Even 14% think it is required. I wonder, if this means just an email address for questions, a mailing list, or commercial support. I would very much appreciate your feedback, either via email or comments below the post.

Charts

About the survey participants

R code

You can access the R code for all charts via github.

Survey: Writing package vignette

2012-01-05T18:01:00.000Z

I am currently co-writing the vignette for the ChainLadder package and wonder what I should be focusing on. I have co-written the vignette of the googleVis package in the past and based it purely and what I thought would work. So, this is an experiment to find out, if user feedback will help me to write a better vignette. Let's see how it develops. I will make the data available once I have at least 10 submission.

Loading...
Thanks!

Is R turning into an operating system?

2011-12-31T14:22:00.001Z

Over the years I convinced my colleagues and IT guys that LaTeX/XeLaTeX is the way forward to produce lots of customer reports with individual data, charts, analysis and text. Success! But of course the operating system in the office is still MS Windows.

With my background in Solaris/Linux/Mac OSX I am still a little bit lost in the Windows world, when I have to do such simple tasks as finding and replacing a string in lots of files. Apparently the acronym is FART (find and replace text).

So, what to you do without your beloved command line tools and admin rights? Eventually you start using R instead. So here is my little work-around for replacing "Merry Christmas" with "Happy New Year" in lots of files:

for( f in filenames ){

  x <- readLines(f)
  y <- gsub( "Merry Christmas", "Happy New Year", x )
  cat(y, file=f, sep="\n")

}

You can find a complete self-contained example on github.

Of course R is not an operating system, but yet it can complement it well, if your other resources are limited.

Last but not least: Happy New Year!

googleVis 0.2.13: new stepped area chart and improved geo charts

2011-12-23T10:06:00.004Z

On 7th December Google published a new version of their Visualisation API. The new version adds a new chart type: Stepped Area Chart and provides improvements to Geo Chart. Now Geo Chart has similar functionality to Geo Map, but while Geo Map requires Flash, Geo Chart doesn't, as it renders SVG/VML graphics. So it also works on your iOS devices.

These new features have been added to the googleVis R package in version 0.2.13, which went live on CRAN a few days ago.

The function gvisSteppedAreaChart works very much in the same way as gvisAreaChart. Here is a little example:

 library(googleVis)
df <- data.frame(country=c("US", "GB", "BR"), 
                 val1=c(1,3,4), val2=c(23,12,32))
SteppedArea <- gvisSteppedAreaChart(df, xvar="country", 
                                        yvar=c("val1", "val2"),
                           options=list(isStacked=TRUE, 
                                        width=400, height=150))
plot(SteppedArea2)

The interface to gvisGeoChart changed slightly to take into account the new version of Geo Chart by Google. The argument numvar has been renamed to colorvar and a new argument sizevar has been added. This allows you to set the size and colour of the bubbles in displayMode='markers' depending on columns in your data frame. Further, you can set far more options than you could before, in particular you can set not only the region, but also the resolution of your map. Although more granular maps are not available for all countries, for more details see the Google documentation.

Here are two examples, plotting the test data CityPopularity with a Geo Chart.The first plot shows the popularity of New York, Boston, Miami, Chicago, Los Angeles and Houston on the US map, with the resolution set to 'metros' and region set to 'US'. The Google Map API makes the correct assumption about which cities we mean.

library(googleVis> ## requires googleVis version >= 0.2.13
gcus <- gvisGeoChart(CityPopularity, 
        locationvar="City", colorvar="Popularity", 
          options=list(displayMode="markers", 
                       region="US", resolution="metros"), 
         chartid="GeoChart_US")
plot(gcus)

In the second example we set the region to 'US-TX', therefore Google will look for cities with the same names in Texas. And what a surprise, there are cities/towns named Chicago, Los Angeles, Miami, Boston and of course Houston in Texas.

gctx <- gvisGeoChart(CityPopularity, 
        locationvar="City", colorvar="Popularity", 
          options=list(displayMode="markers", 
                       region="US-TX", resolution="metros"), 
         chartid="GeoChart_TX")
plot(gctx)

With the new version of the Visualisation API Google introduced also the concept of DataTable Roles. This is an interesting idea, as it allows you to add context to the data, similar to the approach used with annotated time lines. Google classifies the DataTable Roles still experimental, but it is a space to watch and ideas on how this could be translated into R will be much appreciated.

And now the news of the googleVis package since version 0.2.10:

Version 0.2.13 [2011-12-19]
==========================

Changes

    o The list of arguments for gvisGeoChart changed:
      - the argument 'numvar' has been renamed to 'colorvar' to
        reflect the updated Google API. Additionally gvisGeoChart
     gained a new argument 'sizevar'.
    o Updated googleVis vignette with a section on using googleVis 
      output in presentations  
    o Renamed demo EventListner to EventListener

NEW FEATURES

    o Google published a new version of their Visualisation API on 7
      December 2011. Some of the new features have been implemented
      into googleVis already:
      - New stepped area chart function gvisSteppedAreaChart
      - gvisGeoChart has a new marker mode, similar to the mode in
       gvisGeoMap. See example(gvisGeoChart) for the new
        functionalities.

Version 0.2.12 [2011-12-07]
==========================

Bug Fixes

    o gvisMotionChart didn't display data with special characters,
      e.g. spaces, &, %, in column names correctly. 
      Thanks to Alexander Holcroft for reporting this issue.

Version 0.2.11 [2011-11-16]
==========================

Changes

   o Updated vignette and documentation with instructions on changing
     the Flash security settings to display Flash charts locally. 
     Thanks to Tony Breyal.
   o New example to plot weekly data with gvisMotionChart
   o Removed local copies of gadget files to reduce package file
     size. A local copy of the R script to generate the original gadget
     files is still included in inst/gadgets 

Version 0.2.10 [2011-09-24]
==========================

Changes

   o Updated section 'Using googleVis output with Google Sites,
     Blogger, etc.' vignette

   o Updated example for gvisMotionChart, showing how the initial
     chart setting can be changed, e.g to display a line chart.

   o New example for gvisAnnotatedTimeLine, showing how to shade
     areas. Thanks to Mike Silberbauer for providing the initial code.    
   
NEW FEATURES
 
    o New demo WorldBank. It demonstrates how country level data can
      be accessed from the World Bank via their API and displayed with a
      Motion Chart. Inspired by Google's Public Data Explorer, see
      http://www.google.com/publicdata/home

Data is the new gold

2011-12-13T22:41:00.006Z

We need more data journalism. How else will we find the nuggets of data and information worth reading?

Life should become easier for data journalists, as the Guardian, one of the data journalism pioneers, points out in this article about the new open data initiative of the European Union (EU). The aims of the EU's open data strategy are bold. Data is seen as the new gold of the digital age. The EU is estimating that public data is already generating economic value of €32bn each year, with growth potential to €70bn, if more data will be made available. Here is the link to the press statement, which I highly recommend reading:

EUROPA - Press Releases - Neelie Kroes Vice-President of the European Commission responsible for the Digital Agenda, Data is the new gold, Opening Remarks, Press Conference on Open Data Strategy Brussels, 12th December 2011

I am particularly impressed that the EU even aims to harmonise the way data will be published by the various bodies. We know that working with data, open or proprietary, often means spending a lot of time on cleaning, reshaping and transforming it, in order to join it with other sources and to make sense out of it.

Data standards would really help in this respect. And the EU is pushing this as well. I can observe this in the insurance industry already, where new European regulatory requirements (Solvency II) force companies to increase their data management capabilities. This is often a huge investment and has to be seen as a long term project.

Although the press statement doesn't mention anything about open source software projects, I think that they are essential for unfolding the full potential of open data.

Open source projects like R provide a platform to share new ideas. I'd say that R, but equally other languages as well, provide interfaces between minds and hands. Packages, libraries, etc. make it possible to spread ideas and knowledge. Having access to scientific papers is great but being able to test the ideas in practice accelerates the time it takes to embed new developments from academia into the business world.

LondonR, 6 December 2011

2011-12-10T08:54:00.001Z

The London R user group met again last Wednesday at the Shooting Star pub. And it was busy. More than 80 people had turned up. Was it the free beer and food, sponsored by Mango, which attracted the folks or the speakers? Or the venue? James Long, who organises the Chicago R user group meetings and who gave gave the first talk that night, noted that to his knowledge only the London and Chicago R users would meet in a pub.

However, it were the speakers and their talks which attracted me:

James Long: Easy Parallel Stochastic Simulations using Amazon's EC2 & Segue
Chibisi Chima-Okereke: Actuarial Pricing Using General Linear Models In R
Richard Saldanha: Practical Optimisation

You will notice that this London R meeting had a theme around risk pricing. James talked about reinsurance pricing using R in the cloud, while Chibisi focused more on personal lines insurance with generalised linear models and Richard came from the angle of investment management and portfolio optimisation.

Fitting distributions with R

2011-12-01T17:40:00.003Z

Fitting distribution with R is something I have to do once in a while.

A good starting point to learn more about distribution fitting with R is Vito Ricci's tutorial on CRAN. I also find the vignettes of the actuar and fitdistrplus package a good read. I haven't looked into the recently published Handbook of fitting statistical distributions with R, by Z. Karian and E.J. Dudewicz, but it might be worthwhile in certain cases, see Xi'An's review. A more comprehensive overview of the various R packages is given by the CRAN Task View: Probability Distributions, maintained by Christophe Dutang.

How do you decide which distribution might be a good starting point?

I came across the paper Probabilistic approaches to risk by Aswath Damodaran. In Appendix 6.1 Aswath discusses the key characteristics of the most common distributions and in Figure 6A.15 he provides us with a decision tree diagram for choosing a distribution:

Figure 6A.15 from Probabilistic approaches to risk by Aswath Damodaran

JD Long points in his blog entry about Fitting distribution X to data from distribution Y to the Clickable diagram of distribution relationships by John Cook. With those two charts I find it not too difficult anymore to find a reasonable starting point.

Once I have decided which distribution might be a good fit I start usually with the fitdistr function of the MASS package. However, since I discovered the fitdistrplus package I have become very fond of the fitdist function, as it comes with a wonderful plot method. It plots an empirical histogram with a theoretical density curve, a QQ and PP-plot and the empirical cumulative distribution with the theoretical distribution. Further, the package provides also goodness of fit tests via gofstat.

Suppose you have only 50 data points, of which you believe that they follow a log-normal distribution. How much variance can we expect? Well, let's experiment. We draw 50 random numbers from a log-normal distribution, fit the distribution to the sample data and repeat the exercise 50 times and plot the results using the plot function of the fitdistrplus package.

You will notice quite a big variance in the results. For some samples other distributions, e.g. logistic, could provide a better fit. You might argue that 50 data points is not a lot of data, but in real life it often is, and hence this little example already shows us that fitting a distribution to data is not just about applying a algorithm, but requires a sound understanding of the process which generated the data as well.

Interactive presentations with deck.js

2011-11-20T18:39:00.002Z

Data analysis is often an iterative and interactive process. However, when I present about this subject, I feel often limited by the presentation software I use. It doesn't matter if I use LaTeX/PDF, PowerPoint or Keynote. In all cases it is either very difficult or impossible to include interactive charts, such as Flash or SVG charts. As a result I have to switch between various applications during the talk. This can be fun, but quite often it is not.

The other day I came across a presentation by Christopher Gandrud. Christopher had used deck.js, a JavaScript library for building HTML presentations by Caleb Troughton.

This looked like an interesting approach to me and fortunately the learning curve was not too steep, although I am by no means an html or JavaScript expert. So I created my first deck.js presentation based on the content of previous googleVis presentations. For the first time I can embed videos, Flash and SVG charts without using lots of different apps. I am actually quite pleased by the result, see here: Getting started with googleVis

Now imagine a presentation hosted on a server with R installed! You could combine your slides with R using one of the following packages R.rsp, brew, Rook, etc and run live demos, without opening a console.

Stochastic reserving with R: ChainLadder 0.1.5-1 released

2011-11-12T21:55:00.015Z

Today we published version 0.1.5-1 of the ChainLadder package for R. It provides methods which are typically used in insurance claims reserving to forecast future claims payments.

Claims development and chain-ladder forecast of the RAA data set using the Mack method

The package started out of presentations given at the Stochastic Reserving Seminar at the Institute of Actuaries in 2007, 2008 and 2010, followed by talks at CAS meetings in 2008 and 2010.

Initially the package came with implementations of the Mack-, Munich- and Bootstrap Chain-Ladder methods. Since version 0.1.3-3 it also provides general multivariate chain ladder models by Wayne Zhang. Version 0.1.4-0 introduced new functions on loss development factor fitting and Cape Cod by Daniel Murphy following a paper by David Clark. Version 0.1.5-0 has added loss reserving models within the generalized linear model framework following a paper by England P. and Verrall R. (1999) implemented by Wayne Zhang.

For more details see the project web site: http://code.google.com/p/chainladder/ and an early blog entry about R in the insurance industry.

Changes in version 0.1.5-1:

Internal changes to plot.MackChainLadder to pass new checks introduced by R 2.14.0.
Commented out unnecessary creation of 'io' matrix in ClarkCapeCod function. Allows for analysis of very large matrices for CapeCod without running out of RAM. 'io' matrix is an integral part of ClarkLDF, and so remains in that function.
plot.clark method
- Removed "conclusion" stated in QQplot of clark methods.
- Restore 'par' settings upon exit
- Slight change to the title
Reduced the minimum 'theta' boundary for weibull growth function
Added warnings to as.triangle if origin or dev. period are not numeric

Here is a little example using the googleVis package to display the RAA claims development triangle:

library(ChainLadder)
library(googleVis)
data(RAA) # example data set of the ChainLadder package
class(RAA) <- "matrix" # change the class from triangle to matrix
df <- as.data.frame(t(RAA)) # coerce triangle into a data.frame
names(df) <- 1981 : 1990
df$dev <-  1:10
plot(gvisLineChart(df, "dev", options=list(gvis.editor="Edit me!", hAxis.title="dev. period")))

Installing R 2.14.0 on an iBook G4 running Mac OS 10.4.11

2011-10-30T08:00:00.004Z

My 12" iBook G4 is celebrating its 8th birthday today! Time for a little present. How about R 2.14.0?

The iBook is still in daily use, mostly for browsing the web, writing e-mails and this blog; and I still use it for R as well. For a long time it run R 2.10.1, the last PowerPC binary version available on CRAN for Mac OS 10.4.11 (Tiger).

But, R 2.10.1 is a bit dated by now and for the development of my googleVis package I require at least R 2.11.0. So I decided to try to install the most recent version from source, using Xcode 2.5 and TeXLive-2008.

R 2.14.0 is expected to be released on Monday (31st October 2011). The pre-release version is already available on CRAN. I assume that the pre-release version is pretty close to the final version of R 2.14.0, so why wait?

It was actually surprisingly easy to compile the command line version of R from sources. The GUI would be a nice to have, but I am perfectly happy to run R via the Terminal, xterm and Emacs. However, it shouldn't be a surprise that running configure, make, make install on a 800 Mhz G4 with 640MB memory does take its time.

Below you will find the building details. Please feel free to get in touch with me, if you would like access to my Apple Disk Image (dmg) file. You find my e-mail address in the maintainer field of the googleVis package.

Building R from source on Mac OS 10.4 with Xcode 2.5 (gcc-4.0.1)

Before you start, make sure you have all the Apple Developer Tools installed. I have Xcode installed in /Developer/Applications.

From the pre-release directory on CRAN I downloaded the file R-rc_2011-10-28_r57465.tar.gz.

The first steps were to extract the archive and to run the configure scripts to build the various Makefiles. To do this, I opened the Terminal programme (it's in the Utilities folder of Applications), changed into the directory in which I stored the tar.gz-file and typed:

tar xvfz R-rc_2011-10-28_r57465.tar.gz
cd R-rc
./configure

This process took a little while (about 15 minutes) and at the end I received the following statement:

R is now configured for powerpc-apple-darwin8.11.0

  Source directory:          .
  Installation directory:    /Library/Frameworks

  C compiler:                gcc -std=gnu99  -g -O2
  Fortran 77 compiler:       gfortran  -g -O2

  C++ compiler:              g++  -g -O2
  Fortran 90/95 compiler:    gfortran -g -O2
  Obj-C compiler:            gcc -g -O2 -fobjc-exceptions

  Interfaces supported:      X11, aqua, tcltk
  External libraries:        readline, ICU
  Additional capabilities:   NLS
  Options enabled:           framework, shared BLAS, R profiling, Java

  Recommended packages:      yes

With all the relevant Makefiles in place I could start the build process via:

make -j8

Now I had time for a cup of tea, as the build took about one hour. Finally, to finish the installation, I placed the new R version into its place in /Library/Frameworks/ by typing:

sudo make install

Job done. Let's test it:

Grappa:~ Markus$ R

R version 2.14.0 RC (2011-10-28 r57465)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: powerpc-apple-darwin8.11.0 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> for(i in 1:8) print("Happy birthday iBook!")
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"
[1] "Happy birthday iBook!"

The installation of additional packages worked straightforward via install.packages(c("vector of packages")), though it took time, as everything was build from sources. This it what it looks like on my iBook G4 today:

> installed.packages()[,"Version"]
 ChainLadder GillespieSSA        Hmisc     ISOcodes   KernSmooth         MASS 
   "0.1.5-0"      "0.5-4"      "3.8-3" "2011.07.31"     "2.23-6"     "7.3-16" 
      Matrix  R.methodsS3         R.oo        R.rsp      R.utils RColorBrewer 
     "1.0-1"      "1.2.1"      "1.8.2"      "0.6.2"      "1.8.5"      "1.0-5" 
       RCurl      RJSONIO        RUnit         Rook          XML       actuar 
    "1.6-10"     "0.96-0"     "0.4.26"      "1.0-2"      "3.4-3"      "1.1-2" 
        base       bitops         boot         brew          car        class 
    "2.14.0"    "1.0-4.1"      "1.3-3"      "1.0-6"     "2.0-11"      "7.3-3" 
     cluster         coda    codetools         coin   colorspace     compiler 
    "1.14.1"     "0.14-4"      "0.2-8"     "1.0-20"      "1.1-0"     "2.14.0" 
  data.table     datasets       digest    flexclust      foreign          gam 
     "1.7.1"     "2.14.0"      "0.5.1"      "1.3-2"     "0.8-46"     "1.04.1" 
     ggplot2    googleVis    grDevices     graphics         grid    iterators 
     "0.8.9"     "0.2.10"     "2.14.0"     "2.14.0"     "2.14.0"      "1.0.5" 
   itertools      lattice       lmtest       mclust      methods         mgcv 
     "0.1-1"     "0.20-0"     "0.9-29"     "3.4.10"     "2.14.0"      "1.7-9" 
  modeltools      mvtnorm         nlme         nnet     parallel        party 
    "0.2-18"   "0.9-9991"    "3.1-102"      "7.3-1"     "2.14.0"  "0.9-99994" 
        plyr        proto         pscl      reshape        rpart     sandwich 
       "1.6"    "0.3-9.2"     "1.04.1"      "0.8.4"     "3.1-50"      "2.2-8" 
     spatial      splines      statmod        stats       stats4  strucchange 
     "7.3-3"     "2.14.0"     "1.4.13"     "2.14.0"     "2.14.0"      "1.4-6" 
    survival    systemfit        tcltk        tools        utils          vcd 
   "2.36-10"      "1.1-8"     "2.14.0"     "2.14.0"     "2.14.0"     "1.2-12" 
         zoo 
     "1.7-5"

Using Sweave with XeLaTeX

2011-10-23T17:35:00.007+01:00

Using R with LaTeX via Sweave is a great way to create reproducible output. However, using specific fonts, e.g. your corporate fonts, can be painful with pdflatex. Over the last few weeks I have fallen in love with the TeX format XeLaTeX and its XeTeX engine.

With XeLaTeX I had to overcome some hurdles, which I would like to share here:

attaching files,
trimming and clipping images,
learning how to use the tikzDevice package.

The Sweave file of the above document is attached to the PDF-file itself, but you can also find it on github: SweaveXeLaTeXExample.Rnw.

R related books: Traditional vs online publishing

2011-10-12T08:23:00.010+01:00

How many R related books have been published so far? Who is the most popular publisher? How many other manuals, tutorials and books have been published online? Let's find out.

A few years ago I used the publication list on r-project.org as an argument with the IT department that R is an established statistical programming language and that they should allow me to install it on my PC. I believe at the time there were about 20 R related books available.

A recent post on Recology pointed me to a talk given by Ed Goodwin at the Houston R user group meeting about regular expressions in R, something I always wanted to learn properly, but never got around to do.

So let's see, if we can manage to extract the information of published R books and texts from r-project.org, with what we learned from Ed about regular expressions in R.

We will start analysing the bib-file on r-project.org by publisher and then move on to look more closely at the number of titles published over time, including the self-published PDF-files on CRAN.

We read the bib-file into R using the readLines function and start analysing the data with regular expressions. The function regexpr allows us to find for each line the character place where "publisher =" starts, or it will return -1 otherwise, if no entry is found. We continue with the R function strsplit to cut the strings into sub-component for further analysis:

bibfile <- readLines("http://www.r-project.org/doc/bib/R-books.bib")
pub.start.pos <- regexpr("publisher =", bibfile, perl=TRUE)
pub.lines <- which( pub.start.pos > 0 ) 
pub.split <- strsplit(bibfile[pub.lines], "[ =,]", perl=TRUE) 
publishers <- sapply(pub.split, function(x) paste(x[-c(1:5)]))
publishers <- gsub("[{}\"),\\]", "", publishers)
publishers <- gsub("c\\(", "", publishers)
 
s=c("Springer", "Wiley", "Sage", "Chapman & Hall", "CRC press", "Servicio")
r=c("Springer", "Wiley", "Sage", "Chapman & Hall/CRC", 
    "Chapman & Hall/CRC","Universidad de La Rioja")
 
for(i in seq(along=s)){
  publishers[regexpr(s[i], publishers, ignore.case=TRUE) > 0] <- r[i]
}
 
pubs <- table(publishers)
pubs <- pubs[order(pubs)]
opar <- par(mar = c(4, 15, 4, 2))
barplot(pubs, horiz=TRUE, las=1, col="#DC3912",
  xlab=format(Sys.time(), "As at %d %b %Y"), 
  border=NA, main="Number of R books\nby publisher")
par(opar)

We note that Springer is by far the most popular publisher for R related books. Thus, if you are looking for a specific topic around R your safest bet would be to check out Springer's portfolio.

However, although currently Springer is the publisher with the highest appetite for R, you may be able to find the information free online on r-project.org, in particular if your are looking for a tutorial like document/book.

Hence we want to compare the number of published R books in a traditional way, versus the PDF-files contributed online on CRAN: http://cran.r-project.org/doc/contrib/. CRAN is in this respect also a bit of a publisher, as I assume that the guys behind CRAN have some kind of a filtering and QA process.

We use the XML package to read the online directory content of the contributed books into R to get a better understanding of the published PDF-files. In a similar approach as above, we analyse the R books and PDF files published by year. Please find the R code below the charts.

From the two charts below we can see that over the years more and more R texts have been made available, illustrating the increased interest in R. The first chart shows the number of books/documents published in each year, while the second chart shows the same data in a cumulative way.

Today there are 206 R related books available either on CRAN or via your bookshop. Of the 206 texts 113 are published in the traditional sense with a publishing house, and the number is still growing. However the growth has slowed down a bit over the recent years, after a peak of 26 new books in 2009.

From the second chart I can see that I must have had the discussion on R I mentioned earlier around 2004 - 2005.

&amp;lt;br /&amp;gt;

## Continue from the above R code
year.start.pos <- regexpr("year =", bibfile, perl=TRUE)
year.lines <- which( year.start.pos > 0 ) 
year.split <- strsplit(bibfile[year.lines], "[ =,]", perl=TRUE) 
Pub.year <- as.numeric(sapply(year.split, "[[", 6))
tradPub <- as.data.frame(table(Pub.year))
names(tradPub) <- c("year", "traditional book")
 
## Get information about the online published PDF files
library(XML)
webPub <-
readHTMLTable(readLines("http://cran.r-project.org/doc/contrib/"))[[1]]
## look only at the PDF files
pdfs <- regexpr(".pdf", as.character(webPub[,2]))
 
webPub.modified <- webPub[pdfs > 0, 3]
webPub.modified <- strsplit(na.omit(as.character(webPub.modified)), "[ -]")
webPub.year <- as.numeric(sapply(webPub.modified, "[[", 3))
 
webPub <- as.data.frame(table(webPub.year))
names(webPub) <- c("year", "online PDF")
## Merge information about the online and traditional books
totalPub <- merge(tradPub, webPub, all=TRUE)
totalPub[is.na(totalPub)] <- 0 
## Calculate the cumulative statistics
cumPub <- data.frame(year=totalPub$year, apply(totalPub[,2:3], 2, cumsum ))
names(cumPub) <- c("year", "traditional book",
           "online PDF")
## We use the googleVis package to create area charts
library(googleVis)
cumPlot <- gvisAreaChart(cumPub, "year", 
      options=list(title="Number of R related books available", 
    isStacked=TRUE))
 
incPlot <- gvisAreaChart(totalPub, "year", 
          options=list(title="Number of R
       related books published",
       isStacked=FALSE))
 
inccumPlot <- gvisMerge(incPlot, cumPlot)
 
plot(inccumPlot)

Created by Pretty R at inside-R.org

Setting the initial view of a motion chart in R

2011-09-30T08:00:00.004+01:00

Following on from my article about accessing and plotting World Bank data with R I want to talk about how to change the initial view of a motion chart.

Over the last couple of weeks I have been asked a view times how to do this. For instance Stephen O'Grady wanted to create a motion chart, which shows initially a line chart, rather than a bubble chart.

Changing the initial settings of a motion chart is actually quite easy, if you know how to. The trick is to use the state argument in the list of options of gvisMotionChart.

As a case study we will come back to the World Bank data set and try to do some homework givem by Duncan Temple Lang in his course on introduction to statistical computing course by. Duncan asked his students to query the World Bank data base to create a plot, which would show the number of internet users per 1000 in Africa over time as line chart. Further, he would like to see a legend next to the chart to identify which country is which and tooltips for each curve, so we can identify the country.

A motion chart, displayed as a line chart, would do the trick.

Okay, getting the data is easy, thanks to the WDI package, or via a direct download, and so it is to create a motion chart with bubbles. Interactively we can change the bubble chart into a line chart, we can select some countries and change the y-axis to log-scale. However, when we reload the page we are back to square one: a bubble chart. So the idea is to pass the changed chart setting on to the initial plot. You find those settings, of the current view, as a string in the advanced tab of the settings window. Click on the wrench symbol in the bottom right hand corner of a motion chart to access this window.

Screen shot the settings window of a motion chart

Copy this string and paste it into the state argument of the options list.

Here is an example, where I pre-selected Sierra Leone and Seychelles (countries with the lowest and highest number of internet users) together with Africa, North Africa and Sub-Saharan Africa (all income levels). You find the R code below to replicate the plot.

What does the data tell you? Play around with the graph, e.g. change it to a column graph, deselect all countries and change the y-axes to linear again, and hit the play button. How could we improve the plot?

For more details about the googleVis package see the vignette, and for more details about motion charts and its option see also the Google Visualisation API documentation.

## Based on ideas and code from Duncan Temple Lang, see:
## http://eeyore.ucdavis.edu/stat242/Homeworks/XML.html
library(googleVis)
 
url = "http://api.worldbank.org/countries/all/indicators/IT.NET.USER.P3?date=1990:2010&format=json&per_page=12000"
iu = fromJSON(url)
internetUsers = data.frame(year = as.numeric(sapply(iu[[2]], "[[", "date")),
  InternetUsersPerThousands = as.numeric(sapply(iu[[2]],
   function(x) ifelse(is.null(x[["value"]]),  NA, x[["value"]]))),
  country = sapply(iu[[2]], function(x)
    x[["country"]]['value']))
 
## Create a line plot with gvisMotionChart
## Set initial state with a few regions selected and a log y-axes
myState <- '
{"yZoomedIn":false,"yZoomedDataMin":0,"xAxisOption":"_TIME",
 "xZoomedDataMax":1230768000000,"time":"2009","playDuration":15000,
 "yAxisOption":"2","yLambda":0,"duration":{"timeUnit":"Y","multiplier":1},
 "xZoomedDataMin":631152000000,"uniColorForNonSelected":false,
 "showTrails":false,"orderedByY":false,"nonSelectedAlpha":0,
 "orderedByX":false,"sizeOption":"_UNISIZE","xLambda":1,
 "colorOption":"_UNIQUE_COLOR",
 "iconKeySettings":[{"key":{"dim0":"Sub-Saharan Africa (all income levels)"}},
                    {"key":{"dim0":"North Africa"}},
                    {"key":{"dim0":"Seychelles"}},
                    {"key":{"dim0":"Sierra Leone"}},
                    {"key":{"dim0":"Africa"}}],
 "dimensions":{"iconDimensions":["dim0"]},
 "yZoomedDataMax":220,"xZoomedIn":false,"iconType":"LINE"}
'
M <- gvisMotionChart(internetUsers, "country", "year", 
  options=list(width=650, height=350, state=myState))
plot(M)
 
## Session Info
## sessionInfo()
## R version 2.13.1 Patched (2011-08-14 r56742)
## Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] googleVis_0.2.10 RJSONIO_0.95-0 
## 
## loaded via a namespace (and not attached):
## [1] tools_2.13.1

Created by Pretty R at inside-R.org

Accessing and plotting World Bank data with R

2011-09-25T15:41:00.006+01:00

Over the past couple of days I played around with the data sets of the World Bank, and I have to admit that I am blown away by it. It is amazing, to see what is available on their web site. It is worth visiting their Data Visualisation Tools page. It is fantastic that they provide an API to their data. They have used it to build an iPhone App which is pretty cool. You can have the world's data in your pocket.

In this post I will show you how we can access data from the World Bank in R. As an example we create a motion chart, in the Hans Rosling style, as you find it on the Google Public Data Explorer site, which also uses data from the World Bank. Doing this, should give us the confidence that we understand the World Bank's interface. You can find this example as demo WorldBank as part of the googleVis package from version 0.2.10 onwards.

So let's try to replicate the initial plot of the Google Public Data Explorer, which shows fertility rate against life expectancy for each country from 1960 to today, whereby the countries are represented as bubbles, with the size reflecting the population and the colour the region.

Duncan Temple Lang provides us with examples for accessing the World Bank's data using his RJSONIO and RCurl packages. The World Bank data is available via their API either as XML or JSON. We will use JSON as it is straightforward to read the JSON data set into R and to transform it into a data frame with the fromJSON function of the RJSONIO package. In order to query the data base we have to know which indicator variable we want and what its key is. Thankfully, the World Bank provides us with a page which lists all indicator variables. Clicking on any of those reveals the indicator key in the URL. For our example we get the following mappings:

Indicator	Key
fertility rate	SP.DYN.TFRT.IN
life expectancy	SP.DYN.LE00.IN
population	SP.POP.TOTL
GDP per capita (current US$)	NY.GDP.PCAP.CD

That's about it. From Duncan we have learned how to create the URL string to query the data base, and how to transform the query result from JSON into a data frame. The rest is re-arranging the data and combining the various data sets to get the final table. We display it via a motion chart using the gvisMotionChart function of the googleVis package. You find the detailed R code below.

&amp;lt;br /&amp;gt;

## This demo shows how country level data can be accessed from
## the World Bank via their API and displayed with a Motion Chart.
## Inspired by Google's Public Data Explorer, see
## http://www.google.com/publicdata/home
##
## For the World Bank Data terms of use see:
## http://data.worldbank.org/summary-terms-of-use
##
## To run this demo an internet connection and Flash are required.
## This demo is part of the googleVis R package.
##
## Markus Gesmann, 24 September 2011
## Distributed under GPL 2 or later
 
 
getWorldBankData <- function(id='SP.POP.TOTL', date='1960:2010',
                             value="value", per.page=12000){ 
  require(RJSONIO)
  url <- paste("http://api.worldbank.org/countries/all/indicators/", id,
               "?date=", date, "&format=json&per_page=", per.page,
               sep="")
 
  wbData <- fromJSON(url)[[2]]
 
  wbData = data.frame(
    year = as.numeric(sapply(wbData, "[[", "date")), 
    value = as.numeric(sapply(wbData, function(x)
      ifelse(is.null(x[["value"]]),NA, x[["value"]]))),  
    country.name = sapply(wbData, function(x) x[["country"]]['value']),
    country.id = sapply(wbData, function(x) x[["country"]]['id'])    
    )
 
  names(wbData)[2] <- value
 
  return(wbData)
}
 
getWorldBankCountries <- function(){
  require(RJSONIO)
  wbCountries <-
    fromJSON("http://api.worldbank.org/countries?per_page=12000&format=json") 
  wbCountries <- data.frame(t(sapply(wbCountries[[2]], unlist)))
  wbCountries$longitude <- as.numeric(wbCountries$longitude)
  wbCountries$latitude <- as.numeric(wbCountries$latitude)
  levels(wbCountries$region.value) <- gsub(" \\(all income levels\\)",
                                           "", levels(wbCountries$region.value))
  return(wbCountries)
}
 
## Create a string 1960:this year, e.g. 1960:2011
years <- paste("1960:", format(Sys.Date(), "%Y"), sep="")
 
## Fertility rate
fertility.rate <- getWorldBankData(id='SP.DYN.TFRT.IN',
                                   date=years, value="fertility.rate")
 
## Life Expectancy
life.exp <- getWorldBankData(id='SP.DYN.LE00.IN',  date=years,
                             value="life.expectancy") 
 
## Population
population <- getWorldBankData(id='SP.POP.TOTL',  date=years,
                               value="population")
 
## GDP per capita (current US$)
GDP.per.capita <- getWorldBankData(id='NY.GDP.PCAP.CD',
                                   date=years,
                                   value="GDP.per.capita.Current.USD") 
 
## Merge data sets
wbData <- merge(life.exp, fertility.rate)
wbData <- merge(wbData, population)
wbData <- merge(wbData, GDP.per.capita)
 
## Get country mappings
wbCountries <- getWorldBankCountries()
 
## Add regional information
wbData <- merge(wbData, wbCountries[c("iso2Code", "region.value", 
                                      "incomeLevel.value")],
                by.x="country.id", by.y="iso2Code")
 
## Filter out the aggregates and country id column
subData <- subset(wbData, !region.value %in% "Aggregates" , select=
                  -country.id) 
 
## Create a motion chart
M <- gvisMotionChart(subData, idvar="country.name", timevar="year",
                     options=list(width=700, height=600))
 
## Display the chart in your browser
plot(M)

Created by Pretty R at inside-R.org

Addition: We could simplify the code by using the WDI package by Vincent Arel-Bundock, as Diego points out in his comment below.

R in the insurance industry

2011-09-16T20:38:00.008+01:00

Let's talk about R in the insurance industry today. David Smith's blog entry reminded me about our poster at the R user conference in Warwick in August 2011:

Using R in Insurance

We presented examples on how R can be used in the insurance industry. We had a lot of fun presenting our poster. By accident we had printed the poster with quite a bit of access white space to the right. So we asked everyone who came along to sign it and by the end of the evening we had over 100 signatures!

For the historians under the readers, here is my five year old poster from GIRO in Vienna 2006.

Poster session at useR! 2011 in Warwick, UK

Yesterday Wayne Zhang, with whom I collaborate on the ChainLadder package, released the first version of his new cplm package on CRAN. The name cplm is short for compound Poisson linear models. The cplm package is for fitting Tweedie compound Poisson linear models using the Monte Carlo EM algorithm. The form of the models that are handled in the package are generalized linear models, mixed-effect models and Bayesian models. For non-Bayesian models, maximum likelihood estimations are obtained for all parameters in the model, especially for the index parameter. Estimation for the Bayesian model is performed by Markov Chain Monte Carlo simulations. These models find their application in actuarial science, see also his paper.

Here are a few more insurance related packages:

ChainLadder - Reserving methods in R. The package provides Mack-, Munich-, Bootstrap, and Multivariate-chain-ladder methods, as well as the LDF Curve Fitting methods of Dave Clark and GLM-based reserving models.
cplm - Monte Carlo EM algorithms and Bayesian methods for fitting Tweedie compound Poisson linear models.
lossDev - A Bayesian time series loss development model. Features include skewed-t distribution with time-varying scale parameter, Reversible Jump MCMC for determining the functional form of the consumption path, and a structural break in this path.
actuar: Loss distributions modelling, risk theory (including ruin theory), simulation of compound hierarchical models and credibility theory.
fitdistrplus: Help to fit of a parametric distribution to non-censored or censored data
favir: Formatted Actuarial Vignettes in R. FAViR lowers the learning curve of the R environment. It is a series of peer-reviewed Sweave papers that use a consistent style.
mondate: R packackge to keep track of dates in terms of months
lifecontingencies - Package to perform actuarial evaluation of life contingencies

Other useful documents:

Introduction to R for Actuaries by Nigel de Silva
An Interactive Introduction To R by Michael Driscoll and Dan Murphy
If you have trouble with your IT department to get R on your machine this document might help you to put some good arguments forward. The report "An Actuarial Toolkit" was presented at the GIRO convention 2006 in Vienna.

Help! There is a special interest group for R in insurance:

R-SIG-insurance -- Special Interest Group on using R in actuarial science and insurance

LondonR, 7 September 2011

2011-09-11T11:29:00.002+01:00

On 7 September 2011 I attended the London R user group meeting. It was a very good turn out with about 50 attendees at the Shooting Star, a pub close to Liverpool Street Station. The session started at 18:00 with four presentations, followed by drinks sponsored by Mango Solutions. The slides of the presentation are available on londonr.org.

The first presentation was given by Lisa Wainer from UCL Department of Security and Crime Science about crime data analysis using R. Lisa presented about a project with Merseyside police, where she had built software, in R with the gWidgets package, called the Hot Products Early Warning System, that is used to help understand and characterise the acquisitive crime problem in Merseyside on an ongoing basis, detecting emerging trends in hot products.

Chris Wood gave an insightful talk about his research on sediment biogeochemical modelling in the North Sea. His model uses a set differential equations with over 20 parameters. Chris is able to analyse and fit his model to data he gathered on an expedition in the North Sea using R, the deSolve package and having access to the super-computer at the University of Southampton. How cool is this?

Jean-Robert Avettand-Fenoel talked about the Rook package and how R and Rook has helped him to roll out new applications to his colleagues faster than using Excel, VBA and C++ or RExcel. Rook allows you to build web apps with R. The package is maintained by Jeffery Horner, who also brought us the brew package. The brew allows us, in combination with Rapache, to mix html and R code in the same file. This is quite similar to the approach taken by Sweave for LaTeX and R. However, Rook provides a way to run R web applications on your desktop with the new internal R web server named Rhttpd.

The final presentation was actually given by myself talking about the googleVis package and the recent developments in version 0.2.9:

Including googleVis output into a blogger post

2011-09-11T10:39:00.005+01:00

It seems that you cannot include Google Visualisation Charts into a blog post directly.
So, I tried to include the output of a googleVis function as a gadget, but also unsuccessfully.
Although you can include gadgets into your site template, it doesn't seem to work with blog posts. So, here is the trick which works for me: the iframe tag.
The following geo map is included as

<iframe width="100%" height="400px" frameborder="0" src="http://dl.dropbox.com/u/7586336/blogger/AndrewGeoMap.html">
</iframe>

&amp;amp;amp;lt;p&amp;amp;amp;gt;&amp;amp;amp;amp;amp;lt;br /&amp;amp;amp;amp;amp;gt; &amp;amp;amp;lt;/p&amp;amp;amp;gt;
As you can see, the chart itself is actually displayed in a page hosted by Dropbox and only inserted into this post via the iframe-tag.

For those of you, who would like to replicate the plot of Hurricane Andrew, here is the R code:

library(googleVis)
AndrewGeoMap <- gvisGeoMap(Andrew, locationvar='LatLong', numvar='Speed_kt',
                          hovervar='Category',
                          options=list(width=600,height=300,
                          region='US', dataMode='Markers'))
plot(AndrewGeoMap)
print(AndrewGeoMap, file="~/Dropbox/Public/AndrewGeoMap.html")

Created by Pretty R at inside-R.org

Correction (18 October 2011)

I just figured out that we can actually embed a chart into a blogger post directly. You can literately copy and paste the code directly into the post. However, it doesn't seem to be displayed with MS Internet Explorer.

Anyhow, here is the example from above again:


print(AndrewGeoMap, "chart", file="~/Desktop/AndrewGeoMap.js")

Now I copied and pasted the content of that file below:

googleVis 0.2.9

2011-09-04T19:43:00.003+01:00

Today we published googleVis 0.2.9 on CRAN. The new version updates the package for the new features of the Google Visualisation API and brings a new in-page editor option.

Here is a simple example, displaying the participants of the R user Conference 2011 in Warwick by country. Notice the 'Edit me' button in the top left corner of the chart, which allows you to change and customise the graph.

library(XML)
url <- "http://www.warwick.ac.uk/statsdept/useR-2011/participant-list.html"
participants <- readHTMLTable(readLines(url), which=1, stringsAsFactors=F)
names(participants) <- c("Name", "Country", "Organisation")
## Correct typo and shortcut
participants$Country <- gsub("Kngdom","Kingdom",participants$Country)
participants$Country <- gsub("USA","United States",participants$Country)
participants$Country <- factor(participants$Country)
partCountry <- as.data.frame(xtabs( ~ Country, data=participants))
library(googleVis)
## Please note the option gvis.editor requires googleVis version >= 0.2.9
G <- gvisGeoChart(partCountry,"Country", "Freq", options=list(gvis.editor="Edit me") )
plot(G)

Created by Pretty R at inside-R.org