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.
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (May 4, 2015). Working with “large” datasets, with dplyr and data.table. Freakonometrics. Retrieved September 16, 2024 from https://doi.org/10.58079/ouzl
Nice article on the usage of dplyr for data manipulation.
https://honingds.com/blog/data-manipulation-with-dplyr/
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()
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
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!
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