fond-carte

Working with “large” datasets, with dplyr and data.table

A few months ago, I was doing some training on data science for actuaries, and I started to get interesting puzzeling questions. For instance, Fleur was working on telematic data, and she’s been challenging my (rudimentary) knowledge of R. As claimed by Donald Knuth, “we should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil“. So usually, in my courses, and my training, codes are very basic, and easy to understand. But usually poorly efficient. Since I was challenged, to work on very large datasets, we’ve been working on R functions to manipulate those possibly (very) large dataset, and to run some simple functions as fast as possible (with simple filter and aggregation functions).

In order to illustrate, let us generate our “large” telematic dataset. Assume that we have 10,000 drivers, each of them drives about 200 times, and each time, we have, say, 80 locations. That mean around 160 million observations. It is “large”, but not huge.

> rm(list=ls())
> N_id=10000
> N_tr=200
> T_tr=80

In order to have a code as general as possible, assume that we have some kind of randomness,

> set.seed(1)
> N=rpois(N_id,N_tr)
> N_traj=rpois(sum(N),T_tr)

By “observation”, we consider a driver Id., a Trajectory Id., and a location (latitude and longitude) at some specific dates (e.g. every 15 sec.). Again, just because we want some dataset to illustrate, swe will draw drivers’s home randomly (here uniformly on some square)

> origin_lat=runif(N_id,-5,5)
> origin_lon=runif(N_id,-5,5)

And, then, from those locations, we generate a 2-dimensional random walk,

> lat=lon=Traj_Id=rep(NA,sum(N_traj))
> Pers_Id=rep(NA,length(N_traj))
> s=1
> for(i in 1:N_id){Pers_Id[s:(s+N[i])]=i;s=s+N[i]}
> s=1
> for(i in 1:length(N_traj)){lat[s:(s+N_traj[i])]=origin_lat[Pers_Id[i]]+
+  cumsum(c(0,rnorm(N_traj[i]-1,0,sd=.2)));
+  lon[s:(s+N_traj[i])]=origin_lon[Pers_Id[i]]+
+  cumsum(c(0,rnorm(N_traj[i]-1,0,sd=.2)));
+  s=s+N_traj[i]}

We have something which looks like

Then, we create variables for the Driver Id., and the Trajectory Id.,

> Pers_Id=rep(NA,sum(N_traj))
> N0=cumsum(c(0,N))
> s=1
> for(i in 1:N_id){Pers_Id[s:(s+sum(N_traj[(N0[i]+1):N0[i+1]]))]=i;s=s+sum(N_traj[(N0[i]+1):N0[i+1]])} 
> s=1
> for(i in 1:sum(N)){Traj_Id[s:(s+N_traj[i])]=i;s=s+N_traj[i]}

Now, we should store that dataset. Consider for instance a standard data frame,

> df=data.frame(Pers_Id,Traj_Id,lat,lon)

But in order to run faster codes (or supposed to be faster), we can consider a local data frame,

> library(dplyr)
> ldf=data_frame(Pers_Id,Traj_Id,lat,lon)

or a data table,

> library(data.table)
> dt=data.table(Pers_Id,Traj_Id,lat,lon)

Observe that those files are almost the same size,

> object.size(df)
3839908904 bytes

Actually, I have already created those files. The three files are stored on Dropbox, e.g. dt_json_2.RData. From that file, two alternative techniques can be considered

  • store the file on a local drive, and then load it

It will not be that long (even on a basic laptop)

> rm(list=ls())
> library(data.table)
> system.time( load("dt_json_2.RData") )
       user      system     elapsed
      21.53        1.33       27.71
  • load the file directly from Dropbox

That’s a bit tricky, but it is possible

> dropbox_ldf <- "https://dl.dropboxusercontent.com/s/l7rolinwojn37e2/ldf_json_2.RData"
> dropbox_df  <- "https://dl.dropboxusercontent.com/s/wzr19v9pyl7ah0j/df_json_2.RData"
> dropbox_dt  <- "https://dl.dropboxusercontent.com/s/nlwi1mmsxr5f23k/dt_json_2.RData"
 
> source_https <- function(loc,...){
+  curl=getCurlHandle()
+  curlSetOpt(cookiejar="cookies.txt",
+ useragent="Mozilla/5.0", followlocation=TRUE)
+  tmp <- getURLContent(loc, .opts=list(ssl.verifypeer=FALSE),curl=curl)
+  fch <- source(tmp)
+ fch
+ }

Then, you have to be (really) patient (unless you have a very good internet connexion).

Now that we have our dataset, consider the following simple problem. Given a final location, say, close to the (0,0), i.e. final distance to (0,0) smaller than 1, what was the initial location ? One the graph below, given that a trajectory ended in the black circle, what could have been the location of the starting point ?

Here, we can use the three different dataset, that rely on three different packages.

  • Using the (standard) data frame format

Consider here our data frame file,

> system.time(load("df_json_2.RData"))
       user      system     elapsed
      21.76        1.39       29.59

In order to get location of the initial starting point given information on the ending point, let us create two variables

> system.time(n <- nrow(df))
       user      system     elapsed
          0           0           0 
> system.time(df$first<-c(1,df$Traj_Id[2:n]!=
df$Traj_Id[1:(n-1)]))
       user      system     elapsed
       3.12        1.23        4.37 
> system.time(df$last<-c(df$Traj_Id[2:n]!=
df$Traj_Id[1:(n-1)],1))
       user      system     elapsed
       2.90        6.23       31.58

I don’t know why the second one was that long. Maybe because of memory issues. Because our dataset is now quite large

> object.size(df)
6399847720 bytes

On a Windows machine, we might start have problems, because we also want to create a variable to test wether our final location was in the black disk, or not,

> lat_0=0
> lon_0=0
> system.time(df$test <-(lat_0-df$lat)^2+(lon_0-df$lon)^2) <=1)&(df$last==1)
Error: cannot allocate vector of size 1.2 Gb/pre>

Maybe we can skip the creation of the two previous variables

> df$first <- NULL
> df$last <- NULL
> object.size(df)
3839908904 bytes
> system.time(df$test<-(lat_0-df$lat)^2+(lon_0-df$lon)^2) <=1)&(c(df$Traj_Id[2:n]!=
df$Traj_Id[1:(n-1)],1)==1))
      user    system    elapsed 
      9.39     11.10      41.17

Let us create a small vector, with the list of all trajectories

> system.time(list_Traj <- unique(df$Traj_Id[df$test==TRUE]))
       user      system     elapsed
       0.72        0.25        0.98

and then, we can extract initial location of all those trajectories

> system.time(base <- df[(df$Traj_Id %in% list_Traj)& (c(1,df$Traj_Id[2:n]!=df$Traj_Id[1:(n-1)])==1),c("lat","lon")])
       user      system     elapsed
       11.7         2.7        14.4 
> head(base)
             lat      lon
556   -0.9597891 2.469243
2866  -0.9597891 2.469243
4321  -0.9597891 2.469243
5677  -0.9597891 2.469243
9403  -0.9597891 2.469243
10432 -0.9597891 2.469243
> nrow(base0
[1]  63453

Even with 160 million lines, that was not too painful. Based on that table, it is possible to get the distribution of the initial location,

> X <- base[,c("lon","lat")]
> library(KernSmooth)
> kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])),gridsize = c(251L, 251L))
> image(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat,col=
rev(heat.colors(100)))
> contour(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE)

Here it took about one minute to extract the information out of those 160 million lines. What about other formats ?

  • Using the data table format

Let us load our data table

> rm(list=ls())
> library(data.table)
> system.time( load("dt_json_2.RData") )
       user      system     elapsed
      21.53        1.33       27.71

First, we should specify the key,

> system.time( setkey(dt,Traj_Id) )
       user      system     elapsed
       0.38        0.09        0.47

A first idea to get the starting point of all trajectories is to use

> system.time(depart <-dt[!duplicated(Traj_Id)]) 
       user      system     elapsed
       2.48        0.65        3.14

but it is possible to use a more elegant code, to get the starting and the ending point

> system.time( depart <- dt[J(unique(Traj_Id)), mult = "first"])
       user      system     elapsed
       2.64        0.75        3.39 
> system.time( arrivee <- dt[J(unique(Traj_Id)), mult = "last"] )
       user      system     elapsed
       2.81        0.51        3.33

Then we can create a distance to (0,0) variable

> lat_0=0
> lon_0=0
> system.time( arrivee[,dist:=(lat-lat_0)^2+(lon-lon_0)^2] )
       user      system     elapsed
       0.03        0.08        1.60

and finally, create a dataset with all trajectories that ended in the black disk

> system.time( fin <- subset(arrivee,dist <= 1) )
       user      system     elapsed
       0.04        0.00        0.78

We can remove variables that we do not need

> system.time( fin[,Pers_Id:=NULL] )
       user      system     elapsed
       0.00        0.00        0.07 
> system.time( fin[,lat:=NULL] )
       user      system     elapsed
        0.0         0.0         0.2 
> system.time( fin[,lon:=NULL] )
       user      system     elapsed
          0           0           0

set the key (just in case)

> system.time( setkey(fin, Traj_Id) )
       user      system     elapsed
       0.00        0.00        0.42

and then, use a merge, with the starting point dataset

> system.time( base <<- merge(fin,depart,all.x=TRUE) )
       user      system     elapsed
       0.02        0.00        0.17

and we’re done

> system.time( base <<- merge(fin,depart,all.x=TRUE) )
       user      system     elapsed
       0.02        0.00        0.17 
> 
> head(base)
   Traj_Id       dist Pers_Id        lat      lon
1:       8 0.41251163       1 -0.9597891 2.469243
2:      36 0.34545373       1 -0.9597891 2.469243
3:      54 0.24766671       1 -0.9597891 2.469243
4:      71 0.00210023       1 -0.9597891 2.469243
5:     117 0.00755432       1 -0.9597891 2.469243
6:     130 0.82806342       1 -0.9597891 2.469243

Hopefully, we have the same dataset as previously. Here, it was quite fast, less than 10 seconds.

  • Using a local data frame

A finaly alternative is to use our local data frame

> rm(list=ls())
> library(dplyr)
> library(data.table)
> system.time( load("ldf_json_2.RData") )
       user      system     elapsed
      27.53        1.77       69.84

We can adapt our previous code here,

> system.time( ldepart <<- ldf %>% group_by(Traj_Id) %>%
+ summarise(first_lat=head(lat,1),
first_lon=head(lon,1)) )
       user      system     elapsed
      60.82        1.46       80.54 

> system.time( larrive <<- ldf %>% group_by(Traj_Id) %>%
+ summarise(last_lat=tail(lat,1),last_lon=tail(lon,1)) )
       user      system     elapsed
      60.81        0.31       62.15 

> lat_0=0
> lon_0=0
> system.time( system.time( larrive <<- mutate(larrive,dist=(last_lat-lat_0)^2+(last_lon-lon_0)^2) ))
       user      system     elapsed
       0.08        0.00        0.09 
> system.time( lfin <<- filter(larrive,dist<=1) )
       user      system     elapsed
       0.05        0.00        0.04

Here also, we end with a merging function

> system.time( lbase <- left_join(lfin,ldepart) )
Joining by: "Traj_Id"
       user      system     elapsed
       0.53        0.05        0.66

One more time, the output contains the same information

> head(lbase)
Source: local data frame [63,453 x 6]
 
  Traj_Id    last_lat     last_lon        dist  first_lat first_lon
1       8  0.41374639  0.491248980 0.412511638 -0.9597891  2.469243
2      36  0.58774352  0.003360806 0.345453735 -0.9597891  2.469243
3      54  0.34479069 -0.358867800 0.247666719 -0.9597891  2.469243
4      71 -0.04341135  0.014686416 0.002100236 -0.9597891  2.469243
5     117 -0.05103682 -0.070353141 0.007554322 -0.9597891  2.469243
6     130 -0.56196768 -0.715720445 0.828063425 -0.9597891  2.469243

And here, it took about 2 minutes… The longest part was to extract those first and last observations. So far, it looks like data.table is just perfect to deal with those “large” datasets.


9 thoughts on “Working with “large” datasets, with dplyr and data.table”

  1. I have re-written and annotated the author’s code to make it more clear what’s being done and also to reduce the coding by using list operations instead of loops. Doing so eliminates the need to initialize variables like Pers_Id and use auxiliary counters like the variable, s. It is also an example of the usage of the apply family of R functions.

    The lat/long are stored in variables lat_n / lon_n which are “ragged” lists. Each element of the list corresponds to a driver is a complete trajectory. This simplifies plotting a subset of trajectories to be used for illustrations.

    The following is an R-script taken from my R session (I use RStudio).

    ## First Figure
    # This code is a modification of the code which produces the first
    # figure in the reference “Working with ‘large’ datasets, with dplyr and
    # data.table” located in:
    # http://www.r-bloggers.com/working-with-large-datasets-with-dplyr-and-data-table/
    # Got rid of the for loops.
    N_id=100 # Number of Drivers
    N_tr=200 # Number of trips per driver
    T_tr=80 # Number of locations visited
    set.seed(1)
    # Vector of number of trips for each driver. The number of trips
    # driven by the ith driver will be N[i]. These are sampled from
    # a poisson distribution with lambda = N_tr.
    N=rpois(N_id,N_tr)
    # sum(N) is the total number of trips for all drivers
    # Poisson distributed number of locations visited for each trip
    # lambda = T_tr
    N_traj=rpois(sum(N),T_tr)
    box=2 # size of the box where drivers originate
    origin_lat=runif(N_id,-box,box) # starting latitude for each driver
    origin_lon=runif(N_id,-box,box) # starting longitude for each driver
    # Allows one to recover the identity of the driver who made a given
    # trip.
    Pers_Id=unlist(sapply(1:N_id,function(x) rep(x,N[x])))
    # Compute the trajectories
    lat_n=lapply(1:sum(N),function(x) cumsum(c(origin_lat[Pers_Id[x]],rnorm(N_traj[x],0,.2))))
    lon_n=lapply(1:sum(N),function(x) cumsum(c(origin_lon[Pers_Id[x]],rnorm(N_traj[x],0,.2))))
    ## Plot some of the results
    # Print a random subset of trajectories
    #
    # sample size
    sampsiz=60
    # sample subset of trajectories
    cnts=sample(1:length(N_traj),sampsiz)
    # The next statement is used to increase the resolution
    # of the graphical output beyond 72 dpi which RStudio uses
    # by default. You will not see a picture in the plot viewer.
    # There will be a plot in your working directory with the name
    # “fname.tiff”. This statement can be skipped.
    tiff(“fname.tiff”, width = 11, height = 8.5, units = ‘in’, res = 216)
    # Set up plot frame
    # Frame box size
    lim=8
    plot(-lim:lim,-lim:lim,xlab=”lat”,ylab=”lon”,type=”n”,main=”Random Walks”)
    # Plot trajectories
    plts=lapply(cnts,function(x) lines(unlist(lat_n[x]),unlist(lon_n[x]),col=sample(rainbow(1000))))
    # Plot points of origin for each trajectory
    points(origin_lat[Pers_Id[cnts]],origin_lon[Pers_Id[cnts]],pch=16)
    # Turn off device mode. If you skipped the tiff() statement
    # then skip this next statement as well.
    dev.off()

  2. You can substantially reduce the computation time when using dplyr if instead of head() and tail() you use first() and last().

    So instead of:
    system.time(ldepart %
    group_by(Traj_Id) %>%
    summarise(first_lat=head(lat,1), first_lon=head(lon,1)))
    user system elapsed
    63.768 1.229 65.115

    use:
    system.time(ldepart %
    group_by(Traj_Id) %>%
    summarise(first_lat=first(lat), first_lon=first(lon)))
    user system elapsed
    11.045 0.940 11.987

    Nelson

  3. Hi,
    I am trying to reproduce your code and I get an error.
    > warnings()
    Warning messages:
    1: In lat[s:(s + N_traj[i])] = origin_lat[Pers_Id[i]] + cumsum(c(0, … :
    número de items para para sustituir no es un múltiplo de la longitud del reemplazo

    > dim(lat)
    NULL
    > dim(lon)
    NULL

    R version 3.1.3 (2015-03-09) — “Smooth Sidewalk”
    Copyright (C) 2015 The R Foundation for Statistical Computing
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    How do you plot the data (I think it’s not in the post)

    Thanks for share!

    1. Hi,
      You can plot the data by using lines() :
      My advice is to limit it to a small part of you dataset unless you’re very patient 🙂

      This one works with classic data.frame:
      for (i in 1:length(N_traj)) {
      lines(df$lon[df$Traj_Id==i], df$lat[df$Traj_Id==i], type=”l”, col=colors()[i])
      }

      Fleur

Leave a Reply

Your email address will not be published. Required fields are marked *