CYBAEA Data and Analysis

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.

SNA with R: Loading your network data

2009-04-01T16:08:00Z

We are interested in Social Network Analysis using the statistical analysis and computing platform R. As usual with R, the documentation is pretty bad, so this series collects our notes as we learn more about the available packages and how they work. We use here the statnet group of packages, which seems to be the most comprehensive and most actively maintained network analysis packages.

The first task which we consider in this post is to load our data into a network object, which is how all the statnet packages represent a network. Typically for R, the documentation is voluminous but not always as helpful as one could want.

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,...
...

Implementation in the `network` class

In the naive implementation of this data as a network, we would have the sources and destinations (broadly speaking: people) as vertices and the calls as edges. That broadly seems to make sense: people are connected by the calls they make, and that is the social relationship we wish to model.

In the terminology of the network class, that means that our network will be directed (calls and money transfers have a direction from one person to another) and will need to allow multiple edges between the same endpoints (because any one person can, and indeed usually will, make several calls to the same other person).

We could consider dropping the multiple attribute of the network and instead represent the fact that A has called B with a single edge and perhaps have the number of calls and their total duration as an edge attribute. We will investigate this another time, but it is surely a less faithful representation of the data that we have (and we would need to drop information like the time of call).

Mapping customer identifiers to network vertex numbers

One thing they seem to forget to tell you in the documentation is that when you import your data your vertex identifiers (which in our case is customer or account numbers) must be changed to number the vertices and that this numbering must be sequential and start from 1. Being used to an environment where the vertex identifiers are arbitrary (and arrays usually start from 0), this one had me puzzled for a while. The error message that tells you your vertex numbering is not what the package expected is spectacularly unhelpful:

> n <- network(m, matrix.type="edgelist", directed=TRUE, multiple=TRUE)
Error in add.edges(g, as.list(x[, 1]), as.list(x[, 2]), edge.check = edge.check) : 
  (edge check) Illegal vertex reference in addEdges_R.  Exiting.

For the discussion that follows, we will assume that you have changed your identifies externally to R.

Loading the data

The good news is that our data is essentially in a format that the network package calls edge list and which it can import directly.

I say “essentially” because for some strange reason the package expects the destination to come before the source which seems ass-backwards to me. But assume we have our data in a file cdr.csv like this (we only have calls here):

       src,      dest,     start,  duration
         5,         6,1238510028,        52
         6,         9,1238510627,       154
         9,         7,1238511103,        48
         6,         5,1238511145,        49
...

Then we can load the data into R easily:

> library("network")
> m <- matrix(scan(file="cdr.csv", what=integer(0), skip=1, sep=','), ncol=4, byrow=TRUE)
Read 1896 items
> # Swap columns for ass-backward network package
> m[,c(1,2)] <- m[,c(2,1)]

> # Create network
> net <- network(m, matrix.type="edgelist", directed=TRUE, multiple=TRUE)

> summary(net)
Network attributes:
 vertices = 10
 directed = TRUE
 hyper = FALSE
 loops = FALSE
 multiple = TRUE
 bipartite = FALSE
 total edges = 474 
   missing edges = 0 
   non-missing edges = 474 
 density = 5.266667 

Vertex attributes:
 vertex.names:
   character valued attribute
   10 valid vertex names

No edge attributes

Network adjacency matrix:
Error in as.matrix.network.adjacency(x = x, attrname = attrname, ...) : 
  Multigraphs not currently supported in as.matrix.network.adjacency.  Exiting.
In addition: Warning message:
In network.density(x) :
  Network is multiplex - no general way to define density.  Returning value for a non-multiplex network (hope that's what you wanted).

OK, that's a lot of warnings, but it basically worked. We have figured out how to load our network data into the network package in R.

Performance

We can’t do an exhaustive performance review now, but let us at least make sure we can load medium-sized networks. We change our CDR simulator to emit the desitnation before the source just like network likes it and let it run.

The first file has 2,645,288 (simulated) CDR lines from 100k customers and it loads OK on our small development workstation even with the default settings:

> library("network")
> n <- network(matrix(scan(file="cdr.1e5x1e0.csv", 
                           what=integer(0), skip=1, sep=','), 
                      ncol=4, byrow=TRUE), 
               matrix.type="edgelist", directed=TRUE, multiple=TRUE)
Read 10581152 items
> proc.time()
   user  system elapsed 
138.304   1.597 140.878 
> save(n, file="n.RData", ascii=FALSE, compress=FALSE)

The size of the saved network object is 373MB (only 27MB compressed).

We can potentially save some time and memory by not explicitly not performing the edge check (again: the documentation frustrates us and is silent on what the defaults are for the network call we used above) so we try this for our next file with 51,316,641 lines of CDR data (again for 100k customers) which also saves us some column swapping:

> library("network")
> m <- matrix(scan(file="cdr.51M.csv", 
                   what=integer(0), skip=1, sep=','),
              ncol=4, byrow=TRUE)
Read 205266564 items
> num_vert <- max(m[,1], m[,2])
> num_vert
[1] 100000
> n <- network.initialize(n=num_vert, directed=TRUE, multiple=TRUE)
> add.edges(n, tail=m[,2], head=m[,1], edge.check=FALSE)
> proc.time()
(several hours: I’ll let you know when it is done)
> rm(m)
> save(n, file="n.RData", ascii=FALSE, compress=TRUE)

Our attempted optimization did not seem to matter and this network is too big for the machine and the network package. Building the network was painful as I was working on the workstation at the same time. The machine has 16GB installed RAM, but it was clearly running out and swapping extensively.

51 million might be a reasonable size data set for some Financial Services applications but it is clearly a trivial number of records for Telecommunications. I’ll need to do some more digging around.

Does anybody have any SNA benchmarks? I like the KXEN implementation for its simplicity and speed so I might get a copy and try it out. Any R performance experts who could make suggestions in the comments? How big are your networks?

Jump to comments.

R tips: Swapping columns in a matrix

2009-03-31T15:59:00Z

Using R, the statistical analysis and computing platform, swapping two columns in a matrix is really easy: m[ , c(1,2)] <- m[ , c(2,1)].

Note, however, that this does not swap the column names (if you have any) but only the values. You could do something like colnames(m)[c(1,2)] <- colnames(m)[c(2,1)] if you need the names changed as well, but better is perhaps just to assign:

m <- m[ , c(2, 1, 3:ncol(m)) ]

R tips: Eliminating the “save workspace image” prompt on exit

2009-03-26T08:14:00Z

When using R, the statistical analysis and computing platform, I find it really annoying that it always prompts to save the workspace when I exit. This is how I turn it off.

I wish there was an option to change the default of the q/quit functions. I start and stop R frequently and so the exit question which I have to answer every time is really annoying:

Save workspace image? [y/n/c]:

Why is there no R option to disable this prompt? If I want to save the image, I have already saved it. And I don’t like the default name anyhow, preferring to give my own with save.image(file=...). For a while, I had a function defined in my ~/.Rprofile that terminated the session without prompting.

exit <- function() { q("no") }

While this means I can type exit() and avoid the annoying prompt, in practice I normally type Control-D to end the session which still calls the normal q function with its annoying prompt.

So instead I use the alias functionality of my (bash) shell to change the default. In my ~/.bashrc I now have

alias R="$(/usr/bin/which R) --no-save"

And finally I am happy. But I still think R should have an option (accessible through options) to change the default behavior.

Jump to comments.

CYBAEA Data and Analysis

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

Benchmark results

You may also like these posts:

R: Eliminating observed values with zero variance

You may also like these posts:

Beautiful Data

You may also like these posts:

Massively parallel database for analytics

The Knapsack Problem

You may also like these posts:

OECD Statistics

R tips: Determine if function is called from specific package

You may also like these posts:

R tips: Installing Rmpi on Fedora Linux

1. Install Open MPI on Fedora Core

2. Install the Rmpi package for R

2.1. Current versions of the Rmpi package

2.2. Older versions of the Rmpi package

3. Test it

You may also like these posts:

Data Mashups in R from O'Reilly

How to win the KDD Cup Challenge with R and gbm

You may also like these posts:

R used by KDD 2009 cup winner of slow challenge

You may also like these posts:

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

Implementation in the `network` class