CYBAEA Data and Analysis

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

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

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. This is worse, much worse (you only need about 6 times as many workers to reduce the productivity to ¼).

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.

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.
	seconds	index	seconds	index	seconds	index	seconds	index	seconds	index	seconds	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

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?

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.

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

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.

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

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!

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.

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?

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.

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?

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.

R tips: Keep your packages up-to-date

2009-03-25T20:59:00Z

In this entry in a small series of tips for the use of the R statistical analysis and computing tool, we look at how to keep your addon packages up-to-date.

One of the great strengths of R is the many packages available. All the new approaches, as well as some of the best implementations of your old favorites are there. But it can also be a little daunting, and so the CRAN task views are often the best way to get started and download a reasonable “bundle” of packages for your analysis.

First we need a place to store the packages. On Linux (and other Unix-like systems) I use the file ~/.Renviron to set the R_LIBS variable to where I want the files:

## R environment
R_LIBS="~/R"

On Windows, I set the same variable for the user account. Don’t forget to create the directory.

Now your can start R and install the CRAN task view package:

> install.packages("ctv")

Then I have a few things in my ~/.Rprofile startup file. The previous command probably prompted you for a download mirror which is annoying, so let’s exit R and edit the startup file to contain:

## Default CRAN mirror
local({r <- getOption("repos"); r["CRAN"] <- "http://cran.uk.r-project.org"; options(repos=r)})
## Libraries
require("utils", quietly=TRUE)
require("ctv", quietly=TRUE)

Then I define three functions. The first is to install the views I need. I like to try new things, so my list is long. Edit it to suit your needs:

install.myviews <- function() {
  require("ctv", quietly=TRUE)
  my.views = c("Bayesian", "Cluster", "Graphics", "gR", "HighPerformanceComputing", "MachineLearning", "Multivariate", "NaturalLanguageProcessing", "Robust", "SocialSciences", "Spatial", "Survival", "TimeSeries")
  install.views(views=my.views, lib=Sys.getenv("R_LIBS"), dependencies=c("Depends","Suggests"))
}

Try it out! Save the file, start R, and type install.myviews() at the prompt. If your list is as long as mine, then this may take some time and you may get some warnings and errors. We might add a tip on these later, but the main reason for the errors is probably that you are missing the development files for external libraries (or that R just can’t find it).

Now that we have finally got them, we need to make sure they are up-to-date. I add two functions to ~/.Rprofile:

update.local <- function() {
  update.packages(lib.loc=Sys.getenv("R_LIBS"), ask=FALSE)
}

update.myviews <- function() {
  require("ctv", quietly=TRUE)
  my.views = c("Bayesian", "Cluster", "Graphics", "gR", "HighPerformanceComputing", "MachineLearning", "Multivariate", "NaturalLanguageProcessing", "Robust", "SocialSciences", "Spatial", "Survival", "TimeSeries")
  update.views(views=my.views, lib.loc=Sys.getenv("R_LIBS"))
}

The first allows me to easily update all my locally installed libraries (not just these installed from views). The second updates my views which is useful when the view definitions change (rarely, but it happens as the recommended packages evolve).

Now I can of course update from the R command prompt using update.local() or update.myviews(). But that is not the main benefit. I can now update directly from the shell command line using commands like:

echo "update.local()" > /tmp/r.cmd
R CMD BATCH /tmp/r.cmd /tmp/r.out

The beauty of this is that I can add it to my crontab(5) and have it run automatically every night or every week as I feel I need it. This way I always have the latest versions installed.