CYBAEA Data and Analysis

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 ) {
	  print("Solution found")
	  print(c[x,])
	  break
	}
	r <- r + 1L
  }
})

appetizer.solution(15.05)
# [1] "Solution found"
# [1] 3.55 5.80 5.80

Brute force works, it just doesn’t scale well.

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.

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.

First install the openmpi libraries using:

yum install openmpi openmpi-devel openmpi-libs

Then 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

Now it 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

Painful.

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.