High resolution graphics with R

For most purposes PDF or other vector graphic formats such as windows metafile and SVG work just fine. However, if I plot lots of points, say 100k, then those files can get quite large and bitmap formats like PNG can be the better option. I just have to be mindful of the resolution.

As an example I create the following plot:

x <- rnorm(100000)
plot(x, main="100,000 points", col=adjustcolor("black", alpha=0.2))

Saving the plot as a PDF creates a 5.2 MB big file on my computer, while the PNG output is only 62 KB instead. Of course, the PNG doesn't look as crisp as the PDF file.

png("100kPoints72dpi.png", units = "px", width=400, height=400)
plot(x, main="100,000 points", col=adjustcolor("black", alpha=0.2))
dev.off()

Hence, I increase the resolution to 150 dots per pixel.

png("100kHighRes150dpi.png", units="px", width=400, height=400, res=150)
plot(x, main="100,000 points", col=adjustcolor("black", alpha=0.2))
dev.off()

This looks a bit odd. The file size is only 29 KB but the annotations look too big. Well, the file has only 400 x 400 pixels and the size of a pixel is fixed. Thus, I have to provide more pixels, or in other words increase the plot size. Doubling the width and height as I double the resolution makes sense.

png("100kHighRes150dpi2.png", units="px", width=800, height=800, res=150)
plot(x, main="100,000 points", col=adjustcolor("black", alpha=0.2))
dev.off()

Next I increase the resolution further to 300 dpi and the graphic size to 1600 x 1600 pixels. The file is still very crisp. Of course the file size increased. Now it is 654 KB in size, yet sill only about 1/8 of the PDF and I can embed it in LaTeX as well.

png("100kHighRes300dpi.png", units="px", width=1600, height=1600, res=300)
plot(x, main="100,000 points", col=adjustcolor("black", alpha=0.2))
dev.off()

Note, you can click on the charts to access the original files of this post.

Setting axes limits with googleVis

I posted about the various googleVis axis options for base charts, such as line, bar and area charts earlier, but I somehow forgot to mention how to set the axes limits.

Unfortunately, there are no arguments such as ylim and xlim. Instead, the Google Charts axes options are set via hAxes and vAxes, with h and v indicating the horizontal and vertical axis. More precisely, I have to set viewWindowMode : 'explicit' and set the viewWindow to the desired min and max values. Additionally, I have to wrap all of this in [{}] brackets as those settings are sub-options of h/vAxes. There are also options minValue and maxValue, but they only allow you to extend the axes ranges.

Here is a minimal example, setting the y-axis limits from 0 to 10:




With more than one variable to plot I can use the series argument to decide which variables I want on the left and right axes and set the viewWindow accordingly. Again, here is a minimal example:

Session Info

sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/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.4.3

loaded via a namespace (and not attached):
[1] RJSONIO_1.0-3

googleVis tutorial at useR!2013

Today Diego and I will give our googleVis tutorial at useR!2013 in Albacete, Spain.

googleVis Tutorial at useR! 2013

We will cover:

• Introduction and motivation
• Google Chart Tools
• R package googleVis
  
  • Concepts of googleVis
  • Case studies
• googleVis on shiny

R package development

Building R packages is not particular hard, but it can be a bit of a daunting endeavour at the beginning, particularly if you are more of a statistician than a computer scientist or programmer.

Some concepts may appear foreign or like red tape, yet many of them evolved over time for a reason. They help to stay organise, collaborate more effectively with others and write better code.

So, here are my slides of the R package development workshop at Lancaster University.

R package development

For a detailed and authoritative reference on R package development see the Writing R Extensions manual on CRAN.

Interactive slides with googleVis on shiny

Following on from last week's post, here are my slides on using googleVis on shiny from the Advanced R workshop at Lancaster University, 21 May 2013.

googleVis on shiny

Again, I wrote my slides in RMarkdown and I used slidify to create the HTML5 presentation. Unfortunately you may have to reload the slides that use googleVis on shiny as the JavaScript code in the background is potentially not ideal. Any pointers, which could help to improve the performance will be much appreciated.

Many of the examples in my slides are taken from my post First steps of using googleVis on shiny, however the presentation also demonstrates that it is possible to inject JavaScript code into a googleVis chart to trigger a shiny event, see also the example below.
Read more »

Interactive presentation with slidify and googleVis

Last week I was invited to give an introduction to googleVis at Lancaster University. This time I decided to use the R package slidify for my talk. Slidify, like knitr, is built on Markdown and makes it very easy to create beautiful HTML5 presentations.

Introduction to googleVis

Separating content from layout is always a good idea. Markup languages such as TeX/LaTeX or HTML are built on this principle. Ramnath Vaidyanathan has done a fantastic job with slidify, as it is very straightforward to create presentations with R. There are a couple of advantages compared to traditional presentation software packages:

• RMarkdown helps me to focus on the content
• Integration of R code is build in
• HTML5 allows me to embed interactive content, such as
  
  • Videos
  • googleVis and other interactive charts
  • shiny apps (more on this next week)

In the past I have used knitr in combination with pandoc to generate a slidy presentation. However, with slidfiy I can do all this in R directly. And better, Ramnath provides me with a choice of different layout frameworks and syntax highlighting options. Finally to top it all, publishing the slides on Github was only one more R statement: publish('mages', 'Introduction_to_googleVis').

I will give a half-day tutorial on googleVis with Diego de Castillo at useR2013! in Albacete on 9 July 2013. I hope to see some of you there.

How to change the alpha value of colours in R

Often I like to reduce the alpha value (level of transparency) of colours to identify patterns of over-plotting when displaying lots of data points with R. So, here is a tiny function that allows me to add an alpha value to a given vector of colours, e.g. a RColorBrewer palette, using col2rgb and rgb, which has an argument for alpha, in combination with the wonderful apply and sapply functions.


The example below illustrates how this function can be used with colours provided in different formats, thanks to the col2rgb function.

Read more »

How to set axis options in googleVis

Setting axis options in googleVis charts can be a bit tricky. Here I present two examples where I set several options to customise the layout of a line and combo chart with two axes.

The parameters have to be set in line with the Google Chart Tools API, which uses a JavaScript syntax. In googleVis chart options are set via a list in the options argument. Some of the list items can be a bit more complex, often wrapped in {} brackets, e.g. for various formatting options or in [] brackets, if there are multiple series to consider. Within those brackets sub-options are set via argument : value, using the : character for assignments.

There are many other options as part of the Google Chart Tools API, which are not supported by googleVis yet, such as columns roles, controls and dashboards, etc. Please get in touch if you have ideas in this regard and/or would like to collaborate.

In my first example I display two series of dummy data in a line chart with two axes. The left hand scale is in percentages and the right hand scale in amounts. Note in the code below how I set the various parameters and the placements of the different kinds of brackets.
Read more »

How to use optim in R

A friend of mine asked me the other day how she could use the function optim in R to fit data. Of course there are functions for fitting data in R and I wrote about this earlier. However, she wanted to understand how to do this from scratch using optim.

The function optim provides algorithms for general purpose optimisations and the documentation is perfectly reasonable, but I remember that it took me a little while to get my head around how to pass data and parameters to optim. Thus, here are two simple examples.

I start with a linear regression by minimising the residual sum of square and discuss how to carry out a maximum likelihood estimation in the second example.

Minimise residual sum of squares

I start with an x-y data set, which I believe has a linear relationship and therefore I'd like to fit y against x by minimising the residual sum of squares.

dat=data.frame(x=c(1,2,3,4,5,6), 
               y=c(1,3,5,6,8,12))

Next, I create a function that calculates the residual sum of square of my data against a linear model with two parameter. Think of y = par[1] + par[2] * x.

min.RSS <- function(data, par) {
              with(data, sum((par[1] + par[2] * x - y)^2))
             }

Optim minimises a function by varying its parameters. The first argument of optim are the parameters I'd like to vary, par in this case; the second argument is the function to be minimised, min.RSS. The tricky bit is to understand how to apply optim to your data. The solution is the ... argument in optim, which allows me to pass other arguments through to min.RSS, here my data. Therefore I can use the following statement:

result <- optim(par = c(0, 1), min.RSS, data = dat)
# I find the optimised parameters in result$par
# the minimised RSS is stored in result$value
result
## $par
## [1] -1.267  2.029
## 
## $value
## [1] 2.819
## 
## $counts
## function gradient 
##       89       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Let me plot the result:

plot(y ~ x, data = dat)
abline(a = result$par[1], b = result$par[2], col = "red")

Read more »

First steps of using googleVis on shiny

The guys at RStudio have done a fantastic job with shiny. It is really easy to build web apps with R using shiny. With the help of Joe Cheng from RStudio we figured out a way to make googleVis work on shiny as well. This allows you to make use of the Google Charts Tools in your shiny app directly from R. What I present here are three initial examples which seem to work in most browsers. The third example even uses a neat trick to create an animated geo chart.



However, before we upload the next version of googleVis to CRAN we decided to present a preview of version 0.4.0 here, asking for feedback. It would not be fair on the guys behind CRAN to release something into the wild, only to be told by users within a few days that we missed something. Hence, you can get the new version of googleVis only from the download page of our project site for the time being.

You may have read the post on RStudio's blog that the shiny API changed slightly: the reactivePlot and reactiveText functions have been renamed to renderPlot and renderText with simplified input parameters. Thanks to Joe, there is now also a renderGvis function as part of the googleVis package, which works in very much the same way as the other two.

To run the following examples you need shiny version 0.4.0 and googleVis version 0.4.0 or higher. Note that the charts below are not screen shots, but real apps hosted on RStudio's infrastructure.
Read more »

Reserving based on log-incremental payments in R, part III

This is the third post about Christofides' paper on Regression models based on log-incremental payments [1]. The first post covered the fundamentals of Christofides' reserving model in sections A - F, the second focused on a more realistic example and model reduction of sections G - K. Today's post will wrap up the paper with sections L - M and discuss data normalisation and claims inflation.

I will use the same triangle of incremental claims data as introduced in my previous post. The final model had three parameters for origin periods and two parameters for development periods. It is possible to reduce the model further as Christofides illustrates in section L onwards by using an inflation index to bring all claims payments to current value and a claims volume adjustment or weight for each origin period to normalise the triangle.

In his example Christofides uses claims volume adjustments for the origin years and an earning or inflation index for the different payment calendar years. The claims volume adjustments aims to normalise the triangle for similar exposures across origin periods, while the earnings index, which measures largely wages and other forms of compensations, is used as a first proxy for claims inflation. Note that the earnings index shows significant year on year changes from 5% to 9%. Barnett and Zehnwirth [2] would probably recommend to add further parameters for the calendar year effects to the model.

# Page D5.36
ClaimsVolume <- data.frame(origin=0:6, 
  volume.index=c(1.43, 1.45, 1.52, 1.35, 1.29, 1.47, 1.91))
# Page D5.36
EarningIndex <- data.frame(cal=0:6, 
  earning.index=c(1.55, 1.41, 1.3, 1.23, 1.13, 1.05, 1))
# Year on year changes
round((1-EarningIndex$earning.index[-1]/EarningIndex$earning.index[-7]),2)
# [1] 0.09 0.08 0.05 0.08 0.07 0.05

dat <- merge(merge(dat, ClaimsVolume), EarningIndex)

# Normalise data for volume and earnings
dat$logvalue.ind.inf <- with(dat, log(value/volume.index*earning.index))
 
with(dat, interaction.plot(dev, origin, logvalue.ind.inf))
points(1+dat$dev, dat$logvalue.ind.inf, pch=16, cex=0.8)
Indeed, the interaction plot shows the various origin years now to be much more closely grouped. Only the single point of the last origin period stands out now.  Christofides tests several models with different numbers of origin levels, but I am happy with the minimal model using only one parameter for the origin period, namely the intercept:

Read more »

Reserving based on log-incremental payments in R, part II

Following on from last week's post I will continue to go through the paper Regression models based on log-incremental payments by Stavros Christofides [1]. In the previous post I introduced the model from the first 15 pages up to section F. Today I will progress with sections G to K which illustrate the model with a more realistic incremental claims payments triangle from a UK Motor Non-Comprehensive account:

# Page D5.17
tri <- t(matrix(
  c(3511, 3215, 2266, 1712, 1059, 587, 340,
    4001, 3702, 2278, 1180,  956, 629,  NA,
    4355, 3932, 1946, 1522, 1238,  NA,  NA,
    4295, 3455, 2023, 1320,   NA,  NA,  NA,
    4150, 3747, 2320,   NA,   NA,  NA,  NA,
    5102, 4548,   NA,   NA,   NA,  NA,  NA,
    6283,   NA,   NA,   NA,   NA,  NA,  NA), nc=7))

The rows show origin period data, e.g. accident years, underwriting years or years of account and the columns present the development periods or lags. The triangle appears to be fairly well behaved. The last two years in rows 6 and 7 appear to be slightly higher than rows 2 to 5 and the values in row 1 are lower in comparison to the later years. The last payment of £1,238 in the third row stands out a bit as well.

Before I plot the data, I will transform the triangle into a data frame and add extra columns:

m <- dim(tri)[1]; n <- dim(tri)[2]
dat <- data.frame(
  origin=rep(0:(m-1), n),
  dev=rep(0:(n-1), each=m),
  value=as.vector(tri))

## Add dimensions as factors
dat <- with(dat, data.frame(origin, dev, cal=origin+dev, 
                            value, logvalue=log(value),
                            originf=factor(origin),
                            devf=as.factor(dev),
                            calf=as.factor(origin+dev)))

I am particularly interested in the decay of claims payments in the development year direction for each origin year on the original and log-scale. The interaction.plot of the stats package does an excellent job for this:

op <- par(mfrow=c(2,1), mar=c(4,4,2,2))
with(dat, interaction.plot(x.factor=dev, trace.factor=origin, 
          response=value))
points(dat$devf, dat$value, pch=16, cex=0.5)
with(dat, interaction.plot(x.factor=dev, trace.factor=origin, 
          response=logvalue))
points(dat$devf, dat$logvalue, pch=16, cex=0.5)
par(op)

Indeed the origin years 1 to 4 (rows 2 to 5) look quite similar and the decay of claims in development year direction appears to be linear on a log-scale from development year 1 onwards.

Based on those observations Christofides suggests two models; the first one will have a unique level for each origin year and a unique level for the zero development period. The parameters for development periods 1 to 6 are assumed to follow a linear relationship with the same slope \(s\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + d_j + \epsilon_{ij}
&\mbox{for } i,\,j \mbox{ from } 0 \mbox{ to } 6\\
\mbox{where } d_0 &= d,\quad d_j = s \cdot j
&\mbox{for } j > 0
\end{align}and \(\epsilon_{ij} \sim N(0, \sigma^2)\). The second model will be a reduced version of the above with only two levels for the origin years 5 and 6. Hence, I add four more columns to my data frame:Read more »

Reserving based on log-incremental payments in R, part I

A recent post on the PirateGrunt blog on claims reserving inspired me to look into the paper Regression models based on log-incremental payments by Stavros Christofides [1], published as part of the Claims Reserving Manual (Version 2) of the Institute of Actuaries.

The paper is available together with a spread sheet model, illustrating the calculations. It is very much based on ideas by Barnett and Zehnwirth, see [2] for a reference. However, doing statistical analysis in a spread sheet programme is often cumbersome. I will go through the first 15 pages of Christofides' paper today and illustrate how the model can be implemented in R.

Let's start with the example data of an incremental claims triangle:

## Page D5.4
tri <- t(matrix(
  c(11073,  6427, 1839, 766,
    14799,  9357, 2344,  NA,
    15636, 10523,   NA,  NA,
    16913,    NA,   NA,  NA), 
  nc=4, dimnames=list(origin=0:3, dev=0:3)))

The above triangle shows incremental claims payments for four origin (accident) years over time (development years). It is the aim to predict the bottom right triangle of future claims payments, assuming no further claims after four development years.

Christofides model assumes the following structure for the incremental paid claims \(P_{ij}\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + b_j + \epsilon_{ij}
\end{align}where i and j go from 0 to 3, \(b_0=0\) and \(\epsilon_{ij} \sim N(0, \sigma^2)\). Unlike the basic chain-ladder method, this is a stochastic model that allows me to test my assumptions and calculate various statistics, e.g. standards errors of my predictions.
Read more »

Changing colours and legends in lattice plots

Lattice plots are a great way of displaying multivariate data in R. Deepayan Sarkar, the author of lattice, has written a fantastic book about Multivariate Data Visualization with R [1]. However, I often have to refer back to the help pages to remind myself how to set and change the legend and how to ensure that the legend will use the same colours as my plot. Thus, I thought I write down an example for future reference.

Eventually I would like to get the following result:

I start with a simple bar chart, using the Insurance data set of the MASS package.

Please note that if you use R-2.15.2 or higher you may get a warning message with the current version of lattice (0.20-10), but this is a known issue, and I don't think you need to worry about it.

library(MASS) 
library(lattice)
## Plot the claims frequency against age group by engine size and district
barchart(Claims/Holders ~ Age | Group, groups=District,
         data=Insurance, origin=0, auto.key=TRUE)

To change the legend I specify a list with the various options for auto.key. Here I want to set the legend items next to each other, add a title for the legend and change its font size: Read more »

Data.table rocks! Data manipulation the fast way in R

I really should make it a habit of using data.table. The speed and simplicity of this R package are astonishing.

Here is a simple example: I have a data frame showing incremental claims development by line of business and origin year. Now I would like add a column with the cumulative claims position for each line of business and each origin year along the development years.

It's one line with data.table! Here it is:

myData[order(dev), cvalue:=cumsum(value), by=list(origin, lob)]

It is even easy to read! Notice also that I don't have to copy the data. The operator ':=' works by reference and is one of the reasons why data.table is so fast.

And it is getting even better. Suppose you want to get the latest claims development position for each line of business and origin year. Again, it is only one line: Read more »

Simulating neurons or how to solve delay differential equations in R

I discussed earlier how the action potential of a neuron can be modelled via the Hodgkin-Huxely equations. Here I will present a simple model that describes how action potentials can be generated and propagated across neurons. The tricky bit here is that I use delay differential equations (DDE) to take into account the propagation time of the signal across the network.

My model is based on the paper: Epileptiform activity in a neocortical network: a mathematical model by F. Giannakopoulos, U. Bihler, C. Hauptmann and H. J. Luhmann. The article presents a flexible and efficient modelling framework for:

• large populations with arbitrary geometry
• different synaptic connections with individual dynamic characteristics
• cell specific axonal dynamics

Read more »

Connecting the real world to R with an Arduino

If connecting data to the real world is the next sexy job, then how do I do this? And how do I connect the real world to R?

It can be done as Matt Shottwell showed with his home made ECG and a patched version of R at useR! 2011. However, there are other options as well and here I will use an Arduino. The Arduino is an open-source electronics prototyping platform. It has been around for a number of years and is very popular with hardware hackers. So, I had to have a go at the Arduino as well.

My Arduino starter kit from oomlout

The example I will present here is silly - it doesn't do anything meaningful and yet I believe it shows the core building blocks for future projects: Read an analog signal into the computer via the Arduino, transform it with R through Rserve and display it graphically in real time. The video below demonstrates the final result. As I turn the potentiometer random points are displayed on the screen, with the standard deviation set by the analog output (A0) of the Arduino and fed into the rnorm function in R, while at the same time the LED brightness changes.

Read more »

Using R in Insurance, Presentation at GIRO 2012

Every year the UK’s general insurance actuarial community organises a big conference, which they call GIRO, short for General Insurance Research Organising committee.

This year's conference is in Brussels from 18 - 21 September 2012. Despite the fact that Brussels is actually in Belgium the UK actuaries will travel all the way to enjoy good beer and great talks.

On Wednesday morning I will run a session on Using R in insurance. It would be great to see some of you there.

I prepared the slides with R, RStudio, knitr, pandoc and slidy again. My title page shows a word cloud about the GIRO conference. It uses the wordcloud package and was inspired by Ian Fellows' post on FellStats.

The last slide shows the output of sessionInfo(). I am sure it will become helpful one day, when I have to remind myself how I actually created the slides and which packages and versions I used.

Interactive web graphs with R - Overview and googleVis tutorial

Today I feel very lucky, as I have been invited to the Royal Statistical Society conference to give a tutorial on interactive web graphs with R and googleVis.

I prepared my slides with RStudio, knitr, pandoc and slidy, similar to my Cambridge R talk. You can access the RSS slides online here and you find the original R-Markdown file on github. You will notice some HTML code in the file, which I had to use to overcome my knowledge gaps of Markdown or its limitations. However, my output format will always be HTML, so that should be ok. To convert the Rmd-file into a HTML slidy presentation execute the following statements on the command line:

Rscript -e "library(knitr); knit('googleVis_at_RSS_2012.Rmd')"
pandoc -s -S -i -t slidy --mathjax googleVis_at_RSS_2012.md 
  -o googleVis_at_RSS_2012.html

Read more »

Rook rocks! Example with googleVis

What is Rook?

Rook is a web server interface for R, written by Jeffrey Horner, the author of rApache and brew. But unlike other web frameworks for R, such as brew, R.rsp (which I have used in the past¹), Rserve, gWidgetWWWW or sumo (which I haven't used yet) Rook appears incredible lightweight.

Rook doesn't need any configuration. It is an R package, which works out of the box with the R HTTP server (R ≥ 2.13.0 required). That means no configuration files are needed. No files have to be placed in particular folders, and I don't have to worry about access permissions. Instead, I can run web applications on my local desktop.

Screen shot of a Rook app running in a browser

Web applications have the big advantage that they run in a browser and hence are somewhat independent of the operating systems and installed software. All I need is R and a web browser. But here is the catch: I have to learn a little bit about the HTTP protocol to develop Rook apps.

Read more »