CYBAEA Data and Analysis

R code for Chapter 2 of Non-Life Insurance Pricing with GLM

2012-03-13T16:57:00Z

Amazon UK | US

We continue working our way through the examples, case studies, and exercises of what is affectionately known here as “the two bears book” (Swedish björn = bear) and more formally as Non-Life Insurance Pricing with Generalized Linear Models by Esbjörn Ohlsson and Börn Johansson (Amazon UK | US).

At this stage, our purpose is to reproduce the analysis from the book using the R statistical computing and analysis platform, and to answer the data analysis elements of the exercises and case studies. Any critique of the approach and of pricing and modeling in the Insurance industry in general will wait for a later article.

In the following, we will assume that the reader has a copy of the book and a working installation of R. We will be using the data.table, foreach, and ggplot2 packages which are not part of the standard distribution, so the reader should install them first (e.g. by executing install.packages(c("data.table", "foreach", "ggplot2"), dependencies = TRUE) from within an R session).

The beginning of the code is not so important, but here goes:

#!/usr/bin/Rscript
## PricingGLM-2.r - Code for Chapter 2 of "Non-Life Insurance Pricing with GLM"
## Copyright © 2012 CYBAEA Limited (http://www.cybaea.net)

## @book{ohlsson2010non,
##   title={Non-Life Insurance Pricing with Generalized Linear Models},
##   author={Ohlsson, E. and Johansson, B.},
##   isbn={9783642107900},
##   series={Eaa Series: Textbook},
##   url={http://books.google.com/books?id=l4rjeflJ\_bIC},
##   year={2010},
##   publisher={Springer Verlag}
## }

Example 2.5: Moped insurance continued

We continue the moped insurance example, and we use the data that we saved in our Chapter 1 session. The goal is to reproduce Table 2.7 so we start building that as a data frame after loading the data.

## Load the data from last.
if (!exists("table.1.2"))
    load("table.1.2.RData")

library("foreach")

## We are looking to reproduce table 2.7 which we start building here,
## adding columns as we go.
table.2.7 <-
    data.frame(rating.factor =
               c(rep("Vehicle class", nlevels(table.1.2$premiekl)),
                 rep("Vehicle age",   nlevels(table.1.2$moptva)),
                 rep("Zone",          nlevels(table.1.2$zon))),
               class =
               c(levels(table.1.2$premiekl),
                 levels(table.1.2$moptva),
                 levels(table.1.2$zon)),
               stringsAsFactors = FALSE)

We next calculate the duration and number of claims for each level of each rating factor. We also set the contrasts for the levels, using the same idiom as in our Chapter 1 session. The foreach package is convenient to use here, but you can of course do it with a normal loop and a couple of variables.

## Calculate duration per rating factor level and also set the
## contrasts (using the same idiom as in the code for the previous
## chapter). We use foreach here to execute the loop both for its
## side-effect (setting the contrasts) and to accumulate the sums.
new.cols <-
    foreach (rating.factor = c("premiekl", "moptva", "zon"),
             .combine = rbind) %do%
{
    nclaims <- tapply(table.1.2$antskad, table.1.2[[rating.factor]], sum)
    sums <- tapply(table.1.2$dur, table.1.2[[rating.factor]], sum)
    n.levels <- nlevels(table.1.2[[rating.factor]])
    contrasts(table.1.2[[rating.factor]]) <-
        contr.treatment(n.levels)[rank(-sums, ties.method = "first"), ]
    data.frame(duration = sums, n.claims = nclaims)
}
table.2.7 <- cbind(table.2.7, new.cols)
rm(new.cols)

Next, we need to build the frequency and severity models separately, as in the discussion at the beginning of section 2.3.4 on page 34.

Frequency model

Amazon UK | US

Note here the use of the offset() term as opposed to the weights= argument we have used before. The offset is log(dur) because our link function is log.

See help("Insurance", package = "MASS") for a similar example and sections 7.1 (p. 189--190) and 7.3 of Modern Applied Statistics with S (Amazon UK | US ) for a (very brief) discussion.

model.frequency <-
    glm(antskad ~ premiekl + moptva + zon + offset(log(dur)),
        data = table.1.2, family = poisson)

rels <- coef( model.frequency )
rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] )
table.2.7$rels.frequency <-
    c(c(1, rels[1])[rank(-table.2.7$duration[1:2], ties.method = "first")],
      c(1, rels[2])[rank(-table.2.7$duration[3:4], ties.method = "first")],
      c(1, rels[3:8])[rank(-table.2.7$duration[5:11], ties.method = "first")])

Severity model

There are a couple of points to note here:

We will model using a Gamma distribution for the errors, following the discussion on page 20 and also pages 33-34. Note the point that this is only one of several plausible candidate distributions.
Because we are using the Gamma distribution we need to remove the zero values from the data; we do this using the table.1.2[table.1.2$medskad > 0, ] construct.
To reproduce the values from the book, we use the non-canonical "log" link function even though the canonical function ("inverse") gives a slightly better fit (residual deviance 5.9 versus 8.0 on 16 degrees of freedom). This follows the approach discussed in Example 2.3 on page 30.

model.severity <-
    glm(medskad ~ premiekl + moptva + zon,
        data = table.1.2[table.1.2$medskad > 0, ],
        family = Gamma("log"), weights = antskad)

rels <- coef( model.severity )
rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] )
## Aside: For the canonical link function use
## rels <- rels[1] / (rels[1] + rels[-1])

table.2.7$rels.severity <-
    c(c(1, rels[1])[rank(-table.2.7$duration[1:2], ties.method = "first")],
      c(1, rels[2])[rank(-table.2.7$duration[3:4], ties.method = "first")],
      c(1, rels[3:8])[rank(-table.2.7$duration[5:11], ties.method = "first")])

Combining the models

Now it is trivial to combine and display the results.

table.2.7$rels.pure.premium <- with(table.2.7, rels.frequency * rels.severity)
print(table.2.7, digits = 2)

	rating.factor	class	duration	n.claims	rels.frequency	rels.severity	rels.pure.premium
1	Vehicle class	1	9833.20	391	1.00	1.00	1.00
2	Vehicle class	2	8825.10	395	0.78	0.55	0.42
11	Vehicle age	1	1918.40	141	1.55	1.79	2.78
21	Vehicle age	2	16739.90	645	1.00	1.00	1.00
12	Zone	1	1451.40	206	7.10	1.21	8.62
22	Zone	2	2486.30	209	4.17	1.07	4.48
3	Zone	3	2888.70	132	2.23	1.07	2.38
4	Zone	4	10069.10	207	1.00	1.00	1.00
5	Zone	5	246.10	6	1.20	1.21	1.46
6	Zone	6	1369.20	23	0.79	0.98	0.78
7	Zone	7	147.50	3	1.00	1.20	1.20

2.4 Case Study: Motorcycle Insurance

The authors use the term Case Study for larger exercises. There is quite a lot to cover here, and some of the questions are more about discussing the results. We will focus on getting the results.

The Data

Obtaining the data

First we read the (fixed width) data. The column names are the same (based on Swedish language) as in the book.

columns <- c(agarald = 2L, kon = 1L, zon = 1L, mcklass = 1L, fordald = 2L,
             bonuskl = 1L, duration = 8L, antskad = 4L, skadkost = 8L)
column.classes <- c("integer", rep("factor", 3), "integer",
                    "factor", "numeric", rep("integer", 2))
stopifnot(length(columns) == length(column.classes))
con <- url("http://www2.math.su.se/~esbj/GLMbook/mccase.txt")
mccase <- read.fwf(con, widths = columns, header = FALSE,
                   col.names = names(columns),
                   colClasses = column.classes,
                   na.strings = NULL, comment.char = "")
try(close(con), silent = TRUE)
rm(columns, column.classes, con)
mccase$mcklass <- ordered(mccase$mcklass)
mccase$bonuskl <- ordered(mccase$bonuskl)

Adding meta-data information

We are more than a little obsessed with documenting our data, so here goes. See the book for the details.

comment(mccase) <-
    c("Title: Partial casco insurance for motorcycles from Wasa, 1994--1998",
      "Source: http://www2.math.su.se/~esbj/GLMbook/mccase.txt",
      "Copyright: http://www2.math.su.se/~esbj/GLMbook/")
comment(mccase$agarald) <-
    c("The owner's age, between 0 and 99",
      "Name: Age of Owner")
comment(mccase$kon) <-
    c("Name: Gender of Owner",
      "Code: M=Male",
      "Code: K=Female")
comment(mccase$zon) <-
    c("Name: Geographic Zone",
      "Code: 1=Central and semi-central parts of Sweden's three largest cities",
      "Code: 2=suburbs and middle-sized towns",
      "Code: 3=Lesser towns, except those in 5 or 7",
      "Code: 4=Small towns and countryside, except 5--7",
      "Code: 5=Northern towns",
      "Code: 6=Northern countryside",
      "Code: 7=Gotland (Sweden's largest island)")
comment(mccase$mcklass) <-
    c("Name: MC Class",
      "Description: A classification by the EV ratio, defined as (engine power in kW × 100) / (vehicle weight in kg + 75), rounded to the nearest lower integer.",
      "Code: 1=EV ratio <= 5",
      "Code: 2=EV ratio 6--8",
      "Code: 3=EV ratio 9--12",
      "Code: 4=EV ratio 13--15",
      "Code: 5=EV ratio 16--19",
      "Code: 6=EV ratio 20--24",
      "Code: 7=EV ratio >= 25")
comment(mccase$fordald) <-
    c("Vehicle age, between 0 and 99",
      "Name: Vehicle Age")
comment(mccase$bonuskl) <-
    c("Name: Bonus Class",
      "Description: A driver starts with bonus class 1; for each claim-free year the bonus class is increased by 1. After the first claim the bonus is decreased by 2; the driver cannot return to class 7 with less than 6 consecutive claim free years.")
comment(mccase$duration) <-
    c("Name: Duration",
      "Comment: The number of policy years",
      "Unit: year")
comment(mccase$antskad) <-
    c("Name: Number of Claims")
comment(mccase$skadkost) <-
    c("Name: Cost of Claims",
      "Unit: SEK")

Rating factors

We finally add the rating factors from the book.

mccase$rating.1 <- mccase$zon
mccase$rating.2 <- mccase$mcklass
mccase$rating.3 <-
    cut(mccase$fordald, breaks = c(0, 1, 4, 99),
        labels = as.character(1:3), include.lowest = TRUE,
        ordered_result = TRUE)
mccase$rating.4 <- ordered(mccase$bonuskl) # Drop comments
levels(mccase$rating.4) <-                 # Combine levels
    c("1", "1", "2", "2", rep("3", 3))

Save the data

Never forget to save.

save(mccase, file = "mccase.RData")

Now we can start tackling the exercises.

Problem 1: Aggregate to cells of current tariff

Aggregate the data to the cells of the current tariff. Compute the empirical claim frequency and severity at this level.

We really like the data.table package. The syntax is much easier on the eye (and the hand). As a bonus, it scales well to much larger data sets than data.frame. If you are not already using it, now is the time to start.

if (!exists("mccase"))
    load("mccase.RData")

## Conver to data.table
library("data.table")
mccase <- data.table(mccase, key = paste("rating", 1:4, sep = "."))

## Aggregate to levels of current rating factors
mccase.current <-
    mccase[,
           list(duration = sum(duration),
                antskad = sum(antskad),
                skadkost = sum(skadkost),
                num.policies = .N),
           by = key(mccase)]             

## Claim frequency and severity. Change NaN to NA.
mccase.current$claim.freq <-
    with(mccase.current, ifelse(duration != 0, antskad / duration, NA_real_))
mccase.current$severity <-
    with(mccase.current, ifelse(antskad != 0, skadkost / antskad, NA_real_))

## Save
save(mccase.current, file = "mccase.current.RData")

Problem 2: Determine how the duration and number of claims is distributed

Determine how the duration and number of claims is distributed for each of the rating factor classes, as an indication of the accuracy of the statistical analysis.

This is one of those 'there is no right answer' questions, but let us have a look at the data. First load it if needed and also load the libraries we will be using.

## Load data if needed
library("data.table")
if (!exists("mccase"))
    load("mccase.RData")
if (!exists("mccase.current"))
    load("mccase.current.RData")
if (!is(mccase, "data.table"))
    mccase <- data.table(mccase, key = paste("rating", 1:4, sep = "."))

library("grid")
library("ggplot2")

1. Number of claims (antskad)

Let us start with the number of claims for the undelying data. First we plot the number of policies for each number of claims (range of number of claims is {0, 1, 2}); this is one way to do it:

plot.titles <- c("Geo. zone", "MC class", "Vehicle age", "Bonus class")
plots <-
    lapply(1:4,
           function(i)
           ggplot(mccase, aes(antskad))
           + geom_histogram(aes(weight = duration))
           + scale_x_discrete(limits = c(0, 2))
           + scale_y_log10()
           + facet_grid(paste("rating.", i, " ~ .", sep = ""),
                        scales = "fixed")
           ## We drop the axis titles to make more room for the data
           + opts(axis.title.x = theme_blank(), axis.title.y = theme_blank(),
                  axis.text.x = theme_blank(),  axis.text.y = theme_blank(),
                  axis.title.y = theme_blank(), axis.ticks = theme_blank(),
                  title = plot.titles[i])
           )

grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow = 1, ncol = 4)))
## We can ignore the warnings from displaying the plots for now
for (i in 1:4)
    print(plots[[i]], vp = viewport(layout.pos.col = i))

There is certainly some variation there, especially on the first rating factor.

Recall from the discussion on page 18 the assumption that the number of claims for a single policy is a Poisson process and the distribution of the number of claims in a tariff cell is Poisson distributed. Let us first look at the whole data set (and note that the data.table notation is used; if using data.frame be sure to have drop=FALSE):

data <- mccase[order(antskad), list(N = .N, w = sum(duration)), by = antskad]
M <- glm(N ~ antskad, family = poisson(), weights = w, data = data)
data$predicted <- round(predict(M, data[, list(antskad)], type = "response"))
print(data[, list(antskad, N, predicted)], digits = 1)

	antskad	n	predicted
1	0	63878	63878.0
2	1	643	646.0
3	2	27	7.0

Not too bad, but then it can't really be with that heavy weighting to the first value of antskad. We can show the same split by the levels of the first rating factor as an example:

data <-
    mccase[order(rating.1, antskad),
           list(N = .N, w = sum(duration)),
           by = list(rating.1, antskad)]
modelPoisson <- function(data) {
    M <- glm(N ~ antskad, family = poisson(), weights = w, data = data)
    return(M)
}
predictionPoisson <- function (data) {
    M <- modelPoisson(data)
    p <- predict(M, data[, list(antskad)], type = "response")
    return(p)
}
data$predicted <-
    unlist(lapply(levels(data$rating.1),
                  function (l) round(predictionPoisson(data[rating.1 == l]))))
print(data[, list(rating.1, antskad, N, predicted)], digits = 1)

	rating.1	antskad	n	predicted
1	1	0	8409	8409.0
2	1	1	163	164.0
3	1	2	10	3.0
4	2	0	11632	11632.0
5	2	1	157	157.0
6	2	2	5	2.0
7	3	0	12604	12604.0
8	3	1	113	113.0
9	3	2	5	1.0
10	4	0	24626	24626.0
11	4	1	184	185.0
12	4	2	6	1.0
13	5	0	2368	2368.0
14	5	1	9	9.0
15	6	0	3867	3867.0
16	6	1	16	16.0
17	6	2	1	0.0
18	7	0	372	372.0
19	7	1	1	1.0

You get the idea....

If we can't really see very much at the individual claims level, how does it appear when we look at the levels aggregated to the current tariff cells? Here, of course, we have much more data with up to 33 claims in one cell, but we will limit the display to the first few number of claims.

high.limit <- 7L                 # Only display up to this many claims
plots <-
    lapply(1:4,
           function(i)
           ggplot(mccase.current, aes(antskad))
           + geom_histogram(breaks = 0:high.limit, aes(weight = duration))
           + scale_x_discrete(limits = c(0, high.limit), breaks = 0:high.limit)
           ## Following xlim() needed for scale_x_discrete only; see
           ## https://groups.google.com/d/msg/ggplot2/wLWGCUz8K6k/DeVudyfXyKgJ
           + xlim(0, high.limit)
           + scale_y_log10()
           + facet_grid(paste("rating.", i, " ~ .", sep = ""),
                        scales = "fixed")
           ## We drop the axis titles to make more room for the data
           + opts(axis.title.x = theme_blank(), axis.title.y = theme_blank(),
                  axis.text.x = theme_blank(),  axis.text.y = theme_blank(),
                  axis.title.y = theme_blank(), axis.ticks = theme_blank(),
                  title = plot.titles[i])
           )

grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow = 1, ncol = 4)))
for (i in 1:4)
    print(plots[[i]], vp = viewport(layout.pos.col = i))

There is certainly something going on for low bonus classes and high vehicle ages that is not straightforward Poisson distribution. Looking at the two in combination we see the old motorcycle - low bonus class (rating.3 = 3 and rating.4 = 1) is the main culpit, but also note that there is relatively little data in these cells.

ggplot(mccase.current, aes(antskad)) +
    geom_histogram(breaks = 0:high.limit, aes(weight = duration)) +
    scale_x_discrete(limits = c(0, high.limit), breaks = 0:high.limit) +
    xlim(0, high.limit) +
    scale_y_log10() +
    facet_grid(rating.4 ~ rating.3, scales = "fixed", labeller = "label_both")

We can fit the poisson distribution using the same predictionPoisson() function as before. We also look at the fit in a simple way (the - much! - better approach is to use the residuals from the fitted model)

## rating.3 is the vehicle age.
data <-
    mccase.current[order(rating.3, antskad),
                   list(N = .N, w = sum(duration)),
                   by = list(rating.3, antskad)]
data$predicted <-
    unlist(lapply(levels(data$rating.3),
                  function (l) round(predictionPoisson(data[rating.3 == l]))))
## Show the fit
print(data[antskad <= high.limit,
           list(rating.3, antskad, N, predicted)], digits = 1)
## Show simplistic residuals per level of the rating factor
print(data[antskad <= high.limit,
           list(res = sum(abs(predicted - N))),
           by = rating.3])

## rating.4 is the bonus class
data <-
    mccase.current[order(rating.4, antskad),
                   list(N = .N, w = sum(duration)),
                   by = list(rating.4, antskad)]
data$predicted <-
    unlist(lapply(levels(data$rating.4),
                  function (l) round(predictionPoisson(data[rating.4 == l]))))
## Show the fit
print(data[antskad <= high.limit,
           list(rating.4, antskad, N, predicted)], digits = 1)
## Show simplistic residuals per rating factor
print(data[antskad <= high.limit,
           list(res = sum(abs(predicted - N))),
           by = rating.4])

We will not show the tables here (as we said: there is a better way). Instead, we can show the fitted Poisson distribution:

ggplot(data, aes(x = antskad, y = N)) +
    geom_bar(breaks = 0L:high.limit, stat = "identity") +
    facet_wrap( ~ rating.4) +
    geom_line(aes(y = predicted), data = data, colour = "red", size = 1) +
    xlim(-0.5, high.limit + 0.5)

2. Duration

We can plot the distribution of duration as before. Here 440 is the 90% quantile for the duration (412) rounded up to the next multiplier of binwidth.

plots <-
    lapply(1:4,
           function(i)
           ggplot(mccase.current, aes(duration))
           + geom_histogram(binwidth = 40)
           + xlim(0, 440)
           + scale_y_log10()
           + facet_grid(paste("rating.", i, " ~ .", sep = ""),
                        scales = "fixed")
           ## We drop the axis titles to make more room for the data
           + opts(axis.title.x = theme_blank(), axis.title.y = theme_blank(),
                  axis.text.x = theme_blank(),  axis.text.y = theme_blank(),
                  axis.title.y = theme_blank(), axis.ticks = theme_blank(),
                  title = plot.titles[i])
           )

grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow = 1, ncol = 4)))
## We can ignore the warnings from displaying the plots for now
for (i in 1:4)
    print(plots[[i]], vp = viewport(layout.pos.col = i))

However, in this case it is much more illuminating to look at the overall distribution:

ggplot(mccase, aes(duration)) +
    geom_histogram(binwidth = 0.05) +
    xlim(0, 3) +
    opts(title = "Duration of motorcycle policies") +
    annotate("text", 3.0, 1e4, label = "Histogram bin width = 0.05",
             size = 3, hjust = 1, vjust = 1)

Big spikes near 1/2 and 1 year (and in general spikes around 1/2 year multiples) suggests that either (1) this is not a random sample of policies or (2) there is a strong seasonal effect in the time of year people initially sign up for the policy.

Problem 3: Determine the relativities for claim frequency and severity

Determine the relativities for claim frequency and severity separately, by using GLMs; use the results to get relativities for the pure premium.

We first load the data (if needed) and create a data frame to hold the output, following the structure of Table 2.8 in the book.

library("data.table")
if (!exists("mccase.current"))
    load("mccase.current.RData")

case.2.4 <-
    data.frame(rating.factor =
               c(rep("Zone",        nlevels(mccase.current$rating.1)),
                 rep("MC class",    nlevels(mccase.current$rating.2)),
                 rep("Vehicle age", nlevels(mccase.current$rating.3)),
                 rep("Bonus class", nlevels(mccase.current$rating.4))),
               class =
               with(mccase.current,
                    c(levels(rating.1), levels(rating.2),
                      levels(rating.3), levels(rating.4))),
               ## These are the values from Table 2.8 in the book:
               relativity =
               c(7.678, 4.227, 1.336, 1.000, 1.734, 1.402, 1.402,
                 0.625, 0.769, 1.000, 1.406, 1.875, 4.062, 6.873,
                 2.000, 1.200, 1.000,
                 1.250, 1.125, 1.000),
               stringsAsFactors = FALSE)
print(case.2.4, digits = 3)

	rating.factor	class	relativity
1	zone	1	7.678
2	zone	2	4.227
3	zone	3	1.336
4	zone	4	1.000
5	zone	5	1.734
6	zone	6	1.402
7	zone	7	1.402
8	mc class	1	0.625
9	mc class	2	0.769
10	mc class	3	1.000
11	mc class	4	1.406
12	mc class	5	1.875
13	mc class	6	4.062
14	mc class	7	6.873
15	vehicle age	1	2.000
16	vehicle age	2	1.200
17	vehicle age	3	1.000
18	bonus class	1	1.250
19	bonus class	2	1.125
20	bonus class	3	1.000

Close enough to the book.

First we set the contrasts so the baseline for the models is the level with the highest duration. This is the same approach we used before.

library("foreach")
new.cols <-
    foreach (rating.factor = paste("rating", 1:4, sep = "."),
             .combine = rbind) %do%
{
    totals <- mccase.current[, list(D = sum(duration),
                                    N = sum(antskad),
                                    C = sum(skadkost)),
                             by = rating.factor]
    n.levels <- nlevels(mccase.current[[rating.factor]])
    contrasts(mccase.current[[rating.factor]]) <-
        contr.treatment(n.levels)[rank(-totals[["D"]], ties.method = "first"), ]
    data.frame(duration = totals[["D"]],
               n.claims = totals[["N"]],
               skadkost = totals[["C"]])
}
case.2.4 <- cbind(case.2.4, new.cols)
rm(new.cols)

We next fit the models and combine the results as in Example 2.5 above. Here we also make a note of a simple measure of the goodness of fit: the residual deviance and the degrees of freedom. The frequency model is a good fit, the severity model is not.

## Model the frequency

model.frequency <-
    glm(antskad ~
        rating.1 + rating.2 + rating.3 + rating.4 + offset(log(duration)),
        data = mccase.current[duration > 0], family = poisson)
## Res. dev. 360 on 389 dof

rels <- coef( model.frequency )
rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] )
case.2.4$rels.frequency <-
    c(c(1, rels[1:6])[rank(-case.2.4$duration[1:7], ties.method = "first")],
      c(1, rels[7:12])[rank(-case.2.4$duration[8:14], ties.method = "first")],
      c(1, rels[13:14])[rank(-case.2.4$duration[15:17], ties.method = "first")],
      c(1, rels[15:16])[rank(-case.2.4$duration[18:20], ties.method = "first")])

## Model the severity. We stick with the non-canonical link function
## for the time being.

model.severity <-
    glm(skadkost ~ rating.1 + rating.2 + rating.3 + rating.4,
        data = mccase.current[skadkost > 0,],
        family = Gamma("log"), weights = antskad)
## Res.dev. 516 on 164 dof

rels <- coef( model.severity )
rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] )
case.2.4$rels.severity <-
    c(c(1, rels[1:6])[rank(-case.2.4$duration[1:7], ties.method = "first")],
      c(1, rels[7:12])[rank(-case.2.4$duration[8:14], ties.method = "first")],
      c(1, rels[13:14])[rank(-case.2.4$duration[15:17], ties.method = "first")],
      c(1, rels[15:16])[rank(-case.2.4$duration[18:20], ties.method = "first")])

## Combine the frequency and severity
case.2.4$rels.pure.prem <- with(case.2.4, rels.frequency * rels.severity)

Finally we convert, save, and compare with the current values.

## Convert to data.table
library("data.table")
case.2.4 <- data.table(case.2.4)

## Save
save(case.2.4, file = "case.2.4.RData")

## Compare with current values
print(case.2.4[,
               list(rating.factor, class, duration, n.claims,
                    skadkostK = round(skadkost / 1e3),
                    relativity, rels.pure.prem)],
      digits = 3)

	rating.factor	class	duration	n.claims	skadkostK	relativity	rels.pure.prem
1	Zone	1	6205.310	183	5540.000	7.678	6.239
2	Zone	2	10103.090	167	4811.000	4.227	3.179
3	Zone	3	11676.573	123	2523.000	1.336	0.993
4	Zone	4	32628.493	196	3775.000	1.000	1.000
5	Zone	5	1582.112	9	105.000	1.734	0.155
6	Zone	6	2799.945	18	288.000	1.402	0.223
7	Zone	7	241.288	1	1.000	1.402	0.002
8	MC class	1	5190.351	46	993.000	0.625	0.395
9	MC class	2	3990.115	57	883.000	0.769	0.536
10	MC class	3	21665.679	166	5372.000	1.000	1.000
11	MC class	4	11739.882	98	2192.000	1.406	0.574
12	MC class	5	13439.926	149	3297.000	1.875	1.413
13	MC class	6	8880.134	175	4161.000	4.062	4.876
14	MC class	7	330.723	6	145.000	6.873	0.602
15	Vehicle age	1	4955.403	126	4964.000	2.000	4.715
16	Vehicle age	2	9753.811	145	5507.000	1.200	2.021
17	Vehicle age	3	50527.597	426	6570.000	1.000	1.000
18	Bonus class	1	19893.370	207	4558.000	1.250	0.797
19	Bonus class	2	9615.764	121	3627.000	1.125	0.604
20	Bonus class	3	35727.677	369	8857.000	1.000	1.000

It is really rather different from the current relativity, as we will also see below. But remember that the severity model was a poor fit: that would be my first place to start looking if staying within the insurance pricing framework outlined in this book.

Problem 4: Discussions

Just note here that the ratio for the relativity we calculated and for the existing one can be found as

with(case.2.4, max(rels.pure.prem) / min(rels.pure.prem)) # 3552
with(case.2.4, max(relativity) / min(relativity))         #   12

Note the huge range of relativities in our new model (3552 versus 12).

Jump to comments.

R code for Chapter 1 of Non-Life Insurance Pricing with GLM

2012-03-01T18:11:00Z

Amazon UK | US

Insurance pricing is backwards and primitive, harking back to an era before computers. One standard (and good) textbook on the topic is Non-Life Insurance Pricing with Generalized Linear Models by Esbjorn Ohlsson and Born Johansson (Amazon UK | US). We have been doing some work in this area recently. Needing a robust internal training course and documented methodology, we have been working our way through the book again and converting the examples and exercises to R, the statistical computing and analysis platform. This is part of a series of posts containing elements of the R code.

#!/usr/bin/Rscript
## PricingGLM-1.r - Code for Chapter 1 of "Non-Life Insurance Pricing with GLM"
## Copyright © 2012 CYBAEA Limited (http://www.cybaea.net/)

## @book{ohlsson2010non,
##   title={Non-Life Insurance Pricing with Generalized Linear Models},
##   author={Ohlsson, E. and Johansson, B.},
##   isbn={9783642107900},
##   series={Eaa Series: Textbook},
##   url={http://books.google.com/books?id=l4rjeflJ\_bIC},
##   year={2010},
##   publisher={Springer Verlag}
## }

With the preliminaries out of the way, let us get started.

Example 1.2

We grab the data for Table 1.2 from the book's web site and store it as an R object with lots of good meta information.

################
### Example 1.2
con <- url("http://www2.math.su.se/~esbj/GLMbook/moppe.sas")
data <- readLines(con, n = 200L, warn = FALSE, encoding = "unknown")
close(con)
## Find the data range
data.start <- grep("^cards;", data) + 1L
data.end   <- grep("^;", data[data.start:999L]) + data.start - 2L
table.1.2  <- read.table(text = data[data.start:data.end],
                         header = FALSE, sep = "", quote = "",
                         col.names = c("premiekl", "moptva", "zon", "dur",
                             "medskad", "antskad", "riskpre", "helpre", "cell"),
                         na.strings = NULL,
                         colClasses = c(rep("factor", 3), "numeric",
                             rep("integer", 4), "NULL"),
                         comment.char = "")
rm(con, data, data.start, data.end)     # Cleanup
comment(table.1.2) <-
    c("Title: Partial casco moped insurance from Wasa insurance, 1994--1999",
      "Source: http://www2.math.su.se/~esbj/GLMbook/moppe.sas",
      "Copyright: http://www2.math.su.se/~esbj/GLMbook/")
## See the SAS code for this derived field
table.1.2$skadfre = with(table.1.2, antskad / dur)
## English language column names as comments:
comment(table.1.2$premiekl) <-
    c("Name: Class",
      "Code: 1=Weight over 60kg and more than 2 gears",
      "Code: 2=Other")
comment(table.1.2$moptva)   <-
    c("Name: Age",
      "Code: 1=At most 1 year",
      "Code: 2=2 years or more")
comment(table.1.2$zon)      <-
    c("Name: Zone",
      "Code: 1=Central and semi-central parts of Sweden's three largest cities",
      "Code: 2=suburbs and middle-sized towns",
      "Code: 3=Lesser towns, except those in 5 or 7",
      "Code: 4=Small towns and countryside, except 5--7",
      "Code: 5=Northern towns",
      "Code: 6=Northern countryside",
      "Code: 7=Gotland (Sweden's largest island)")
comment(table.1.2$dur)      <-
    c("Name: Duration",
      "Unit: year")
comment(table.1.2$medskad)  <-
    c("Name: Claim severity",
      "Unit: SEK")
comment(table.1.2$antskad)  <- "Name: No. claims"
comment(table.1.2$riskpre)  <-
    c("Name: Pure premium",
      "Unit: SEK")
comment(table.1.2$helpre)   <-
    c("Name: Actual premium",
      "Note: The premium for one year according to the tariff in force 1999",
      "Unit: SEK")
comment(table.1.2$skadfre)  <-
    c("Name: Claim frequency",
      "Unit: /year")
## Save results for later
save(table.1.2, file = "table.1.2.RData")
## Print the table (not as pretty as the book)
print(table.1.2)
################

	premiekl	moptva	zon	dur	medskad	antskad	riskpre	helpre	skadfre
1	1	1	1	62.9000	18256	17	4936	2049	0.2703
2	1	1	2	112.9000	13632	7	845	1230	0.0620
3	1	1	3	133.1000	20877	9	1411	762	0.0676
4	1	1	4	376.6000	13045	7	242	396	0.0186
5	1	1	5	9.4000	0	0	0	990	0.0000
6	1	1	6	70.8000	15000	1	212	594	0.0141
7	1	1	7	4.4000	8018	1	1829	396	0.2273
8	1	2	1	352.1000	8232	52	1216	1229	0.1477
9	1	2	2	840.1000	7418	69	609	738	0.0821
10	1	2	3	1378.3000	7318	75	398	457	0.0544
11	1	2	4	5505.3000	6922	136	171	238	0.0247
12	1	2	5	114.1000	11131	2	195	594	0.0175
13	1	2	6	810.9000	5970	14	103	356	0.0173
14	1	2	7	62.3000	6500	1	104	238	0.0161
15	2	1	1	191.6000	7754	43	1740	1024	0.2244
16	2	1	2	237.3000	6933	34	993	615	0.1433
17	2	1	3	162.4000	4402	11	298	381	0.0677
18	2	1	4	446.5000	8214	8	147	198	0.0179
19	2	1	5	13.2000	0	0	0	495	0.0000
20	2	1	6	82.8000	5830	3	211	297	0.0362
21	2	1	7	14.5000	0	0	0	198	0.0000
22	2	2	1	844.8000	4728	94	526	614	0.1113
23	2	2	2	1296.0000	4252	99	325	369	0.0764
24	2	2	3	1214.9000	4212	37	128	229	0.0305
25	2	2	4	3740.7000	3846	56	58	119	0.0150
26	2	2	5	109.4000	3925	4	144	297	0.0366
27	2	2	6	404.7000	5280	5	65	178	0.0124
28	2	2	7	66.3000	7795	1	118	119	0.0151

That was easy. Now for something a little harder.

Example 1.3

Here we are concerned with replicating Table 1.4. We do it slowly, step-by-step, for pedagogical reasons.

################
### Example 1.3
if (!exists("table.1.2"))
    load("table.1.2.RData")
## We calculate each of the columns individually and slowly here
## to show each step

## First we have simply the labels of the table
rating.factor <-
    with(table.1.2,
         c(rep("Vehicle class", nlevels(premiekl)),
           rep("Vehicle age", nlevels(moptva)),
           rep("Zone", nlevels(zon))))

## The Class column
class.num <- with(table.1.2, c(levels(premiekl), levels(moptva), levels(zon)))

## The Duration is the sum of durations within each class
duration.total <-
    c(with(table.1.2, tapply(dur, premiekl, sum)),
      with(table.1.2, tapply(dur, moptva, sum)),
      with(table.1.2, tapply(dur, zon, sum)))

## Calculate relativities in the tariff
## The denominator of the fraction is the class with the highest exposure
## (i.e. the maximum total duration): we make that explicit with the
## which.max() construct.  We also set the contrasts to use this as the base,
## which will be useful for the glm() model later.
class.base <- which.max(duration.total[1:2])
age.base   <- which.max(duration.total[3:4])
zone.base  <- which.max(duration.total[5:11])

rt.class <- with(table.1.2, tapply(helpre, premiekl, sum))
rt.class <- rt.class / rt.class[class.base]
rt.age   <- with(table.1.2, tapply(helpre, moptva, sum))
rt.age   <- rt.age / rt.age[age.base]
rt.zone  <- with(table.1.2, tapply(helpre, zon, sum))
rt.zone  <- rt.zone / rt.zone[zone.base]

contrasts(table.1.2$premiekl) <-
    contr.treatment(nlevels(table.1.2$premiekl))[rank(-duration.total[1:2],
                                                      ties.method = "first"), ]
contrasts(table.1.2$moptva) <-
    contr.treatment(nlevels(table.1.2$moptva))[rank(-duration.total[3:4],
                                                    ties.method = "first"), ]
contrasts(table.1.2$zon) <-
    contr.treatment(nlevels(table.1.2$zon))[rank(-duration.total[5:11],
                                                 ties.method = "first"), ]

The contrasts could also have been set with the base= argument, e.g. contrasts(table.1.2$zon) <- contr.treatment(nlevels(table.1.2$zon), base = zone.base), which would be closer in spirit to the SAS code. But I like the idiom presented here where we follow the duration order; it also extends well to other (i.e. not treatment) contrasts. I just wish rank() had an decreasing= argument like order() which I think would be clearer than using rank(-x) to get a decreasing sort order.

That was the easy part. At this stage in the book you are not really expected to understand the next step so do not despair! We just show how easy it is to replicate the SAS code in R. An alternative approach using direct optimization is outlined in Exercise 1.3 below.

## Relativities of MMT; we use the glm approach here as per the book's
## SAS code at http://www2.math.su.se/~esbj/GLMbook/moppe.sas
m <- glm(riskpre ~ premiekl + moptva + zon, data = table.1.2,
         family = poisson("log"), weights = dur)

## If the next line is a mystery then you need to
## (1) read up on contrasts or
## (2) remember that the link function is log() which is why we use exp here
rels <- exp( coef(m)[1] + coef(m)[-1] ) / exp(coef(m)[1])

rm.class <- c(1, rels[1])               # See rm.zone below for the
rm.age   <- c(rels[2], 1)               # general approach
rm.zone  <- c(1, rels[3:8])[rank(-duration.total[5:11], ties.method = "first")]

## Create and save the data frame
table.1.4 <-
    data.frame(Rating.factor = rating.factor, Class = class.num,
               Duration = duration.total,
               Rel.tariff = c(rt.class, rt.age, rt.zone),
               Rel.MMT    = c(rm.class, rm.age, rm.zone))
save(table.1.4, file = "table.1.4.RData")
print(table.1.4, digits = 3)
rm(rating.factor, class.num, duration.total, class.base, age.base, zone.base,
   rt.class, rt.age, rt.zone, rm.class, rm.age, rm.zone, m, rels)
################

The result is something like this:

	Rating.factor	Class	Duration	Rel.tariff	Rel.MMT
1	Vehicle class	1	9833.20	1.00	1.00
2	Vehicle class	2	8825.10	0.50	0.43
3	Vehicle age	1	1918.40	1.67	2.73
4	Vehicle age	2	16739.90	1.00	1.00
5	Zone	1	1451.40	5.17	8.97
6	Zone	2	2486.30	3.10	4.19
7	Zone	3	2888.70	1.92	2.52
8	Zone	4	10069.10	1.00	1.00
9	Zone	5	246.10	2.50	1.24
10	Zone	6	1369.20	1.50	0.74
11	Zone	7	147.50	1.00	1.23

Note the rather unusual and apparently inconsistent rounding in the book: 147, 1.66, and 5.16 would be better as 148 (the value is 147.5), 1.67, and 5.17.

Exercise 1.3

Here it gets interesting as we get a different value from the authors. Possibly a small bug on our part but at least we provide the code for you to check. So if you spot a problem let us know in the comments.

################
## Exercise 1.3

## The values from the book
g0  <- 0.03305
g12 <- 2.01231
g22 <- 0.74288
dim.names <- list(Milage = c("Low", "High"),
                  Age = c("New", "Old"))
pyears <- matrix(c(47039, 56455, 190513, 28612), nrow = 2,
                 dimnames = dim.names)
claims <- matrix(c(0.033, 0.067, 0.025, 0.049), nrow = 2,
                 dimnames = dim.names)

## Function to calculate the error of the estimate
GvalsError <- function (gvals) {
    ## The current estimates
    g0  <- gvals[1]
    g12 <- gvals[2]
    g22 <- gvals[3]
    ## The current estimates in convenient matrix form
    G  <- matrix(c(1, 1, g12, g22), nrow = 2)
    G1 <- matrix(c(1, g12), nrow = 2, ncol = 2)
    G2 <- matrix(c(1, g22), nrow = 2, ncol = 2, byrow = TRUE)
    ## The calculated values
    G0  <- addmargins(claims * pyears)["Sum", "Sum"] / ( sum(pyears * G1 * G2) )
    G12 <- addmargins(claims * pyears)["High", "Sum"] /
        ( g0 * addmargins(pyears * G2)["High", "Sum"] )
    G22 <- addmargins(claims * pyears)["Sum", "Old"] /
        ( g0 * addmargins(pyears * G1)["Sum", "Old"] )
    ## The sum of squared errors
    error <- (g0 - G0)^2 + (g12 - G12)^2 + (g22 - G22)^2
    return(error)
}

## Minimize the error function to obtain our estimate
gamma <- optim(c(g0, g12, g22), GvalsError)
stopifnot(gamma$convergence == 0)
gamma <- gamma$par

values <- data.frame(legend = c("Our calculation", "Book value"),
                     g0  = c(gamma[1], g0),
                     g12 = c(gamma[2], g12),
                     g22 = c(gamma[3], g22),
                     row.names = "legend")
print(values, digits = 4)

## Close, but not the same.

rm(g0, g12, g22, dim.names, pyears, claims, gamma, values)
################

The resulting table is something like:

	g0	g12	g22
Our calculation	0.0334	1.9951	0.7452
Book value	0.0331	2.0123	0.7429

Close, but not the same. Perhaps they used a different error function.

Jump to comments.

doSMP pulled

2012-03-01T09:16:00Z

They have finally pulled that buggy unreliable piece of code that was doSMP from the CRAN mirrors while (I hear) Revolutions are re-writing it. To use all your cores for analysis on the Windows platform, you can try doSNOW instead; my code is something like the fragment below. Neither option is as attractive as doMC on anything-but-Windows platforms, but sometimes you have to work with legacy systems.


library("foreach")
if (.Platform$OS.type != "windows" && require("multicore")) {
    registerDoMC()
} else if (FALSE &&                     # doSMP is buggy
           require("doSMP")) {
    w <- startWorkers()
    on.exit(stopWorkers(w), add = TRUE)
    registerDoSMP(w)
} else if (require("doSNOW")) {
    cl <- snow::makeCluster(4, type = "SOCK")
    on.exit(snow::stopCluster(cl), add = TRUE)
    registerDoSNOW(cl)
} else {
    registerDoSEQ()
}

Change the number 4 to the number of cores that you want to use on the machine. The explicit name space (snow::) is to avoid confusion if you load the "parallel" package or any of the other packages that also define a makeCluster() function.

I hope Revolutions does a good job on the new version: it needs some love.

R versus SAS/SPSS in corporations

2011-10-28T11:10:00Z

A recent question on one of the LinkedIn groups about the advantages of using R over commercial tools like SAS or IBM SPSS Modeller drew lots of comments for R. We like R a lot and we use it extensively, but I also wanted to balance the discussion. R is great, but looking at commercial organizations near the end of 2011 it is not necessarily the right choice to make.

Background

We have created and managed analytics teams in commercial organizations (mainly telecommunications) across Europe. The teams were using SAS or SPSS. Our company now has a commercial analytics as a service offering and we mainly use R.

The benefits of R is productivity. We want to spend time on the actions from the analytical insights, not the coding, and we choose our tool accordingly. Being a consulting type organization it is easier for us to attract and retain talent.

The advantages of SAS/SPSS in a commercial environment

1. You can buy the tool for money.

Big corporations have procurement departments who do not have a process for free software. Also software spend goes on the balance sheet in a way that the CFO prefers to people but something like R will take a little talent to set up initally. (And yes, we know the Revolutions guys well, but they are not really credible in Europe yet.)

This will change as (a) companies become more mature in their procurement and as (b) commercial support for R improves. (On the latter point, Oracle’s R integration to the database is great news.)

2. You can recruit for the commercial tool

Recruiters are familiar with SAS and SPSS but not with R so it is easier to brief them and to get good quality CVs. This will change and R becomes ever more popular and prevalent. [And yes, we could in theory change recruiters to someone clued in, but again in large corporations there are procurement processes to be followed and existing agreements to be honoured so it will all take months or years.]
There are recognised training programmes for SPSS and (especially) SAS which makes it easier to recruit the technical skills. How do you know what somebody knows when they say they “know R”? How do you even begin to quantify it from a CV? How do you separate the guy who downloaded the tool and just read “An Introduction to R” from the Frank Harrells of this world?

Yes, I would argue (and in fact have argued in Commercial Analytics: The Capabilities) that technical skill is not the most important in an analyst (and can be learned anyhow) but it does help filter the CVs and, you guessed it, fits well with the corporation’s processes.

(One reason we use R internally is that we find that it is, on average, a more interesting type of analyst who is proficient in that tool. It seems to encourage curiosity and love or learning in a way that menu-based tools do not.)

I think the commercial R companies are really missing a trick here to provide recognised certification.
You can’t search for R. Seriously: try searching for R on LinkedIn (tip: there is another way). Much easier to find SAS / SPSS skills in a large CV database (like LinkedIn where this discussion started).

3. You can recruit for the commercial tool.

Yes I know I already said that but there is another reason why this is critical.

R takes talent to use. (That is kind of why we like it.) It takes talent to maintain.

My problem as the manager of a commercial analytical insights team is that it is very hard for me to retain that talent. Think about it: what can I offer in terms of career progression? If you are an analyst you might become a senior analyst but you will always be an analyst. There are no examples of a way up the organization (except perhaps out through IT and then up to CIO). [This too will change with time.] And new challenges: yes, some, but we are not a research university and it tends to be the same few problem types that we are always working on. So if you are an analyst looking for new challenges and more pay, the best thing – the logical and rational thing to do – is to get a new job. And your time with Big Corporation will look good on your CV and you will probably land the job easily.

We can help

Our commercial analytics capabilities model.

If you want to set up a commercial analytical group we can help you get it right first time. The right people, the right processes, the right infrastructure and most importantly the right results. We have done it before and are not tied to any specific tool or vendor.
If you want to improve or enhance your existing analytical teams, then we can Reboot your Analytics to deliver both rapid and sustained commercial results.
And if you just want the results we can provide commercial analytics as a service where we provide the insights and then work with you to turn those insights into commercial actions and better understanding of your business, markets, and customers, leaving you to focus on what you do best.

Jump to comments.

Friday quote: what is the question to which this number is the answer?

2011-08-26T09:05:00Z

John Kay muses on interpreting statistical data:

Always ask of such data “what is the question to which this number is the answer?”. “Earnings before interest, tax, depreciation and amortisation on a like-for-like basis before allowance for exceptional restructuring costs” is the answer to the question “what is the highest profit number we can present without attracting flat disbelief?”.

And on the pitfalls of powerful data analysis tools:

When the data seem to point to an unexpected finding, always consider the possibility that the problem is a feature of the data, rather than a feature of the world. […] It is now easy to import data into a computer program without thought. The unwarranted precision of the projected growth in rail traffic – a 96 per cent increase, rather than a doubling – is a clue that the number was generated by a computer, not a skilled interpreter of evidence.

Statistics are only as valid as the sources from which they are drawn and the abilities of those who use them. When I discover something surprising in data, the most common explanation is that I made a mistake.

A warning on the R save format

2011-08-23T07:20:00Z

The save() function in the R platform for statistical computing is very convenient and I suspect many of us use it a lot. But I was recently bitten by a “feature” of the format which meant I could not recover my data.

How to lose your data with `save()`

I am using Windows on my travel laptop and Linux on my workstation. To speed things up on the latter and make use of my many (well, four) cores, I use the ‘multicore’ package, which I do not have available on the Windows machine.

To illustrate the problem with the save file format, I created a file on the Linux machine simply as:

library("multicore")
a <- list(data = 1:10, fun = mclapply)
save(a, file = "a.RData")

What could be simpler? The mclapply is a function from the ‘multicore’ package but it clearly has no impact on the stored data. (We will show a more realistic example below – work with me here.)

But try to open the save file on a machine without the package installed, like my Windows laptop, and you get:

Error in loadNamespace(name) : there is no package called 'multicore'

There is no way of getting to your precious data without installing the missing package.

If the package has been withdrawn or is no longer available then your data is basically lost.

What can you do?

Some suggestions from the helpful people on R-help:

(Uwe Ligges): You could try to rewrite ./src/main/saveload.R and serialize.R to extract only the parts you need. “This is probably not worth the effort.”
(Prof. Brian Ripley): You could try installing the missing package; R CMD INSTALL --fake should be sufficient to let you load the data. Also suggests that the proposal above would be very hard indeed.
(Martin Morgan): Don't store package functions with your code.

That is three good answers from three of the heavy-weights in the R community. Thank you all!

Martin’s comment is worth expanding. We can change the above example to:

library("multicore")
computeFunction <- function(...) {
    if (require(multicore)) mclapply(...)
    else lapply(...) 
}
a <- list(data = 1:10, fun = computeFunction)
save(a, file = "a.RData")

Now everything works fine! No data is horribly lost: the file loads fine on the ‘multicore’-less machine.

And for the more realistic example, I had been using caret::rfe as Martin knew in the example he provided:

library("caret")
data(BloodBrain)

x <- scale(bbbDescr[,-nearZeroVar(bbbDescr)])
x <- x[, -findCorrelation(cor(x), .8)]
x <- as.data.frame(x)

set.seed(1)
lmProfile <- rfe(x, logBBB,
                 sizes = c(2:25, 30, 35, 40, 45, 50, 55, 60, 65),
                 rfeControl = rfeControl(functions = lmFuncs,
                   number = 5,
                   computeFunction=mclapply))
save(lmProfile, file = "lmProfile.RData")

Slightly less obvious that there is a reference to the external namespace in this code, but easy enough to see if you know what to look for.

For old files I will use the R CMD INSTALL --fake suggestion, but for new data I am going with the last approach and using a computeFunction like this:

### MCCompute: A computeFunction for caret::rfeControl and caret::trainControl 
### that does not leave a reference to the multicore package in the save file
MCCompute <- function(X, FUN, ...) {
    FUN <- match.fun(FUN)
    if (!is.vector(X) || is.object(X)) 
        X <- as.list(X)
    if (require("multicore")) mclapply(X, FUN, ...)
    else lapply(X, FUN, ...)
}

I know that Max Kuhn is rewriting the caret package which should make this a moot point in the near future for that specific case. But the indirection approach is generally useful and will also be relevant in other situations.

Recommendations

My recommendations:

Save data in a data format, not using the save() function which is really for objects (data and code). Suitable formats include CSV and variants, HDF5, and CDF, as well as others.
Avoid references to packages in your objects by using the one level indirection trick exemplified by the MCCompute function shown.

What is your approach? Suggestions in the comments below, please.

Jump to comments.

Friday quote: the handmaiden and the whore

2011-08-19T12:04:00Z

Because it is Friday and because we collect quotes.

If mathematics is the handmaiden of science, statistics is the whore: all that scientists are looking for is a quick fix without the encumbrance of a meaningful relationship. Statisticians are second-class mathematicians, third-rate scientists and fourth-rate thinkers. They are the hyenas, jackals and vultures of the scientific ecology: picking over the bones and carcasses of the game that the big cats, the biologists, the physicists and the chemists, have brought down.

Statistics is a wonderful discipline. It has it all: mathematics and philosophy, analysis and empiricism, as well as applicability, relevance and the fascination of data. It demands clear thinking, good judgement and flair. Statisticians are engaged in an exhausting but exhilarating struggle with the biggest challenge that philosophy makes to science: how do we translate information into knowledge?

―Stephen Senn: Dicing with Death: Chance, Risk and Health

Which one of the two views are closest to your opinion?

Spreadsheet errors

2011-04-20T11:19:00Z

For my sins, I have done more than my fair share of analysis in Excel. I am quite capable of building and maintaining 130Mb spreadsheets (I had a dozen of them for one client). Excel is pretty much installed everywhere, so it is sometimes the only way to get started getting commercial value of the data in the organisation. But I don’t like it and let’s have a look at one reason why. In order not to always pick on Microsoft, we use another application, but you get the same results with Excel.

Y	X1	X2	X3	X4
5.88	1	1	1	1
2.56	6	1	1	1
11.11	1	1	1	1
0.79	6	1	1	1
0.00	6	1	1	1
0.00	0	1	1	1
15.6	8	1	1	1
3.7	4	1	1	1
8.49	3	1	1	1
51.2	6	1	1	1
14.2	7	1	1	1
7.14	5	1	1	1
4.2	7	1	1	1
6.15	4	1	1	1
10.46	6	1	1	1
0.00	8	1	1	1
10.42	2	1	1	1
17.36	5	1	1	1
13.41	8	1	1	1
41.67	0	1	1	1
2.78	0	1	1	1
2.98	8	1	1	1
9.62	7	1	1	1
0.00	0	1	1	1
4.65	5	1	0	2
3.13	3	1	0	2
24.58	6	1	0	2
0.00	1	1	0	2
5.56	4	1	0	2
9.26	3	1	0	2
0.00	0	1	0	2
0.00	0	1	0	2
3.13	1	1	0	2
0.00	0	1	0	2
7.56	5	0	1	3
9.93	6	0	1	3
0.00	8	0	1	3
16.67	6	0	1	3
16.89	7	0	1	3
13.71	6	0	1	3
6.35	5	0	1	3
2.5	3	0	1	3
2.47	7	0	1	3
21.74	3	0	1	3
23.6	8	0	0	4
11.11	8	0	0	4
0.00	7	0	0	4
3.57	8	0	0	4
2.9	5	0	0	4
2.94	3	0	0	4
2.42	8	0	0	4
18.75	4	0	0	4
0.00	5	0	0	4
2.27	3	0	0	4

Spreadsheets are good for some things, but analysing data is not one of them. The example data in the table on the right is from Jeffrey S. Simonoff, “Statistical analysis using Microsoft Excel” (2008), and looks at first (and maybe even second) glance like a reasonable set of observations.

However, the predictors are (accidentally) collinear so no meaningful fit is possible, unless one of them are dropped. We see that very easily if we try to do the analysis using the R statistical computing and analysis platform:

> d <- read.delim("clipboard")  # Read DATA range from clipboard
> summary(lm(Y ~ ., data = d))

Call:
lm(formula = Y ~ ., data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.222  -5.821  -2.546   3.171  40.750 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.1945     3.9749   1.055    0.296
X1            0.3862     0.5652   0.683    0.497
X2            0.2308     3.1590   0.073    0.942
X3            3.7072     2.9922   1.239    0.221
X4                NA         NA      NA       NA

Residual standard error: 10.14 on 50 degrees of freedom
Multiple R-squared: 0.04767,	Adjusted R-squared: -0.009466 
F-statistic: 0.8343 on 3 and 50 DF,  p-value: 0.4814

We have highlighted the message that R has automatically dropped one of the predictors.

Everybody likes to pick on Excel, so let us load the data into version 3.3.2 of LibreOffice, the free Open Source personal productivity suite, instead. It faithfully implements many of the worst features of Excel. You can grab a copy of the spreadsheet GS-spreadsheet-error.ods yourself and see the results. The relevant function in both Excel and LibreOffice for linear regression is LINEST and applying it to the data set give us:

Of the 16 values returned by the function, fully 12 of them are incorrect (highlighted in red), and the '#VALUE!' entries are the only thing that suggests we may have a problem. (The '#N/A' values are a feature of the function and not a problem.) Excluding the X4 values from the function call gives meaningful (and correct) results:

There is so much wrong with doing even this trivial analysis in a spreadsheet that it is hard to know where to start. Some of the problems:

Garbage results instead of errors: Instead of giving meaningful errors or warnings, the spreadsheets simply produce garbage results. This is nearly impossible to debug.
No help on how to correct the problem: In the erroneous results of the first figure, there is no clue, no hint, no help to figure out how to correct the problem. You could argue about R correcting the issue ”automagically”, but at least it finds a solution to the problem and tells you about it.
Error prone output formats: I put in the row and column headings because otherwise it is just too hard to read the data. Where does the function stuff the F statistics again?

And don’t get me started on version control and documentation. Don’t even mention that the maths in Excel are wrong. Remember: Friends do not let friends do data analysis in spreadsheets.

Jump to comments.

Getting started with the Heritage Health Price competition

2011-04-08T08:39:00Z

The US$ 3 million Heritage Health Price competition is on so we take a look at how to get started using the R statistical computing and analysis platform.

We do not have the full set of data yet, so this is a simple warm-up session to predict the days in hospital in year 2 based on the year 1 data.

Prerequisites

Obviously you need to have R installed, and you should also have signed up for the competition (be sure to read the terms carefully) and downloaded and extracted the release 1 data file.

Data preparation

Let’s load the data into R and do some basic housekeeping:

#!/usr/bin/Rscript
## example001.R - simple benchmarks for the HHP
## Copyright © 2011 CYBAEA Limited - http://www.cybaea.net/

##############################
#### DATA PREPARATION

##++++
## Members
members   <- read.csv(file = "HHP_release1/Members_Y1.csv",
                      colClasses = rep("factor", 3),
                      comment.char = "")
##----
##++++
## Claims
claims.Y1 <- read.csv(file = "HHP_release1/Claims_Y1.csv",
                      colClasses = c(
                          rep("factor", 7),
                          "integer",    # paydelay
                          "character",  # LengthOfStay
                          "character",  # dsfs
                          "factor",     # PrimaryConditionGroup
                          "character"   # CharlsonIndex
                          ),
                      comment.char = "")
## Utility function
make.numeric <- function (cv, FUN = mean) {
### make a character vector numeric by splitting on '-'
    sapply(strsplit(gsub("[^[:digit:]]+",
                         " ",
                         cv,
                         perl = TRUE),
                    " ",
                    fixed = TRUE),
           function (x) FUN(as.numeric(x)))
}
## Length of stay as days
{
    z <- make.numeric(claims.Y1$LengthOfStay)
    z.week <- grepl("week", claims.Y1$LengthOfStay, fixed = TRUE)
    z[z.week] <- z[z.week] * 7          # Weeks are 7 days
    z[is.nan(z)] <- 0
    claims.Y1$LengthOfStay.days <- z
}
los.levels <- c("", "1 day", sprintf("%d days", 2:6),
                "1- 2 weeks", "2- 4 weeks", "4- 8 weeks", "8-12 weeks",
                "12-26 weeks", "26+ weeks")
stopifnot(all(claims.Y1$LengthOfStay %in% los.levels))
claims.Y1$LengthOfStay <- factor(claims.Y1$LengthOfStay,
                                 levels = los.levels,
                                 labels = c("0 days", los.levels[-1]),
                                 ordered = TRUE)
## Months since first claim
claims.Y1$dsfs.months <- make.numeric(claims.Y1$dsfs)
## dsfs is an ordered factor and gives the ordering of the claims
dsfs.levels <- c("0- 1 month", sprintf("%d-%2d months", 1:11, 2:12))
claims.Y1$dsfs <- factor(claims.Y1$dsfs, levels = dsfs.levels, ordered = TRUE)
## Index as numeric
claims.Y1$CharlsonIndex.numeric <- make.numeric(claims.Y1$CharlsonIndex)
claims.Y1$CharlsonIndex <- factor(claims.Y1$CharlsonIndex, ordered = TRUE)
##----
##++++
## Days in hospital
dih.Y2    <- read.csv(file = "HHP_release1/DayInHospital_Y2.csv",
                      colClasses = c("factor", "integer"),
                      comment.char = "")
names(dih.Y2)[1] <- "MemberID"          # Fix broken file
##----
save(members, claims.Y1, dih.Y2,
     file = "HHPR1.RData")
##############################

Scoring

We will need a function to score our predictions p against the actual values a. The formula is on the evaluation page and we implement it as:

#!/usr/bin/Rscript
## example001.R - simple benchmarks for the HHP
## Copyright © 2011 CYBAEA Limited - http://www.cybaea.net/

##############################
#### FUNCTION TO CALCULATE SCORE
HPPScore <- function (p, a) {
### Scorng function after
### http://www.heritagehealthprize.com/c/hhp/Details/Evaluation
### Base 10 log from http://www.heritagehealthprize.com/forums/default.aspx?g=posts&m=2226#post2226
    sqrt(mean((log(1+p, 10) - log(1+a, 10))^2))
}
##############################

The simplest benchmarks

The simplest models don’t really model at all: they just use the average and are simple benchmarks.

#!/usr/bin/Rscript
## example001.R - simple benchmarks for the HHP
## Copyright © 2011 CYBAEA Limited - http://www.cybaea.net/

y <- dih.Y2$DaysInHospital_Y2           # Actual
p <- rep(mean(y), NROW(dih.Y2))
cat(sprintf("Score using mean  : %8.6f\n", HPPScore(p, y)))
# Score using mean  : 0.278725

p <- rep(median(y), NROW(dih.Y2))
cat(sprintf("Score using median: %8.6f\n", HPPScore(p, y)))
# Score using median: 0.267969

Simple single-variable linear models

OK, a model that doesn’t use past data isn’t much of a model, so let’s improve on that:

#!/usr/bin/Rscript
## example001.R - simple benchmarks for the HHP
## Copyright © 2011 CYBAEA Limited - http://www.cybaea.net/
library("reshape2")

vars <- dcast(claims.Y1, MemberID ~ ., sum, value_var = "LengthOfStay.days")
names(vars)[2] <- "LengthOfStay"
data <- merge(vars, dih.Y2)

model <- lm(DaysInHospital_Y2 ~ LengthOfStay, data = data)
p <- predict(model)
cat(sprintf("Score using lm(LengthOfStay): %8.6f\n", HPPScore(p, y)))
# Score using lm(LengthOfStay): 0.279062

model <- glm(DaysInHospital_Y2 ~ LengthOfStay,
             family = quasipoisson(),
             data = data)
p <- predict(model, type="response")
cat(sprintf("Score using glm(LengthOfStay): %8.6f\n", HPPScore(p, y)))
# Score using glm(LengthOfStay): 0.278914

Let the competition begin.

Jump to comments.

Benchmarking feature selection with Boruta and caret

2010-11-25T13:43:00Z

Feature selection is the data mining process of selecting the variables from our data set that may have an impact on the outcome we are considering. For commercial data mining, which is often characterised by having too many variables for model building, this is an important step in the analysis process. And since we often work on very large data sets the performance of our process is very important to us.

Having looked at feature selection using the Boruta package and feature selection using the caret package separately, we now consider the performance of the two approaches.

For our tests we will use an artificially constructed trivial data sets that the automated process should have no problems with (but we will be disappointed later on this expectation, as we will see). The data set has an equal number of normal and uniform random variables with mean 0 and variance 1 of which 20% are used for the target variable. There are 10 time as many observations as variables. We create a function to set this up:

make.data <- function (n.var, m.rand = 5, m.obs = 10) {
    n.col <- n.var * m.rand
    n.obs <- n.col * m.obs * 2
    x <- data.frame(N = matrix(rnorm(n = n.col*n.obs),
                        nrow = n.obs, ncol = n.col),
                    U = matrix(runif(n = n.col*n.obs,
                        min = -sqrt(3), max = sqrt(3)), n.obs, n.col))
    deps.n <- 1:n.var
    deps.u <- (1+n.col):(n.var+n.col)
    y <- rowSums(as.matrix(x[, c(deps.n, deps.u)]))
    x <- cbind(x, Y = factor(y >= 0, labels=c("N", "P")))
    attr(x, "vars") <- names(x)[c(deps.n, deps.u)]
    return(x)
}

The Boruta package

Then we run a test using the Boruta package for different sizes:

#!/usr/bin/Rscript
## bench.R - benchmark Boruta package
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)
run.name <- "bench-1"
library("Boruta")

set.seed(1)

sizes <- c(1:10, 10*(2:10), 100*(2:10), 1e3*(2:10))
n.sizes <- length(sizes)
bench <- data.frame(n.vars = sizes, elapsed = NA, right = NA, wrong = NA)
file.name <- paste(run.name, "RData", sep = ".")

for (n in 1:length(sizes)) {
    size <- sizes[n]
    cat(sprintf("[%s] Size = %3d: ", as.character(Sys.time()), size))
    tries <- max(3, round(10/size, 0))
    n.right <- 0
    n.wrong <- 0
    elapsed <- 0
    for (try in 1:tries) {
        cat(tries-try, ".", sep = "")
        x <- make.data(size)
        x.vars <- attr(x, "vars")
        elapsed <- elapsed +
            system.time({b <- Boruta(x[,-NCOL(x)], x[,NCOL(x)])}
                        )["elapsed"]
        b.vars <- names(b$finalDecision)[b$finalDecision!="Rejected"]
        n.right <- n.right + length(intersect(b.vars, x.vars))
        n.wrong <- n.wrong + length(setdiff(b.vars, x.vars))
    }
    elapsed <- elapsed / tries
    cat(" Elapsed = ", round(elapsed, 0), " seconds\n", sep = "")
    n.right <- n.right / tries
    n.wrong <- n.wrong / tries
    bench[n, ] <- c(size, elapsed, n.right, n.wrong)
    save(bench, file = file.name, ascii = FALSE, compress = FALSE)
}

print(bench)

As it turned out, our expectations for the size of data set we could handle were wildly optimistic and we killed the process at size 30. We add to the data set a field with the total number of variables in the x data set and plot the results.

load(file = "bench-1.RData")
bench <- na.omit(bench)
bench$n.elem <- bench$n.var^2 * 1e3
plot(elapsed ~ n.elem, data = bench, type = "b",
     main = "Feature selections with Boruta",
     sub = "Elapsed time versus number of data elements",
     log = "xy",
     xlab = "Elements in data set", ylab = "Elapsed time (seconds)")

Benchmarking results for feature selection with Boruta package shows linear scaling (slope is 1.01 with standard error 0.025 and adjusted R² 0.993)

A quick check using summary(lm(log(elapsed) ~ log(n.elem), data = bench)) shows us a linear scaling with the number of elements (slope is 1.01 with standard error 0.025 and adjusted R² 0.993). The algorithm selects all the right features up to n.vars = 10 when it starts to miss some of them:

Benchmark results for Boruta package
n.vars	right	wrong
1	2.00000	1.1000000
2	4.00000	1.2000000
3	6.00000	1.6666667
4	8.00000	1.3333333
5	10.00000	1.6666667
6	12.00000	1.3333333
7	14.00000	1.0000000
8	16.00000	1.3333333
9	18.00000	0.6666667
10	20.00000	0.3333333
20	39.33333	0.0000000
30	56.33333	0.0000000

A higher accuracy in the feature selection for the larger problems could presumably be achieved by adjusting the maxRuns and perhaps confidence parameters on the Boruta call.

In summary, the Boruta package performs well up to about 20 features out of 100 (n.vars = 10) which runs in about 11 minutes on my machine. If we changed the technical implementation to support multicore, MPI, and other parallel frameworks, then the out of the box settings would be useful up to n.vars of 20 or 30 (40-60 features out of 200-300) which an 8-core machine should be able to complete in 20 minutes or so.

This is still a lot less than the size of data sets we normally work with. (Our usual benchmark is 15,000 variables and 50,000 observations.)

The caret package

One of the nice features of the caret package is that is supports most parallel processing frameworks out of the box, but for comparison with the previous analysis we will (somewhat unfairly) not use that here. The setup is then quite simple, using the same make.data function as before.

#!/usr/bin/Rscript
## bench.R - benchmark caret package
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)
run.name <- "bench-2"
library("caret")
library("randomForest")
set.seed(1)

control <- rfeControl(functions = rfFuncs, verbose = FALSE,
                      returnResamp = "final")

## if ( require("multicore", quietly = TRUE, warn.conflicts = FALSE) ) {
##     control$workers <- multicore:::detectCores()
##     control$computeFunction <- mclapply
##     control$computeArgs <- list(mc.preschedule = FALSE, mc.set.seed = FALSE)
## }

our.sizes <- c(2:10, 10*(2:10), 100*(2:10), 1e3*(2:10))
n.sizes <- length(our.sizes)
bench <- data.frame(n.vars = our.sizes, elapsed = NA, right = NA, wrong = NA)
file.name <- paste(run.name, "RData", sep = ".")

for (n in 1:length(our.sizes)) {
    size <- our.sizes[n]
    cat(sprintf("[%s] Size = %3d: ", as.character(Sys.time()), size))
    tries <- max(3, round(10/size, 0))
    n.right <- 0
    n.wrong <- 0
    elapsed <- 0
    for (try in 1:tries) {
        cat(tries-try, ".", sep = "")
        x <- make.data(size)
        x.vars <- attr(x, "vars")
        elapsed <- elapsed + 
            system.time({p <- rfe(x[,-NCOL(x)], x[,NCOL(x)],
                                  sizes = 1:(2*size), rfeControl = control)}
                        )["elapsed"]
        p.vars <- predictors(p)
        n.right <- n.right + length(intersect(p.vars, x.vars))
        n.wrong <- n.wrong + length(setdiff(p.vars, x.vars))
    }
    elapsed <- elapsed / tries
    cat(" Elapsed = ", round(elapsed, 0), " seconds\n", sep = "")
    n.right <- n.right / tries
    n.wrong <- n.wrong / tries
    bench[n, ] <- c(size, elapsed, n.right, n.wrong)
    save(bench, file = file.name, ascii = FALSE, compress = FALSE)
}

print(bench)

This uses the randomForest classifier from the package of the same name. To use the ipredbagg bagging classifier from Andrea Peters and Torsten Hothorn's ipred: Improved Predictors package we simply change the control object to:

control <- rfeControl(functions = treebagFuncs, verbose = FALSE,
                      returnResamp = "final")

As usual, we were widely optimistic in our guesses for the size of problems we could handle, and had to abort the run.

Benchmarking results for feature selection with caret package using randomForest classifier (slope is 1.17 with standard error 0.024 and adjusted R² 0.996)

Benchmarking results for feature selection with caret package using treebag classifier shows non-power behaviour (nevertheless, a linear log-log fit gives a slope of 1.12 with standard error 0.067 and adjusted R² 0.96)

Benchmark results for caret package using randomForest classifier
n.vars	right	wrong
2	3.20000	3.200000
3	5.00000	0.000000
4	7.00000	0.000000
5	9.00000	0.000000
6	11.00000	0.000000
7	13.00000	0.000000
8	14.66667	0.000000
9	16.66667	0.000000
10	19.00000	0.000000
20	38.66667	1.333333
30	54.00000	86.000000

Benchmark results for caret package using ipredbagg classifier
n.vars	right	wrong
2	3.00000	0
3	5.00000	0
4	7.00000	0
5	9.00000	0
6	10.33333	0
7	13.00000	0
8	14.33333	0
9	16.00000	0
10	18.66667	0
20	35.33333	0
30	54.33333	0
40	69.66667	0

Remember that the right number of significant features are 2 * n.vars and we see that the caret package apparently always miss one feature in its selection, which is very odd and possibly a bug. It is less likely to select the wrong features than Boruta, but that could be partially due to "Tentative" data in Boruta. Timing-wise, performance is a little worse in the non-parallel situation but realistically of course a lot better than Boruta depending on the number of cores on your processor or nodes in your cluster.

Neither approach is suitable out of the box for the sizes of data sets that we normally work with.

Jump to comments.

Feature selection: Using the caret package

2010-11-16T19:35:00Z

Feature selection is an important step for practical commercial data mining which is often characterised by data sets with far too many variables for model building. In a previous post we looked at all-relevant feature selection using the Boruta package while in this post we consider the same (artificial, toy) examples using the caret package. Max Kuhn kindly listed me as a contributor for some performance enhancements I submitted, but the genius behind the package is all his.

The caret package provides a very flexible framework for the analysis as we shall see, but first we set up the artificial test data set as in the previous article.

## Feature-bc.R - Compare Boruta and caret feature selection
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)
run.name <- "feature-bc"
library("caret")

## Load early to get the warnings out of the way:
library("randomForest")
library("ipred")
library("gbm")

set.seed(1)

## Set up artificial test data for our analysis
n.var <- 20
n.obs <- 200
x <- data.frame(V = matrix(rnorm(n.var*n.obs), n.obs, n.var))
n.dep <- floor(n.var/5)
cat( "Number of dependent variables is", n.dep, "\n")
m <- diag(n.dep:1)

## These are our four test targets
y.1 <- factor( ifelse( x$V.1 >= 0, 'A', 'B' ) )
y.2 <- ifelse( rowSums(as.matrix(x[, 1:n.dep]) %*% m) >= 0, "A", "B" )
y.2 <- factor(y.2)
y.3 <- factor(rowSums(x[, 1:n.dep] >= 0))
y.4 <- factor(rowSums(x[, 1:n.dep] >= 0) %% 2)

The flexibility of the caret package is to a large extent implemented by using control objects. Here we specify to use the randomForest classification algorithm (which is also what Boruta uses) and if the multicore package is available then we use that for extra perfomance (you can also use MPI etc – see the documentation):

control <- rfeControl(functions = rfFuncs, method = "boot", verbose = FALSE,
                      returnResamp = "final", number = 50)

if ( require("multicore", quietly = TRUE, warn.conflicts = FALSE) ) {
    control$workers <- multicore:::detectCores()
    control$computeFunction <- mclapply
    control$computeArgs <- list(mc.preschedule = FALSE, mc.set.seed = FALSE)
}

We will consider from one to six features (using the sizes variable) and then we simply let it lose:

sizes <- 1:6

## Use randomForest for prediction
profile.1 <- rfe(x, y.1, sizes = sizes, rfeControl = control)
cat( "rf     : Profile 1 predictors:", predictors(profile.1), fill = TRUE )
profile.2 <- rfe(x, y.2, sizes = sizes, rfeControl = control)
cat( "rf     : Profile 2 predictors:", predictors(profile.2), fill = TRUE )
profile.3 <- rfe(x, y.3, sizes = sizes, rfeControl = control)
cat( "rf     : Profile 3 predictors:", predictors(profile.3), fill = TRUE )
profile.4 <- rfe(x, y.4, sizes = sizes, rfeControl = control)
cat( "rf     : Profile 4 predictors:", predictors(profile.4), fill = TRUE )

The results are:

rf     : Profile 1 predictors: V.1 V.16 V.6
rf     : Profile 2 predictors: V.1 V.2
rf     : Profile 3 predictors: V.4 V.1 V.2
rf     : Profile 4 predictors: V.10 V.11 V.7

If you recall the feature selection with Boruta article, then the results there were

Profile 1: V.1, (V.16, V.17)
Profile 2: V.1, V.2, V,3, (V.8, V.9, V.4)
Profile 3: V.1, V.4, V.3, V.2, (V.7, V.6)
Profile 4: V.10, (V.11, V.13)

To show the flexibility of caret, we can run the analysis with another of the built-in classifiers:

## Use ipred::ipredbag for prediction
control$functions <- treebagFuncs
profile.1 <- rfe(x, y.1, sizes = sizes, rfeControl = control)
cat( "treebag: Profile 1 predictors:", predictors(profile.1), fill = TRUE )
profile.2 <- rfe(x, y.2, sizes = sizes, rfeControl = control)
cat( "treebag: Profile 2 predictors:", predictors(profile.2), fill = TRUE )
profile.3 <- rfe(x, y.3, sizes = sizes, rfeControl = control)
cat( "treebag: Profile 3 predictors:", predictors(profile.3), fill = TRUE )
profile.4 <- rfe(x, y.4, sizes = sizes, rfeControl = control)
cat( "treebag: Profile 4 predictors:", predictors(profile.4), fill = TRUE )

This gives:

treebag: Profile 1 predictors: V.1 V.16
treebag: Profile 2 predictors: V.2 V.1
treebag: Profile 3 predictors: V.1 V.3 V.2
treebag: Profile 4 predictors: V.10 V.11 V.1 V.7 V.13

And of course, if you have your own favourite model class that is not already implemented, then you can easily do that yourself. We like gbm from the package of the same name, which is kind of silly to use here because it provides variable importance automatically as part of the fitting process, but may still be useful. It needs numeric predictors so we do:

## Use gbm for prediction
y.1 <- as.numeric(y.1)-1
y.2 <- as.numeric(y.2)-1
y.3 <- as.numeric(y.3)-1
y.4 <- as.numeric(y.4)-1

gbmFuncs <- treebagFuncs
gbmFuncs$fit <- function (x, y, first, last, ...) {
    library("gbm")
    n.levels <- length(unique(y))
    if ( n.levels == 2 ) {
        distribution = "bernoulli"
    } else {
        distribution = "gaussian"
    }
    gbm.fit(x, y, distribution = distribution, ...)
}
gbmFuncs$pred <- function (object, x) {
    n.trees <- suppressWarnings(gbm.perf(object,
                                         plot.it = FALSE,
                                         method = "OOB"))
    if ( n.trees <= 0 ) n.trees <- object$n.trees
    predict(object, x, n.trees = n.trees, type = "link")
}
control$functions <- gbmFuncs

n.trees <- 1e2                          # Default value for gbm is 100

profile.1 <- rfe(x, y.1, sizes = sizes, rfeControl = control, verbose = FALSE,
                 n.trees = n.trees)
cat( "gbm    : Profile 1 predictors:", predictors(profile.1), fill = TRUE )
profile.2 <- rfe(x, y.2, sizes = sizes, rfeControl = control, verbose = FALSE,
                 n.trees = n.trees)
cat( "gbm    : Profile 2 predictors:", predictors(profile.2), fill = TRUE )
profile.3 <- rfe(x, y.3, sizes = sizes, rfeControl = control, verbose = FALSE,
                 n.trees = n.trees)
cat( "gbm    : Profile 3 predictors:", predictors(profile.3), fill = TRUE )
profile.4 <- rfe(x, y.4, sizes = sizes, rfeControl = control, verbose = FALSE,
                 n.trees = n.trees)
cat( "gbm    : Profile 4 predictors:", predictors(profile.4), fill = TRUE )

And we get the results below:

gbm    : Profile 1 predictors: V.1 V.10 V.11 V.12 V.13
gbm    : Profile 2 predictors: V.1 V.2
gbm    : Profile 3 predictors: V.4 V.1 V.2 V.3 V.7
gbm    : Profile 4 predictors: V.11 V.10 V.1 V.6 V.7 V.18

It is all good and very flexible, for sure, but I can’t really say it is better than the Boruta approach for these simple examples.

Jump to comments.

Feature selection: All-relevant selection with the Boruta package

2010-11-15T10:04:00Z

Feature selection is an important step for practical commercial data mining which is often characterised by data sets with far too many variables for model building. There are two main approaches to selecting the features (variables) we will use for the analysis: the minimal-optimal feature selection which identifies a small (ideally minimal) set of variables that gives the best possible classification result (for a class of classification models) and the all-relevant feature selection which identifies all variables that are in some circumstances relevant for the classification.

In this article we take a first look at the problem of all-relevant feature selection using the Boruta package by Miron B. Kursa and Witold R. Rudnicki. This package is developed for the R statistical computing and analysis platform.

Background

All-relevant feature selection is extremely useful for commercial data miners. We deploy it when we want to understand the mechanisms behind the behaviour or subject of interest, rather than just building a black-box predictive model. This understanding leads us to a better appreciation of our customers (or other subject under investigation) and not just how, but why they behave as they do, which is useful for all areas of the business, including strategy and product development. More narrowly, it also help us define the variables that we want to observe which is what will really make a difference in our ability to predict behaviour (as opposed to, say, run the data mining application a little longer).

I really like the theoretical approach that the Boruta package tries to implement. It is based on the more general idea that by adding randomness to a system and then collecting results from random samples of the bigger system, one can actually reduce the misleading impact of randomness in the original sample.

For the implementation, the Boruta package relies on a random forest classification algorithm. This provides an intrinsic measure of the importance of each feature, known as the Z score. While this score is not directly a statistical measure of the significance of the feature, we can compare it to random permutations of (a selection of) the variables to test if it is higher than the scores from random variables. This is the essence of the implementation in Boruta.

The tests

This article is a first investigation into the performance of the Boruta package. For this initial examination we will use a test data sample that we can control so we know what is important and what is not. We will consider 200 observations of 20 normally distributed random variables:

run.name <- "feature-1"
library("Boruta")
set.seed(1)
## Set up artificial test data for our analysis
n.var <- 20
n.obs <- 200
x <- data.frame(V=matrix(rnorm(n.var*n.obs), n.obs, n.var))

Normal distribution has the advantage of simplicity, but for commercial application where highly non-normally distributed features like money spent are important may not be the best test. Nevertheless, we will use it for now and define a simple utility function before we get on to the tests:

## Utility function to make plots of Boruta test results
make.plots <- function(b, num,
                       true.var = NA,
                       main = paste("Boruta feature selection for test", num)) {
    write.text <- function(b, true.var) {
        if ( !is.na(true.var) ) {
            text(1, max(attStats(b)$meanZ), pos = 4,
                 labels = paste("True vars are V.1-V.",
                     true.var, sep = ""))        
        }
    }
    plot(b, main = main, las = 3, xlab = "")
    write.text(b, true.var)
    png(paste(run.name, num, "png", sep = "."), width = 8, height = 8,
        units = "cm", res = 300, pointsize = 4)
    plot(b, main = main, lwd = 0.5, las = 3, xlab = "")
    write.text(b, true.var)
    dev.off()
}

Test 1: Simple test of single significant variable

For a simple classification based on a single variable, Boruta performs well: while it identifies three variables as being potentially important, this does include the true variable (V.1) and the plot clearly shows it as being by far the most significant.

## 1. Simple test of single variable
y.1 <- factor( ifelse( x$V.1 >= 0, 'A', 'B' ) )

b.1 <- Boruta(x, y.1, doTrace = 2)
make.plots(b.1, 1)

Figure 1: Simple test of Boruta feature selection with single variable.

Test 2: Simple test of linear combination of variables

With a test of a linear combination of the first four variables where the weights are decreasing from 4 to 1, we begin to get closer to the limitations of the approach.

## 2. Simple test of linear combination
n.dep <- floor(n.var/5)
print(n.dep)

m <- diag(n.dep:1)

y.2 <- ifelse( rowSums(as.matrix(x[, 1:n.dep]) %*% m) >= 0, "A", "B" )
y.2 <- factor(y.2)

b.2 <- Boruta(x, y.2, doTrace = 2)
make.plots(b.2, 2, n.dep)

Figure 2: Simple test of Boruta feature selection with linear combination of four variables.

The implementation correctly identified the first three variables (with weights 4, 3, and 2, respectively) as being important, but it had the fourth variable as possible along with the two random variables V.8 and V.9. Still, six variables are more approachable than twenty.

Test 3: Simple test of less-linear combination of four variables

For this text and the following we consider less obvious combinations of the first four variables. If we just count how many of them are positive, then we get to a situation where Boruta excels (because random forests excel at this type of problem).

## 3. Simple test of less-linear combination
y.3 <- factor(rowSums(x[, 1:n.dep] >= 0))
print(summary(y.3))
b.3 <- Boruta(x, y.3, doTrace = 2)
print(b.3)
make.plots(b.3, 3, n.dep)

Figure 3: Simple test of Boruta feature selection counting the positives of four variables.

Test 4: Simple test of non-linear combination

For a spectacular fail of the Boruta approach we will have to consider a classification in the hyperplane of the four variables. For this simple example, we simply count if there are an even or odd number of positive values among the first four variables:

## 4. Simple test of non-linear combination
y.4 <- factor(rowSums(x[, 1:n.dep] >= 0) %% 2)
b.4 <- Boruta(x, y.4, doTrace = 2)
print(b.4)
make.plots(b.4, 4, n.dep)

Figure 4: Simple test of Boruta feature selection with non-linear combination of four variables

Ouch. The package rejects the four known significant variables. It is too hard for the random forest approach. Increasing the number of observations to 1,000 does not help though at 5,000 observations Boruta identifies the four variables right.

Limitations

Some limitations of the Boruta package are worth highlighting:

It only works with classification (factor) target variables. I am not sure why: as far as I remember, the random forest algorithm also provides a variable significance score when it is used as a predictor, not just when it is run as a classifier.
It does not handle missing (NA) values at all. This is quite a problem when working with real data sets, and a shame as random forests are in principle very good at handling missing values. A simple re-write of the package using the party package instead of randomForest should be able to fix this issue.
It does not seem to be completely stable. I have crashed it on several real-world data sets and am working on a minimal set to send to the authors.

But this is a really promising approach, if somewhat slow on large sets. I will have a look at some real-world data in a future post.

Jump to comments.

Big data for R

2010-08-05T08:22:00Z

Revolutions Analytics recently announced their "big data" solution for R. This is great news and a lovely piece of work by the team at Revolutions.

However, if you want to replicate their analysis in standard R, then you can absolutely do so and we show you how.

Data preparation

First you need to prepare the rather large data set that they use in the Revolutions white paper. The preparation script shown below does two passes over alal the files which is not needed: changing it to a single pass is left as an exercise for the reader.... Note that the following script will take a while to run and will need some 30-odd gig of free disk space (another exercise: get rid of the airlines.csv file), but once it is done the analysis is fast.


#!/usr/bin/Rscript
## big.R - Preprocess the airline data
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)

## Install the packages we will use
install.packages("bigmemory",
                 dependencies = c("Depends", "Suggests", "Enhances"))

## Data sets are downloaded from the Data Expo '09 web site at
## http://stat-computing.org/dataexpo/2009/the-data.html
for (year in 1987:2008) {
    file.name <- paste(year, "csv.bz2", sep = ".")
    if ( !file.exists(file.name) ) {
        url.text <- paste("http://stat-computing.org/dataexpo/2009/",
                          year, ".csv.bz2", sep = "")
        cat("Downloading missing data file ", file.name, "\n", sep = "")
        download.file(url.text, file.name)
    }
}

## Read sample file to get column names and types
d <- read.csv("2008.csv.bz2")
integer.columns <- sapply(d, is.integer)
factor.columns  <- sapply(d, is.factor)
factor.levels   <- lapply(d[, factor.columns], levels)
n.rows <- 0L

## Process each file determining the factor levels
## TODO: Combine with next loop
for (year in 1987:2008) {
    file.name <- paste(year, "csv.bz2", sep = ".")
    cat("Processing ", file.name, "\n", sep = "")
    d <- read.csv(file.name)
    n.rows <- n.rows + NROWS(d)
    new.levels <- lapply(d[, factor.columns], levels)
    for ( i in seq(1, length(factor.levels)) ) {
        factor.levels[[i]] <- c(factor.levels[[i]], new.levels[[i]])
    }
    rm(d)
}
save(integer.columns, factor.columns, factor.levels, file = "factors.RData")

## Now convert all factors to integers so we can create a bigmatrix of the data
col.classes <- rep("integer", length(integer.columns))
col.classes[factor.columns] <- "character"
cols  <- which(factor.columns)
first <- TRUE
csv.file <- "airlines.csv"   # Write combined integer-only data to this file
csv.con  <- file(csv.file, open = "w")

for (year in 1987:2008) {
    file.name <- paste(year, "csv.bz2", sep = ".")
    cat("Processing ", file.name, "\n", sep = "")
    d <- read.csv(file.name, colClasses = col.classes)
    ## Convert the strings to integers
    for ( i in seq(1, length(factor.levels)) ) {
        col <- cols[i]
        d[, col] <- match(d[, col], factor.levels[[i]])
    }
    write.table(d, file = csv.con, sep = ",", 
                row.names = FALSE, col.names = first)
    first <- FALSE
}
close(csv.con)

## Now convert to a big.matrix
library("bigmemory")
backing.file    <- "airlines.bin"
descriptor.file <- "airlines.des"
data <- read.big.matrix(csv.file, header = TRUE,
                        type = "integer",
                        backingfile = backing.file,
                        descriptorfile = descriptor.file,
                        extraCols = c("age"))

Sample analysis

All done now. Sample analysis:


## bigScale.R - Replicate the analysis from http://bit.ly/aTFXeN with normal R
##   http://info.revolutionanalytics.com/bigdata.html
## See big.R for the preprocessing of the data

## Load required libraries
library("biglm")
library("bigmemory")
library("biganalytics")
library("bigtabulate")

## Use parallel processing if available
## (Multicore is for "anything-but-Windows" platforms)
if ( require("multicore") ) {
    library("doMC")
    registerDoMC()
} else {
    warning("Consider registering a multi-core 'foreach' processor.")
}

day.names <- c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday",
               "Saturday", "Sunday")

## Attach to the data
descriptor.file <- "airlines.des"
data <- attach.big.matrix(dget(descriptor.file))

## Replicate Table 5 in the Revolutions document:
## Table 5
t.5 <- bigtabulate(data,
                   ccols = "DayOfWeek",
                   summary.cols = "ArrDelay", summary.na.rm = TRUE)
## Pretty-fy the outout
stat.names <- dimnames(t.5.2$summary[[1]])[2][[1]]
t.5.p <- cbind(matrix(unlist(t.5$summary), byrow = TRUE,
                      nrow = length(t.5$summary),
                      ncol = length(stat.names),
                      dimnames = list(day.names, stat.names)),
               ValidObs = t.5$table)
print(t.5.p)
#             min  max     mean       sd    NAs ValidObs
# Monday    -1410 1879 6.669515 30.17812 385262 18136111
# Tuesday   -1426 2137 5.960421 29.06076 417965 18061938
# Wednesday -1405 2598 7.091502 30.37856 405286 18103222
# Thursday  -1395 2453 8.945047 32.30101 400077 18083800
# Friday    -1437 1808 9.606953 33.07271 384009 18091338
# Saturday  -1280 1942 4.187419 28.29972 298328 15915382
# Sunday    -1295 2461 6.525040 31.11353 296602 17143178

## Figure 1
plot(t.5.p[, "mean"], type = "l", ylab="Average arrival delay")

Just like the Revolutions paper. You can now use biglm.big.matrix and bigglm.big.matrix for basic regression and there are also k-means clustering and other functions.

I must admit here that I do not understand the Revolutions regression example, so I have not attempted to replicate it here. It seems kind of sad if they change the syntax to be incompatible with standard R formulas, which is what appears to be happening.

Credit to Michael Kane and Jay Emerson of Yale who showed much of this in their poster The Airline Data Set... What's the big deal?.

Jump to comments.

Area Plots with Intensity Coloring

2010-07-13T07:47:00Z

I am not sure apeescape’s ggplot2 area plot with intensity colouring is really the best way of presenting the information, but it had me intrigued enough to replicate it using base R graphics.

The key technique is to draw a gradient line which R does not support natively so we have to roll our own code for that. Unfortunately, lines(..., type="l") does not recycle the colour col= argument, so we end up with rather more loops than I thought would be necessary.

(The answer is not to use lines(..., type="h") which, confusingly, does recycle the colour col= argument. This one had me for a while, but the type=h lines always start from zero so you do not get the gradient feature.)

We also get a nice opportunity to use the under-appreciated read.fwf function.

##!/usr/bin/Rscript
## nino.R - another version of http://bit.ly/9P9Gh1
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)

## Get the data from the NOAA server
nino <- read.fwf("http://www.cpc.noaa.gov/data/indices/wksst.for",
                 widths=c(-1, 9, rep(c(-5, 4, 4), 4)),
                 skip=4,
                 col.names=c("Week",
                     paste(rep(c("Nino12","Nino3","Nino34","Nino4"), rep(2, 4)),
                           c("SST", "SSTA"), sep=".")))

## Make the date column something useful
nino$Week <- as.Date(nino$Week, format="%d%b%Y")

## Make colour gradients
ncol <- 50
grad.neg <- hsv(4/6, seq(0, 1, length.out=ncol), 1) # Blue gradient
grad.pos <- hsv(  0, seq(0, 1, length.out=ncol), 1) # Red gradient

## Make plot
plot(Nino34.SSTA ~ Week, data=nino, type="n",
     main="Nino34", xlab="Date", ylab="SSTA", axes=FALSE)
do.call(function (...) rect(..., col="gray85", border=NA),
        as.list(par("usr")[c(1, 3, 2, 4)]))

y <- nino$Nino34.SSTA                   # The values we will plot
x <- nino$Week

axis.Date(1, x=x, tck=1, col="white")
axis(2, tck=1, col="white")
box()

idx <- integer(NROW(nino))
idx[y >= 0] <- 1 + round( y[y >= 0] * (ncol - 1) / max( y[y >= 0]), 0)
idx[y <  0] <- 1 + round(-y[y <  0] * (ncol - 1) / max(-y[y <  0]), 0)

draw.gradient <- function(x, ys, cols) {
    xs <- rep(x, 2)
    for (i in seq(1, length(ys)-1))
        plot.xy(list(x=xs, y=c(ys[i], ys[i+1])), type="l", col=cols[i])
}

for (i in 1:length(x)) {
    ys <- seq(0, y[i], length.out=idx[i]+1)
    cols <- (if (y[i] >=0) grad.pos else grad.neg)
    draw.gradient(x[i], ys, cols)
}

The result is a decent gradient:

I deliberately omitted the scale legend on the right hand side following Allan’s First Law of Happy Graphics: Thou shall not present the same information twice.

For less dense information, you should increase the line width. That is left to the reader. (Hint: it is hard to get just right in base graphics, but lwd <- ceiling(par("pin")[1] / dev.size("in")[1] * dev.size("px")[1] / length(x)) could be a starting point for an approximation. We really need gradient-filled polygons in base R.)

Jump to comments.

Employee productivity as function of number of workers revisited

2010-06-22T11:20:00Z

We have a mild obsession with employee productivity and how that declines as companies get bigger. We have previously found that when you treble the number of workers, you halve their individual productivity which is mildly scary.

Let’s try the FTSE-100 index of leading UK companies to see if they are significantly different from the S&P 500 leading American companies that we analyzed four years ago.

We will of course use the R statistical computing and analysis platform for our analysis, and once again we are grateful to Yahoo Finance for providing the data.

The analysis script is available as ftse100.R and is really simple:

## ftse100.R - Display employee productivity for FTSE-100 consitituents
## Copyright © 2010 Allan Engelhardt <http://www.cybaea.net/>
## All Rights Reserved.

## Get the index constituents.
ftse.100 <- read.csv(file = "http://uk.old.finance.yahoo.com/d/quotes.csv?s=@%5EFTSE&f=s&e=.csv", header = FALSE)
names(ftse.100) <- c("symbol")
data <- data.frame(symbol=NULL, employees=NULL, profit=NULL, sector=NULL)

## For each stock symbol, get employees, profit, and sector
for (symbol in ftse.100$symbol) {
    profile.url <- paste("http://uk.finance.yahoo.com/q/pr?s=", symbol, sep="")
    con <- url(profile.url, open = "r")
    text <- readChar(con, 2^24)     # enough bytes
    close(con)
    x <- sub('.*Number of employees:</td><td.*?>[[:space:]]*([[:digit:],]+).*', "\\1", text, ignore.case = TRUE)
    x <- gsub(',', '', x)
    empl <- tryCatch(as.integer(x), warning = function(x) NA)
    x <- sub('.*Net Profit.*?</td><td.*?>[[:space:]]*([+-]?[[:digit:],]+).*', '\\1', text)
    x <- gsub(',', '', x)
    profit <- tryCatch(as.integer(x)*1e6, warning = function(x) NA)
    sector <- sub('.*Sector:</td><td.*?>(.*?)</td>.*', '\\1', text)
    if (any(c(empl, profit) <= 0, is.na(c(empl, profit)))) {
        cat("Error parsing symbol", symbol, "see", profile.url, "\n")
    } else {
        data <- rbind(data, data.frame(symbol=symbol, employees=empl, profit=profit, sector=sector))
    }
    Sys.sleep(1)
}

## Save the data so we don't have to hit Yahoo all the time.
save(data, file = "data.RData")

## Save plot to file:
#png(filename="ftse100.png", width=800, height=800, pointsize=14, bg="white", res=100)

opar <- par(cex.sub = sqrt(sqrt(2)), font.sub = 3, font.lab = 2)

## x and y coordinates of plot and plot limits
x <- with(data, employees)
y <- with(data, profit/employees)
xlim <- c(10^floor(log10(min(x))), 10^ceiling(log10(max(x))))
ylim <- c(10^floor(log10(min(y))), 10^ceiling(log10(max(y))))

## Set up to display different color and symbols
plot_col <- 1
plot_pch <- 1
markers <- 21:25
pchs <- rep(markers, ceiling(length(levels(data$sector))/length(markers)))
palette(rainbow(length(levels(data$sector)), start=3/6, end=6/6))

# Make empty plot:
plot.new()
plot(profit/employees ~ employees, data = data[FALSE, ], 
     type = "p", pch = pchs[plot_pch], col = plot_col,
     log="xy", xaxp = c(xlim, 1), yaxp = c(ylim, 1), xlim = xlim, ylim = ylim,
     main = "Profit per employee (FTSE 100)", xlab = "Employees", ylab = "Profit per employees (GBP)")

## Plot each sector
for (sector in levels(data$sector)) {
    plot.xy(xy.coords(with(data[data$sector == sector,], employees),
                      with(data[data$sector == sector,], profit/employees),
                      log = "xy", xlab = "", ylab = ""),
            type = "p", pch = pchs[plot_pch], col = plot_col, bg = plot_col)
    plot_pch <- plot_pch + 1
    plot_col <- plot_col + 1
}
legend(x = "bottomleft", legend = levels(data$sector), title = "Industry Sectors", 
       col = palette(), pt.bg = palette(), pch = pchs, cex = 2/3, pt.cex = 1, ncol = 2)

## Fit a linear model to the log-log data:
m <- lm(log10(y) ~ log10(x))
xl <- c(xlim[1]*5, xlim[2]/5)
yl <- 10^predict(m, data.frame(x = xl))
lines(xl, yl, col = "darkred", lty = "dashed", lwd = 2)
t <- sprintf("Power = %0.3g", m$coefficients[2])
text(xl[2], yl[2], t, adj = c(0.25, -1.5), col = "darkred", font = 2)

## All done.
par(opar)
dev.off()

Leave it to run and this is what you get:

The power law still broadly holds. In a large company, the productivity of the individual employee is only ¼ of the productivity in a company with one-tenth of the number of workers.

The analysis for the FTSE All-Share index is easy (ftse-all.R) and gives a slope of -0.7605541 for the 301 companies with the required information, which is much worse. More convincingly, fitting the companies with more than 1,000 employees (to avoid some bias of smaller companies needing to have large profits per employee in order to be big enough to afford a stock market listing) gives a slope of -0.2838.

Jump to comments.

Comparing standard R with Revoutions for performance

2010-06-17T09:05:00Z

Following on from my previous post about improving performance of R by linking with optimized linear algebra libraries, I thought it would be useful to try out the five benchmarks Revolutions Analytics have on their Revolutionary Performance pages.

For convenience I collected their tests into a single script revolution_benchmark.R that I can simply run with Rscript --vanilla revolution_benchmark.R.

The results, compared with the speed-up factors Revolution claims for their version:

Revolutions benchmarks compared with R on x86_64 system
	R	R + ATLAS	Speed-up	Revolution’s claimed speed-up
Matrix Multiply	360.96	9.30	37.8	41.0
Cholesky Factorization	27.28	5.65	3.8	21.0
Singular Value Decomposition	98.73	23.57	3.2	12.6
Principal Components Analysis	454.55	40.92	10.1	15.2
Linear Discriminant Analysis	271.44	79.61	2.4	4.4

In all instances Revolution’s claimed speed-up is greater, though probably not significantly so for the Matrix Multiply test and hardly so for the Principal Components Analysis. (Of course, I do not have a copy of Revolution Analytics’ product, so I can’t verify their claims or make a comparable test.)

Whether saving 48 seconds on a linear discriminant analysis is enough to justify buying the product is a decision I leave to you: you know what analysis you do. For me, there are (many) orders of magnitudes to be gained by better algorithms and better variable selections so I am not too worried about factors of 2 or even 10. For extra raw power, I run R on a cloud service like AWS which scales well for many problems and is easy to do with stock R while I guess there are some sort of license implications if you wanted to do the same with Revolution’s product. (But I like Revolution and am still trying to find an excuse to use their product.)

Your mileage may vary.

Jump to comments.

Faster R through better BLAS

2010-06-15T10:21:00Z

Can we make our analysis using the R statistical computing and analysis platform run faster? Usually the answer is yes, and the best way is to improve your algorithm and variable selection.

But recently David Smith was suggesting that a big benefit of their (commercial) version of R was that it was linked to a to a better linear algebra library. So I decided to investigate.

The quick summary is that it only really makes a difference for fairly artificial benchmark tests. For “normal” work you are unlikely to see a difference most of the time.

The environment

I use R on a 64-bit Fedora 12 Linux system. Fortunately, it is very easy to rebuild R using different libraries on this platform. For the following, I will assume that you have a working rpmbuild environment. The test system has a quad core Intel Xeon E5420 CPU with each core running at 2.50 GHz.

Benchmarks

Benchmarking R is complex. Very complex. But for this simple test we use two tests from the R Benchmarks page: MASS-ex.R and R-benchmark-25.R. The first is a simple benchmark using the examples from the MASS package, and has the advantage that it reflects real-world problems and real-world analysis, albeit small problems and short analysis. The second is a much more artificial example and primarily test matrix operations.

We run the MASS benchmark as:

/usr/bin/time -p R --vanilla CMD BATCH MASS-ex.R /dev/null

While the R-benchmark-25 is simply:

Rscript --vanilla R-benchmark-25.R

For the MASS benchmark we simply capture the real elapsed time while the R benchmark 2.5 provides more detailed output for the three classes of tests (matrix calculation, -functions, and program execution) as well as overall summaries. They are all shown in the table below.

Compiler-optimized R

For the experiments that follow the first thing to do is to grab copies of the source RPMs for R and for ATLAS:

cd ~/rpmbuild/SRPMS
yumdownloader --source atlas R
cd ..

At the time I did this, I got R-2.11.0-1.fc12.src.rpm and atlas-3.8.3-12.fc12.src.rpm. I crank up the level of optimization that I do when building from source so the first thing is to edit ~/.rpmrc to include the line optflags: x86_64 -O3 -march=native -m64 -g. With that in place we can simply do:

rpmbuild --rebuild SRPMS/R-2.11.0-1.fc12.src.rpm  #  Change version numbers as needed
su -c 'rpm -Uhv --force RPMS/x86_64/R*2.11.0-1*.rpm RPMS/x86_64/libRmath*2.11.0-1*.rpm'

We now have a compiler-optimized version of R and we can re-run our tests. It doesn't make much difference, but that is also good to know.

ATLAS BLAS libraries

Now let's try linking to the ATLAS BLAS libraries instead. I assume you have them installed (yum install atlas if not) so you can just grab a copy of R-atlas.diff to change the spec file like this:

rpm -ihv SRPMS/R-2.11.0-1.fc12.src.rpm   # Install to your rpmbuild environment
cd SPECS
wget http://static.cybaea.net/files/R-atlas.diff
patch -o R-atlas.spec R.spec R-atlas.diff
cd ..
rpmbuild -bb SPECS/R-atlas.spec
su -c 'rpm -Uhv --force RPMS/x86_64/R*2.11.0-1*.rpm RPMS/x86_64/libRmath*2.11.0-1*.rpm'

You now have a version of R that uses the ATLAS BLAS libraries, so you can re-run the tests. The results are in the table below in the “Optimized R + Standard ATLAS” row.

As expected, the matrix operations from the R-benchmark-25.R runs a lot faster: they complete in about 30-40% of the time, much of which comes from the multi-threading so all four CPU cores are used.

However, for the analysis-heavy code is MASS-ex.R there is little difference. If anything, we see a tiny increase in running time.

Multi-threaded BLAS libraries make no significant difference to real-world analysis problems using R.

Other BLAS libraries

For good measure we also try an optimized version of ATLAS, but it does not make much difference on the x86_64 architecture:

rpmbuild -D "enable_native_atlas 1" --rebuild SRPMS/atlas-3.8.3-12.fc12.src.rpm
su -c 'rpm -Uhv --force RPMS/x86_64/atlas*3.8.3-12*.rpm'

And (only) for completeness, we also try the standard Netlib BLAS and LAPACK libraries (yum install blas lapack) by the same method as the ATLAS library above but with a slightly different change to the SPEC file: R-blas.diff. It performs a little better than vanilla R.

For more information about rebuilding R with different BLAS libraries, see the linear algebra section in the R Installation and Administration manual.

Benchmark results

Benchmark results for various optimizations of R and the BLAS library
R version	MASS-ex.R		R benchmark 2.5
	Real		Total time		Overall mean		Ⅰ. Matrix calc.		Ⅱ. Matrix functions		Ⅲ. Program.
	secs	index	secs	index	secs	index	secs	index	secs	index	secs	index
Base install	19.00	1.00	78.49	1.00	2.11	1.00	2.32	1.00	3.86	1.00	1.05	1.00
Optimized R	18.98	1.00	76.11	0.97	2.02	0.96	2.36	1.02	3.46	0.90	1.02	0.97
Optimized R + Netlib BLAS	18.56	0.98	73.22	0.93	1.81	0.86	2.36	1.02	2.41	0.62	1.04	0.99
Optimized R + Standard ATLAS	19.43	1.02	16.74	0.21	0.97	0.46	0.90	0.39	1.04	0.27	0.99	0.95
Optimized R + Optimized ATLAS	19.31	1.02	16.36	0.21	0.95	0.45	0.84	0.36	1.02	0.26	1.00	0.95

Jump to comments.

R: Eliminating observed values with zero variance

2010-03-08T14:46:00Z

I needed a fast way of eliminating observed values with zero variance from large data sets using the R statistical computing and analysis platform. In other words, I want to find the columns in a data frame that has zero variance. And as fast as possible, because my data sets are large, many, and changing fast. The final result surprised me a little.

I use the KDD Cup 2009 data sets as my reference for this experiment. (You will need to register to download the data.) It is a realistic example of the type of customer data that I usually work with. It has 50,000 observations of 15,000 variables. To load it into R you'll need a reasonably beefy machine. My workstation has 16GB of memory; if yours have less then use a sample of the data.

We load the data into R and propose a few ways in which we may identify the columns we need:

#!/usr/bin/Rscript
## zero-var.R - find the fastest way of eliminating observations with zero variance
## © 2010 Allan Engelhardt, http://www.cybaea.net

## Read the data file.
## We have already converted it to R format and saved it, so we can do
load("train.RData")
## instead of something like
# train <- read.delim(file="../orange_large_train.data.bz2")

## Some suggestions for zero variance functions:
zv.1 <- function(x) {
    ## The literal approach
    y <- var(x, na.rm = TRUE)
    return(is.na(y) || y == 0)
}
zv.2 <- function(x) {
    ## As before, but avoiding direct comparison with zero
    y <- var(x, na.rm = TRUE)
    return(is.na(y) || y < .Machine$double.eps ^ 0.5)
}
zv.3 <- function(x) {
    ## Maybe it is faster to check for equality than to compute?
    y <- x[!is.na(x)]
    return(all(y == y[1]))
}
zv.4 <- function(x) {
    ## Taking out the special case may speed things up?
    ## (At least for this data set where this case is common.)
    z <- is.na(x)
    if ( all(z) ) return(TRUE);
    y <- x[!z]
    return(all(y == y[1]))
}

Now we just have to load the very useful rbenchmark package and let the machine figure it out:

library("rbenchmark")

cat("Running benchmarks:\n")
benchmark(
          zv1 = { sapply(train, zv.1) },
          zv2 = { sapply(train, zv.2) },
          zv3 = { sapply(train, zv.3) },
          zv4 = { sapply(train, zv.4) },
          replications = 5,
          columns = c("test", "elapsed", "relative", "sys.self"),
          order = "elapsed"
          )

The answer (on my machine) is that it is faster to calculate than to check for equality:

Running benchmarks:
  test elapsed relative sys.self
1  zv1  78.619 1.000000    6.395
2  zv2  79.276 1.008357    6.586
3  zv3 113.024 1.437617    1.735
4  zv4 118.579 1.508274    1.716

The two functions based on the core variance function are easily the fastest (despite having to do arithmetic) while taking out the special case in the equality functions is a Bad Idea.

Can you think of an even faster way to do it?

Jump to comments.

Beautiful Data

2009-07-27T19:38:00Z

O'Reilly's recent publication Beautiful Data has a chapter by Jeff Jonas which is enough reason in itself for me to recommend it. The chapter, Data Finds Data, is also available as a PDF download.

I met Jeff a couple of year ago at an ETech conference, and he is easily one of the smartest people I have ever met who is thinking about data.

Jump to comments.

Massively parallel database for analytics

2009-07-22T13:37:00Z

This is by far the best description of why traditional parallel databases (like Teradata, Greenplum et al.) is a evolutionary dead end. But much more than a theoretical discussion, they have built a solution which they call HadoopDB. It is based on Hadoop, PostgreSQL, and Hive and is completely Open Source. Alternative, column-based, backends to PostgreSQL are being implemented now. Read: Announcing release of HadoopDB.

The Knapsack Problem

2009-07-10T20:30:00Z

David posts a question about how to solve this knapsack problem using the R statistical computing and analysis platform. My reply in the comments seems to have disappeared for a while so here is my proposed solution. See David’s blog for my earlier proposed solution with a very common error.


## http://blog.revolution-computing.com/2009/07/because-its-friday-the-knapsack-problem.html
appetizer.solution <- local (
function (target) {
  app <- c(2.15, 2.75, 3.35, 3.55, 4.20, 5.80)
  r <- 2L
  repeat {
	c <- gtools::combinations(length(app), r=r, v=app, repeats.allowed=TRUE)
	s <- rowSums(c)
	if ( all(s > target) ) {
	  print("No solution found")
	  break
	}
	x <- which( abs(s-target) < 1e-4 )
	if ( length(x) > 0L ) {
	  cat("Solution found: ", c[x,], "\n")
	  break
	}
	r <- r + 1L
  }
})

appetizer.solution(15.05)
# Solution found:  2.15 3.55 3.55 5.8

Brute force works, it just doesn’t scale well. (Note that 7×2.15 is another solution.)

Jump to comments.

OECD Statistics

2009-07-02T20:33:00Z

I am a sucker for good quality data. I wrote about data.gov, the US Government data site before, and now I find OECD Statistics which has some 300 data sets, many of which seems to be readily accessible (though some may require subscription)

Exports in multiple formats, including Excel, CSV, and SDMX.

R tips: Determine if function is called from specific package

2009-06-16T10:27:00Z

I like the "multicore" library for a particular task. I can easily write a combination of if(require("multicore",...)) that means that my function will automatically use the parallel mclapply() instead of lapply() where it is available. Which is grand 99% of the time, except when my function is called from mclapply() (or one of the lower level functions) in which case much CPU trashing and grinding of teeth will result.

So, I needed a function to determine if my function was called from any function in the "multicore" library. Here it is.

First define a generally useful function:


is.in.namespace <-
function (ns) {
  for ( frame in seq(1, sys.nframe(), 1) ) {
	fun <- sys.function(frame);
	env <- environment(fun)
	n   <- environmentName(env)
	if ( n == ns ) return(TRUE);
  }
  return(FALSE);
}

Then we use it for our purpose:


is.in.multicore <- function (...) { return(is.in.namespace("multicore")) }
library("multicore")
stopifnot( mclapply(as.list(1), is.in.multicore)[[1]] == TRUE )
stopifnot(   lapply(as.list(1), is.in.multicore)[[1]] == FALSE )
stopifnot( local( {mclapply <- function(x) return(x); mclapply(is.in.multicore())} ) == FALSE )

Easy when you know how.

Jump to comments.

R tips: Installing Rmpi on Fedora Linux

2009-06-12T10:23:00Z

Somebody on the R-help mailing list asked how to get Rmpi working on his Fedora Linux machine so he could do high-performance computing on a cluster of machines (or a single multicore machine) using the R statistical computing and analysis platform. Since it is unusually painful to get working, I might as well copy the instructions here.

1. Install Open MPI on Fedora Core

First install the openmpi libraries using:

yum install openmpi openmpi-devel openmpi-libs

The default installation on Fedora still doesn’t quite work, so you need to execute the following command as root (only once is required, after installation of the package):

ldconfig /usr/lib64/openmpi/lib/

You are not quite done: for R to work right with the libraries, you need to modify the LD_LIBRARY_PATH environment variable to include the path to the Open MPI libraries. I have the following in my ~/.bash_profile:

export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}${LD_LIBRARY_PATH:+:}/usr/lib64/openmpi/lib/"

Edit your file to contain the same, and execute that line at the command prompt and you are ready to continue.

2. Install the `Rmpi` package for `R`

Now that your Open MPI libraries are set up, and what you do next depends on what version of Rmpi you are installing. Most likely you are installing the latest version in which case the following section applies. The instructions for older versions are retained in a later section for reference.

2.1. Current versions of the `Rmpi` package

Make sure you have executed the ldconfig command and set the LD_LIBRARY_PATH environment variables as described in the previous section before you continue.

Since at least version 0.5-8 of the Rmpi library you can install it from the R command line after you have fixed the Open MPI install. At the R prompt do:

install.packages("Rmpi",
                 configure.args =
                 c("--with-Rmpi-include=/usr/include/openmpi-x86_64/",
                   "--with-Rmpi-libpath=/usr/lib64/openmpi/lib/",
                   "--with-Rmpi-type=OPENMPI"))

It should work and install OK. This is obviously quite a mouthful to remember, but help is at hand through the options() mechanism in R. In your ~/.Rprofile you can add something like:

local({
    my.configure.args <-
        list("Rmpi" =
             c("--with-Rmpi-include=/usr/include/openmpi-x86_64/",
               "--with-Rmpi-libpath=/usr/lib64/openmpi/lib/",
               "--with-Rmpi-type=OPENMPI"),
             ## Not needed for Rmpi but shown to illustrate the format
             "ncdf" =
             c("-with-netcdf_incdir=/usr/include/netcdf",
               "-with-netcdf_libdir=/usr/lib64/")
             );
    options("configure.args" = my.configure.args)
})

Then you can just type install.packages("Rmpi") at the R command prompt to install the package.

2.2. Older versions of the `Rmpi` package

The problem is the configuration file configure.ac which is, unfortunately, completely brain-damaged with hard-coded assumptions about which subdirectories should contain header and library files and no way of overriding it.

Download the latest Rmpi package from CRAN and unpack it using tar zxvf Rmpi_0.5-7.tar.gz. Go to the new Rmpi directory and replace the file configure.ac with the one below (for a x86_64 system; for 32 bit you probably need to change -64 to -32):

 Process this file with autoconf to produce a configure script.

AC_INIT(DESCRIPTION)

AC_PROG_CC

MPI_LIBS=`pkg-config --libs openmpi-1.3.1-gcc-64`
MPI_INCLUDE=`pkg-config --cflags openmpi-1.3.1-gcc-64`
MPITYPE="OPENMPI"
MPI_DEPS="-DMPI2"

AC_CHECK_LIB(util, openpty, [ MPI_LIBS="$MPI_LIBS -lutil" ])
AC_CHECK_LIB(pthread, main, [ MPI_LIBS="$MPI_LIBS -lpthread" ])

PKG_LIBS="${MPI_LIBS} -fPIC"
PKG_CPPFLAGS="${MPI_INCLUDE} ${MPI_DEPS} -D${MPITYPE} -fPIC"

AC_SUBST(PKG_LIBS)
AC_SUBST(PKG_CPPFLAGS)
AC_SUBST(DEFS)

AC_OUTPUT(src/Makevars)

The number 1.3.1 may change in future releases of Fedora: see /usr/lib64/pkgconfig/openmpi-*.pc for the current value.

Still in the Rmpi directory do the following in your shell:

autoconf
cd ..
tar zcvf Rmpi_0.5-7-F11.tar.gz Rmpi
R CMD INSTALL Rmpi_0.5-7-F11.tar.gz

3. Test it

Now Rmpi should be working in R:

> library("Rmpi")
> mpi.spawn.Rslaves(nslaves=2)
    2 slaves are spawned successfully. 0 failed.
master (rank 0, comm 1) of size 3 is running on: server
slave1 (rank 1, comm 1) of size 3 is running on: server
slave2 (rank 2, comm 1) of size 3 is running on: server
> x <- c(10,20)
> mpi.apply(x,runif)
[[1]]
 [1] 0.25142616 0.93505554 0.03162852 0.71783194 0.35916139 0.85082154
 [7] 0.35404191 0.14221315 0.60063773 0.71805190

[[2]]
 [1] 0.84157864 0.63481773 0.38217188 0.67839089 0.27827728 0.35429266
 [7] 0.04898744 0.96601584 0.25687905 0.77381186 0.69011927 0.37391028
[13] 0.19017369 0.51196594 0.51970563 0.15791524 0.21358237 0.69642478
[19] 0.12690207 0.44177656

Jump to comments.

Data Mashups in R from O'Reilly

2009-06-09T11:23:00Z

O’Reilly has published Data Mashups in R as a $4.99 PDF download in their Short Cut series. In 27 pages it takes you through an example of how to combine foreclosure information with maps and geographical information to produce plots like the one below. This is all done with the R statistical computing and analysis platform.

They show how to:

Use regular expressions to parse HTML files
Use the XML package to parse XML data from a web service (Yahoo! Geocode)
Find ERSI shape files for your maps
Use PBSmapping to process and display geographical data (GIS)
Importing and using US Census data with your maps

How to win the KDD Cup Challenge with R and gbm

2009-06-01T07:07:00Z

Hugh Miller, the team leader of the winner of the KDD Cup 2009 Slow Challenge (which we wrote about recently) kindly provides more information about how to win this public challenge using the R statistical computing and analysis platform on a laptop (!).

As a reminder of what we wrote before, the challenge provided two anonymized data set each of 50,000 mobile teleco customers and each entry having 15,000 variables. The task was to find the best churn, up-, and cross-sell models.

Hugh summarizes his team’s approach:

Feature selection was an important first step [we mentioned before that this is key for all successful data mining projects – AE]. We looked at how effective each individual variable was as a predictor, which also allowed us to reading parts of the data only, as the whole dataset didn’t fit in memory [my emphasis – AE]. The assessment here was homebrew, making a simple predictor on half the data and measuring performance (by the AUC measure) on the other half:

For categorical variables we just took the average number of 1's in the response for each category and used this as a predictor

For continuous variables we split the variable up into "bins", as you would a histogram, and again took the average number of 1's in the response for each bin as the predictor.

From this we came up with a set of about 200 variables for each model, which we continued to tinker with. The main model was a gradient boosted machine which used the "gbm" package in R. This basically fits a series of small decision trees, up-weighting the observations that are predicted poorly at each iteration. We used Bernoulli loss and also up-weighted the "1" response class. A fair amount of time was spent optimising the number of trees, how big they should be etc, but a fit of 5,000 trees only took a bit over an hour to fit. The package itself is quite powerful as it gives some useful diagnostics such as relative variable importance, allowing us to exclude some and include others.
We used trees to avoid doing much data cleaning – they automatically allow for extreme results, non-linearity, missing values and handle both categorical and continuous variables. The main adjustment we had to make was to aggregate the smaller categories in the categorical variables, as they tended to distort the fits.

They did this on standard Windows laptops (Intel Core 2 Duo 2.66GHz processor, 2GB RAM, 120Gb hard drive) against a competition that had more computing clusters available than Imelda Marcos had shoes. It is not what you’ve got, it’s how you use it :-).

Congratulations to Hugh and his team!

Jump to comments.

R used by KDD 2009 cup winner of slow challenge

2009-05-31T13:17:00Z

The results from the KDD Cup 2009 challenge (which we wrote about before) are in, and the winner of the slow challenge used the R statistical computing and analysis platform for their winning submission.

The write up (username/password may be required) from Hugh Miller and team at the University of Melbourne includes these points:

Decision tree, stub, or Random Forest as base classifiers with Logistic loss or cross-entropy loss function
Models fit in an hour or so
Used the R statistical package
Most of models run on Windows laptop with Intel Core 2 Duo 2.66GHz processor, 2GB RAM, 120Gb hard drive.

Impressive hardware selection! Well done R. Weka was another popular tool among the top entrants. Key for all of them were clever data preparation and variable substitution. The fast track winners from IBM document this in some detail:

We normalized the numerical variables by range, keeping the sparsity. For the categorical variables, we coded them using at most 11 binary columns for each variable. For each categorical variable, we generated a binary feature for each of the ten most common values, encoding whether the instance had this value or not. The eleventh column encoded whether the instance had a value that was not among the top ten most common values. We removed constant attributes, as well as duplicate attributes.
We replaced the missing values by mean for numerical attributes, and coded them as a separate value for discrete attributes. We also added a separate column for each numeric attribute with missing values, indicating wether the value was missing or not. We also tried another approach for imputing missing values based on KNN.
On the large data set we discretized the 100 numerical variables that had the highest mutual information with the target into 10 bins, and added them as extra features.
We tried PCA on the large data set, but it did not seem to help.
Because we noticed that some of the most predictive attributes were not linearly correlated with the targets, we build shallow decision trees (2-4 levels deep) using single numerical attributes and used their predictions as extra features. We also build shallow decision trees using two features at a time and used their prediction as an extra feature in the hope of capturing some non-additive interactions among features.

Jump to comments.

R tips: Use read.table instead of strsplit to split a text column into multiple columns

2009-05-29T10:53:00Z

Someone on the R-help mailing list had a data frame with a column containing IP addresses in quad-dot format (e.g. 1.10.100.200). He wanted to sort by this column and I proposed a solution involving strsplit. But Peter Dalgaard comes up with a much nicer method using read.table on a textConnection object:

> a <- data.frame(cbind(color=c("yellow","red","blue","red"),
                        status=c("no","yes","yes","no"),
                        ip=c("162.131.58.26","2.131.58.16","2.2.58.10","162.131.58.17")))
> con <- textConnection(as.character(a$ip))
> o <- do.call(order,read.table(con, sep="."))
> close(con)
> a[o,]
   color status            ip
3   blue    yes     2.2.58.10
2    red    yes   2.131.58.16
4    red     no 162.131.58.17
1 yellow     no 162.131.58.26

That is very, very neat! Thank you Peter.

Data.gov

2009-05-22T02:23:00Z

I am always on the lookout for useful data sources for training in statistics, so I am excited that Data.gov has opened for business. The purpose of Data.gov is to increase public access to high value, machine readable datasets generated by the US Government.

This is a great initiative which I look forward to explore when I am not in a tiny airport at 3 am (but hey: they have free wifi) and which I hope other countries will take up.

Are there other catalogues of data sets that you use?

Jump to comments.

SNA with R: Loading large networks using the igraph library

2009-05-06T15:33:00Z

We are interested in Social Network Analysis using the statistical analysis and computing platform R. The documentation for R is voluminous but typically not very good, so this entry is part of a series where we document what we learn as we explore the tool and the packages.

In our previous post on SNA we gave up on using the statnet package because it was not able to handle our data volumes. In this entry we have better success with the igraph package.

The task we are considering is still how to load the network data into the R package’s internal representation. We will assume that the raw data for our analysis is in a transactional format that is typical at least in the Telecommunications and Finance industries. In the former the terminology is Call Detail Record (CDR) and an extract may look a little like the following:


          src,         dest,     start,  duration,type,...
+447000000005,+447000000006,1238510028,        52,call,...
+447000000006,+447000000009,1238510627,       154,call,...
+447000000009,+447000000007,1238511103,        48,call,...
+447000000006,+447000000005,1238511145,        49,call,...
+447000000006,+447000000005,1238511678,        12,call,...
+447000000001,+447000000006,1238511735,       147,call,...
+447000000007,+447000000009,1238511806,        26,call,...
+447000000000,+447000000008,1238511825,        19,call,...
+447000000009,+447000000008,1238511900,        28,call,...
...

Here a record indicates that the customer identified as src called (type=call) the customer dest at the given time start and the call lasted duration seconds. In general, there will be (many) more attributes describing the transaction which are represented by the .... In a Financial Services example, the records may be money transfers between accounts.

Loading the data in the `igraph` package

We are able to load the previous test data with 51 million records easily:

> library("igraph")
> m <- matrix(scan(bzfile("cdr.51M.csv.bz2", open="r"), 
+                  what=integer(0), skip=1, sep=','), 
+             ncol=4, byrow=TRUE)
Read 205266564 items
> ### Vertices are numbered from zero in the igraph library
> m[,1] <- m[,1]-1; m[,2] <- m[,2]-1
> g <- graph.edgelist(m[,c(2,1)])
> E(g)$start    <- as.POSIXct(m[,3], origin="1970-01-01", tz="UTC")
> E(g)$duration <- m[,4]
> ns <- neighborhood.size(g, 1)

Time to up the ante! We have a file with simulated call data records containing over 700 million entries where we suspect the algorithm used is under-estimating nodes with small connections. Let’s check on the first ½ billion records (which seems to more-or-less fit in our available memory on this workstation):

> library("igraph")
### Note that R can only handle 2^31-1 elements in a vector (on any
### platform, including 64-bit), so we need to read this file as a
### list.
> s <- scan("cdr.1e6x1e1.csv", what=rep(list(integer(0)),4), skip=1, sep=',', multi.line=FALSE)
Read 700466826 records
> m <- as.vector(rbind(s[[2]], s[[1]]))
> print(length(m))
[1] 1400933652
> length(m) <- 1e9
> g <- graph(m, directed=TRUE)
> ns <- neighborhood.size(g, 1)
> summary(ns)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00   35.00   40.00   42.92   47.00  101.00 
> hist(ns, xlab="Neighborhood size", main="Distribution of neighborhood size", 
       sub="From cdr.1e6x1e1.1e9")

As we suspected, the Monte Carlo algorithm does not provide enough customers with low calling circle sizes. Fortunately it is very easy to add these separately: the hard part is modelling the larger calling circles. A mix of these two algorithms provide a reasonably good fit to actual customer behaviour. (The cut-off at 100 is a parameter to our Monte Carlo simulation program which indeed was 100 for this run.)

Problems

However, it is not all perfect. When we attempt to add the edge parameters in the obvious way it fails:

> length(s[[3]]) <- 0.5e9
> length(s[[4]]) <- 0.5e9
> E(g)$start     <- s[[3]]
Error: cannot allocate vector of size 3.7 Gb
Execution halted
> E(g)$duration  <- s[[4]]

So we are just at the limit. Probably 100 million records is OK in this environment. But the core igraph library is accessible from C so better performance can probably be achieved this way and certainly pointers are 8 byte structures on this machine so we should not have the silly limits that R imposes on us.

Jump to comments.

CYBAEA Data and Analysis

R code for Chapter 2 of Non-Life Insurance Pricing with GLM

Example 2.5: Moped insurance continued

Frequency model

Severity model

Combining the models

2.4 Case Study: Motorcycle Insurance

The Data

Obtaining the data

Adding meta-data information

Rating factors

Save the data

Problem 1: Aggregate to cells of current tariff

Problem 2: Determine how the duration and number of claims is distributed

1. Number of claims (antskad)

2. Duration

Problem 3: Determine the relativities for claim frequency and severity

Problem 4: Discussions

You may also like these posts:

R code for Chapter 1 of Non-Life Insurance Pricing with GLM

Example 1.2

Example 1.3

Exercise 1.3

You may also like these posts:

doSMP pulled

R versus SAS/SPSS in corporations

Background

The advantages of SAS/SPSS in a commercial environment

1. You can buy the tool for money.

2. You can recruit for the commercial tool

3. You can recruit for the commercial tool.

We can help

You may also like these posts:

Friday quote: what is the question to which this number is the answer?

A warning on the R save format

How to lose your data with save()

What can you do?

Recommendations

You may also like these posts:

Friday quote: the handmaiden and the whore

Spreadsheet errors

You may also like these posts:

Getting started with the Heritage Health Price competition

Prerequisites

Data preparation

Scoring

The simplest benchmarks

Simple single-variable linear models

You may also like these posts:

Benchmarking feature selection with Boruta and caret

The Boruta package

The caret package

You may also like these posts:

Feature selection: Using the caret package

You may also like these posts:

Feature selection: All-relevant selection with the Boruta package

Background

The tests

Test 1: Simple test of single significant variable

Test 2: Simple test of linear combination of variables

Test 3: Simple test of less-linear combination of four variables

Test 4: Simple test of non-linear combination

Limitations

You may also like these posts:

Big data for R

Data preparation

Sample analysis

You may also like these posts:

Area Plots with Intensity Coloring

You may also like these posts:

Employee productivity as function of number of workers revisited

You may also like these posts:

Comparing standard R with Revoutions for performance

You may also like these posts:

Faster R through better BLAS

The environment

Benchmarks

Compiler-optimized R

ATLAS BLAS libraries

Other BLAS libraries

How to lose your data with `save()`

2. Install the `Rmpi` package for `R`

2.1. Current versions of the `Rmpi` package

2.2. Older versions of the `Rmpi` package

Loading the data in the `igraph` package