Shell-Based ETL for Algorithm Change Tests

Introduction

As promised, I will be adding a little ETL into the mix of my posts. What is ETL you ask? Well, it stands for "Extract Transform Load" and is shorthand for what some might call "data munging".

Throughout the course of data science/analytics work, we deal with a lot of data. And it is not always in the right format or the right place for us to distill insights from it. This is where ETL comes in: it is the process of 'extracting' data from whatever location and form it came, 'transforming it' and then 'loading' it into the form and place you need it to perform your analysis.

There are several software packages available that are devoted solely to the process of ETL. I'll get into some of those in the future. But using pretty minimal tools available at the shell (yes, I'm assuming you are using some version of Linux/Unix/OSX) can be pretty powerful. In this post, I'll show a short example of a couple of shell commands that are key to doing simple one-off data-munging tasks.

The scenario

I had been working on a custom data analytics tool that relies on an iterative optimization method. I was writing custom code to improve the way that the tool found solutions to the optimization problem. Thus, I could use the code to give me information about the 'before' and 'after': that is, 'before I added the new functionality' and 'after I added the new functionality'. I wanted to verify that the change resulted in a significant improvement in the analytics tool by adding the new functionality.

In the end, the functionality was trying to 'keep a value low' (sorry for being vague, but I cannot reveal exactly what I was doing due to confidentiality constraints). Since it was iterative, I knew that the tool's output included the iteration number (which we will call I) as well as the quantity that we were 'trying to keep low', which we'll call Y for brevity. If my custom update was a success, the value that we are trying to keep low would be significantly lower in the end of the process than if my added code were not there. The tool also output a bunch of irrelevant information that I wanted to filter out.

Using 'grep'

I would have to identify text in the output that would act as a signal for my munging method to extract only the pertinent lines and not all of the other mess. This is a job for the shell utility 'grep'. It simply either selects or filters out lines of text that match a regular expression. In my case, I noticed that the characters '>>>' were on each line that output both I (the iteration number value) and Y. Therefore a simple grep expression would suffice to eliminate all of the lines that do not contain '>>>', like so:

./run_analytics_tool | grep '>>>'

Note that I am running the analytics tool with './run_{analyticstool'} and then I am 'piping' the output into 'grep' to select only the lines that contain the string '>>>'. What this gives me is a line that looks something like this:

>>>I,X1,X2,X3,X4,Y
>>>1,18.837,18.279,1,1,1
>>>2,18.557,18.279,1.94,1,1.04
>>>3,18.232,14.239,2.8,2,1.04
>>>4,18.079,14.239,3.18,2,1.28
>>>5,16.312,12.705,3.48,3,2.12
>>>6,14.266,12.705,3.48,3,3.28
>>>7,13.368,10.339,3.44,4,4.28
>>>8,12.761,10.339,3.78,4,4.86
>>>9,12.24,10.339,4.42,4,4.72
>>>10,11.234,10.339,4.74,4,4
>>>11,10.732,8.851,5.3,5,3.96
>>>12,10.972,8.851,6.4,5,4.04
>>>13,10.46,8.851,5.9,5,4.94
>>>14,9.537,8.851,4.84,5,5.8
>>>15,9.378,8.851,5.46,5,5.76
>>>16,9.401,8.851,6.06,6,5.76
>>>17,9.412,8.851,6.86,6,5.66
>>>18,9.386,8.851,7.52,6,5.7
>>>19,9.444,8.851,9.06,6,5.66
>>>20,9.444,8.851,9.06,6,5.66

We don't care about the X1-X4 values but just the I and Y values. So we need some way to select only specific columns of comma-separated values. This is where 'awk' comes in.

Using 'awk'

Awk is a tool that allows us to select delimited columns and rows and operated on them. It is an extremely useful ETL tool.

In our case, we need to simply select certain columns (columns 1 and 6) from a stream of comma-delimited values. The way that we can do this is as follows:

./run_analytics_tool | grep '>>>' | awk -F "," '{OFS=","} {print $1,$6;}'

This gives us the following:

>>>I,Y
>>>1,1
>>>2,1.04
>>>3,1.04
>>>4,1.28
>>>5,2.12
>>>6,3.28
>>>7,4.28
>>>8,4.86
>>>9,4.72
>>>1,4
>>>1,3.96
>>>12,4.04
>>>13,4.94
>>>14,5.8
>>>15,5.76
>>>16,5.76
>>>17,5.66
>>>18,5.7
>>>19,5.66
>>>20,5.66

Here I am using the awk switch "-F" to say that we want to define the incoming delimiter as a comma rather than the default of a space. Similarly I am defining that we want a comma-delimiter as our "output file separator" (thus OFS) with the first section of the awk command in single quotes '{OFS=","}'. Then I am telling awk to print just the first and sixth values for each line in the stream ( '{print $1,$6;}' ). Note that again, I am 'piping' values from grep in a cascading manner to be fed into awk. Now there is a problem: our first column has a pesky string at the beginning; the one we used to our advantage for filtering with grep: '>>>'. We'd like to remove that. For this task, 'sed' is our tool…

Using 'sed'

Sed, the 'stream editor', is a replacement tool. It allows us to operate on a stream of text and transform it by matching with a regular expression and replacing what is matched with something else. We need to perform the simple task of matching '>>>' on each line and replacing it with nothing, which is itself a string: ''. We do this as follows:

./run_analytics_tool | grep '>>>' | awk -F "," '{OFS=","} {print $1,$6;}' | sed 's/>>>//'

Again, I am piping the text from the output of awk to the input of sed. What this sed expression says is to 'substitute' (thus 's') the value between the first two forward slashes (>>>) with the value in the second and third forward slashes (i.e. the empty string). What this gives us is we wanted:

I,Y
1,1
2,1.04
3,1.04
4,1.28
5,2.12
6,3.28
7,4.28
8,4.86
9,4.72
1,4
1,3.96
12,4.04
13,4.94
14,5.8
15,5.76
16,5.76
17,5.66
18,5.7
19,5.66
20,5.66

Now to see if there is a difference

What we ultimately wanted was a way to test whether the change made a difference in our analytics tool. As such, we would like to generate a set of independent 'before' and 'after' samples that we can test statistically. The way we do this is by wrapping all of what we've done until now into a loop and output the results to file. This is done in the following way:

for i in {1..100}
do
 run_analytics_tool | grep '^>>>' | awk -F "," '{OFS=","} {print $1,$6;}' | sed 's/>>>//' > test${i}_before.csv;
done

## make the change to the analytics tool codebase here
for i in {1..100}
do
 run_analytics_tool | grep '^>>>' | awk -F "," '{OFS=","} {print $1,$6;}' | sed 's/>>>//' > test${i}_after.csv;
done

This gives us 100 samples of the analytics tool run. We do this with the analytics tool before adding the new functionality and after, giving us a total of 200 samples with two treatments. Since this is an iterative tool, we cannot assume that the values found in the iterations themselves are independent, so we will take the last value of Y for each of our independent runs (each iteration of the for loop). Then we'll do a statistical test on those using R to see if the change ended up doing what we wanted to, which was to keep the value Y lower at the end of the iteration runs than if my change weren't implemented in the code.

To do this we will import the data into R from the files that we created:

library(plyr)
df.before <- ldply(list.files(pattern="test.*_before.csv"),
       function(filename) {
        dum = read.csv(filename)
        dum$filename = filename
        return(dum)
       })
df.after <-  ldply(list.files(pattern="test.*_after.csv"),
       function(filename) {
        dum = read.csv(filename)
        dum$filename = filename
        return(dum)
       })

## get it into a single data frame
## but preserve an indicator of which 'treatment' it came from
df.all <- rbind(cbind(df.before, treatment="before"),
    cbind(df.after,  treatment="after"))

Since we are only comparing the end iterations, let's just extract those and drop the rest:

df.all.last <- df.all[df.all$generation==max(df.all$generation),]

And we will assume normality for purposes of illustration and do a t-test to determine whether the change resulted in an improvement in the analytics tool:

t.test(ave.complexity~treatment, df.all.last)

 Welch Two Sample t-test

data:  ave.complexity by treatment
t = 22.301, df = 118.55, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 6.718873 8.028327
sample estimates:
mean in group before  mean in group after 
             10.6382               3.2646

Conclusion

In this post, we demonstrated basic usage of shell-based 'ETL tools': grep, awk and sed. We used the scenario in which we made a change to an algorithm within an analytics tool. This tool had a good deal of useful information in its output as well as a lot of junk that we didn't care about. We filtered out the junk and obtained information that was useful using munging methods. From this information we compared the analytics tool before the algorithm change and after and demonstrated that the change resulted in a significant improvement to the tool.

Data Tables in R

Data Tables in R

Introduction

My posts have been pretty machine learning and data mining heavy until now. And although machine learning and data mining are arguably the most important aspect of what we have come to call 'data science', there are a couple other components that we cannot ignore: how we store data and how we 'munge' or manipulate data. I promise I will get to these aspects of data science in future posts. In a segue into these topics, today's post will be on a package that I've recently discovered that is a useful way to store data in R (well, at least in memory) and get useful aggregate information from it.

Usually we work with the 'data frame' class in R to store data that we are working with in memory. Data frames are one of the reasons that I choose to do my analytics in R versus another environment such as Matlab (although the folks at Mathworks - the makers of Matlab - have added a data frame-like type to their matrices in recent years but it doesn't match the power of the R version). Data frames are useful because - amongst other reasons - we can work with data in a natural way. Columns are named and typed into numeric, character or factor values. And we can store mixed data in data frames - not all columns are required to be the same type (as is the case for matrices).

However, I've recently discovered the 'data.table' library. I was looking for a quick way to operate on grouped data in a data frame and data tables seemed to fit the bill. I use plyr all the time for this type of operation; and plyr is very good at it (if you don't use it, I would recommend starting). But I was doing something very simple (finding min and max by a factor) and the data.table library came up in a search. And it looks like a package that I will probably start incorporating into my set of go-to libraries in R.

Oh, and I know there is a newer version of plyr (dplyr), but I haven't had time to make the switch. Someday…

Example Data

I'll generate some simple example data for demonstration purposes. The data will start its life as a data frame. It is just two columns: a numeric column of uniformly-distributed random values and a column of groups, randomly chosen:

df <- data.frame(x=runif(20),
     y=sample(c("A","B"), size=20, replace=TRUE))
print(df[order(df$y),])

0.970504403812811	A
0.64111499930732	A
0.915523204486817	A
0.510243171826005	A
0.0864250543527305	A
0.104028517380357	A
0.230543906101957	A
0.274059027433395	A
0.825932029169053	A
0.151608517859131	B
0.0283026716206223	B
0.606306979199871	B
0.0992447617463768	B
0.313594420906156	B
0.610857902327552	B
0.588633166160434	B
0.141868675360456	B
0.140498693101108	B
0.219244148349389	B
0.664096125168726	B

The plyr version

So what I would usually do in plyr to find, say, the min value for each group is something like this:

library(plyr)
df.agg <- ddply(df, c("y"), summarize, min.val.per.group=min(x))

A	0.0864250543527305
B	0.0283026716206223

As you can see, this gives us a new data frame with the minimum values for each group, A and B. A similar thing could be done for the max values.

Don't get me wrong, I love plyr and will continue to use it. But data.tables are a different approach to working with grouped frames of data that is intriguing.

The data.table version

In contrast to data frames, data tables maintain a "key" that can be made up of the columns of the table.

First we'll convert our original data frame to a data table. Then we set the key to be the "y" column (remember, the key can actually be a combination of columns if desired). Lastly, we use the tables() function to print out the tables that we have access to in memory (just our dt table right now) and some info on that table, such as its key.

library(data.table)
dt <- as.data.table(df)
setkey(dt, y)
tables()

dt

20

2

1

x,y

y

Having a key gives us access to a simple and concise means by which we can group information on the table as well as filter it. If we only want the "A" group values on our table, we can filter it in this way:

dt["A",]

0.970504403812811	A
0.64111499930732	A
0.915523204486817	A
0.510243171826005	A
0.0864250543527305	A
0.104028517380357	A
0.230543906101957	A
0.274059027433395	A
0.825932029169053	A

Also note that the comma is optional:

dt["A"]

0.970504403812811	A
0.64111499930732	A
0.915523204486817	A
0.510243171826005	A
0.0864250543527305	A
0.104028517380357	A
0.230543906101957	A
0.274059027433395	A
0.825932029169053	A

And we can perform actions on the groups using the second index into our table. So if we want to take the min of the "A" group, we can do it like this:

dt["A", min(x)]

0.0864250543527305

But we can also get the minimum value for each key (or in our case since we defined the key as just the y column, which delineates our groups) succinctly by using the additional index into the data table with the by keyword:

dt[, min(x), by=y]

A	0.0864250543527305
B	0.0283026716206223

And notice that we don't have to address our columns explicitly as columns of the table (e.g. dt$x or dt$y), making things even more concise.

Conclusion

We only went through a brief and simple introduction of the data.table library in R. But hopefully that is enough to make you interested in pursuing this fairly simple but useful library further. I, for one, will be continuing to use it in my daily routine of analytics work thanks to its simplicity and ability to operate on groups of rows in such a succinct way.

K-Fold Cross Validation with Decision Trees in R

1 K-Fold Cross Validation with Decisions Trees in R   decision_trees machine_learning

1.1 Overview

We are going to go through an example of a k-fold cross validation experiment using a decision tree classifier in R.

K-fold cross validation is a method for ensuring a robust error estimate on a trained classification model.

When we train a predictive model, we want that model to not only be accurate on the data that we used to train the model but also generalize to other samples that the model has not yet been presented with. A common technique for ensuring this generalizability is to split data into training data and test data sets. The model is trained on the training data split and then tested on the test dataset to ensure that the model did not only learn to be accurate on the training dataset (overfit).

Similarly, in k-fold cross validation we split the data into k equally-partitioned subsamples. Then for each of the k partitions, we hold out the \(i^{th}\) partition and train our model on the other \(k-1\) partitions and test on the \(i^{th}\) partition. We then average the error over the testing results of all of our k rounds of training/testing.

1.2 Naive Training/Testing

To begin, we will show a naive implementation of a train/test process. In this example, we don't split between the training data and the testing data.

Here, we are training the model on the full dataset.

library(rpart)
data(iris)
rpart.model <- rpart(Species~., data=iris, method="class")
print(rpart.model)

n= 150 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)  
  2) Petal.Length< 2.45 50   0 setosa (1.00000000 0.00000000 0.00000000) *
  3) Petal.Length>=2.45 100  50 versicolor (0.00000000 0.50000000 0.50000000)  
    6) Petal.Width< 1.75 54   5 versicolor (0.00000000 0.90740741 0.09259259) *
    7) Petal.Width>=1.75 46   1 virginica (0.00000000 0.02173913 0.97826087) *

And now we test on the same dataset. From this, we obtain a confusion matrix.

rcart.prediction <- predict(rpart.model, newdata=iris, type="class")
confusion.matrix <- table(iris$Species, rcart.prediction)
print(confusion.matrix)

          rcart.prediction
           setosa versicolor virginica
setosa         50          0         0
versicolor      0         49         1
virginica       0          5        45

The resulting error is as follows:

accuracy.percent <- 100*sum(diag(confusion.matrix))/sum(confusion.matrix)
print(paste("accuracy:",accuracy.percent,"%"))

[1] "accuracy: 96 %"

Pretty good. But our model could very well be overfit. If we were to obtain new measurements for each of these species our accuracy might not be very good because our model fits the data that we trained on well but does not generalize for new data.

1.3 Using k-fold cross-validation to train and test the model

So let's use k-fold cross-validation to obtain a more generalizable model.

library(plyr)
library(rpart)
set.seed(123)
form <- "Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width"
folds <- split(iris, cut(sample(1:nrow(iris)),10))
errs <- rep(NA, length(folds))

for (i in 1:length(folds)) {
 test <- ldply(folds[i], data.frame)
 train <- ldply(folds[-i], data.frame)
 tmp.model <- rpart(form , train, method = "class")
 tmp.predict <- predict(tmp.model, newdata = test, type = "class")
 conf.mat <- table(test$Species, tmp.predict)
 errs[i] <- 1-sum(diag(conf.mat))/sum(conf.mat)
}
print(sprintf("average error using k-fold cross-validation: %.3f percent", 100*mean(errs)))

[1] "average error using k-fold cross-validation: 7.333 percent"

So there we have it. K-fold cross-validation in action. We can see that the error increased when using k-fold cross-validation over simply training and then testing on the same data, which indicates that there may have been bias introduced by overfitting the model in the latter case.

1.4 k-fold cross-validation with C5.0

Let's do the same thing but with a different decision tree algorithm, C5.0. This is an update of J. Ross Quinlan's popular C4.5 algorithm.

library(C50)
library(plyr)
errs.c50 <- rep(NA, length(folds))
form <- "Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width"
folds <- split(iris, cut(sample(1:nrow(iris)),10))
for (i in 1:length(folds)) {
 test <- ldply(folds[i], data.frame)
 train <- ldply(folds[-i], data.frame)
 tmp.model <- C5.0(as.formula(form), train)
 tmp.predict <- predict(tmp.model, newdata=test)
 conf.mat <- table(test$Species, tmp.predict)
 errs.c50[i] <- 1 - sum(diag(conf.mat))/sum(conf.mat)
}

print(sprintf("average error using k-fold cross validation and C5.0 decision tree algorithm: %.3f percent", 100*mean(errs.c50)))

[1] "average error using k-fold cross validation and C5.0 decision tree algorithm: 6.000 percent"