Logistic Regression with R

Logistic Regression

Regression is performed when you want to produce a function that will predict the value of something you don't know (the dependent variable) on the basis of a collection of things you do know (the independent variables).

The problem is that regression is typically done with a linear function and very few real world processes are linear. Hence, a great deal of statistics and machine learning research concerns methods for fitting non-linear functions, but controlling for the explosion in complexity that comes with it.

Logistic Regression is one of the methods that tries to solve this problem. In particular, Logistic Regression produces an output between 0 and 1 which can be interpreted as the probability of your target event happening.

Let's look at the form of Logistic Regression to get a better understanding:

You start with the goal of a function that approximates the probability of the target T for any input vector X :

p(T) = F(X)

In order to assure that F(X) takes the form of a valid probability (i.e. always between 0 and 1) we make us of the logistic function 1/(1+e^-K). If K is a big number the e^-K approaches 0 and hence the output of the logistic function approaches 1. If on the other hand K is a very small number the e^-K becomes very large and the output of the logistic function approaches 0.


So we are fitting the following function:

p(T) = 1 / [ 1 + e^-g(X) ]

We have added the function g(X) to afford us some flexibility in how we feed the input vector X into the logistic function. Here is where we place our usual linear regression function. We say

g(X) = B_0 + B_1 * X_1 + B_2 * X_2 + ...... + B_N * X_N

i.e. a linear function over all the dimensions of the input X.

Now, in order to perform our linear regression, we need to transform the function definition. You can do the transformation yourself if you like. What you will find is that with some re-arrangement you find that the function g(X) is equal to:

g(X) = - ln [ (1-p) / p ]

And by exploiting the properties of the logarithm you can further re-arrange to get the log odds ratio.

g(X) = ln [ p / (1-p) ]

An astute reader might notice a problem. For a target value of 1 (i.e p=1) the fraction is undefined. Luckily we can use the properties of the logarithm again and define our target as

 ln [ p / (1-p) ] =  ln(p) - ln(1-p)

...and this is the target value onto which you perform the linear regression.

In other words you fit the value of the parameters (the Bs) so that

B_0 + B_1*X_1 + B_2*X_2 + ...... + B_N*X_N = ln(p) - ln(1-p)


That is all well and good, how can we do that with R ? you might ask.

Well, I have gone ahead and converted some code from a bunch of different tutorials into a little R workbook that will take you through applied Logistic Regression in R. You can find the Logistic Regression Code Example in my GitHub account right here.

It all boils down to using the Generalised Linear Model function.

This R function will fit your Logistic Regression for you.

If you follow that code example to the end you will get a plot like the one below, which shows you the original data in green, the model fitted to that data in black, and some predictions for unseen parts of the input space in red.

Logistic Regression allows you a great deal of flexibility in your model. The parameterized linear model can be changed how you want, adding or removing independent variables. You can even add higher order combinations of the independent variables.

A common Machine Learning process is to experiment with different forms of this model and examine how the statistical significance of the fit changes.

Just be wary of the pernicious problem of over-fitting.

Maintaining Constant Probability of Data Loss with Increased Cluster Size in HDFS

In a conversation with a colleague some months ago I was asked if I knew how to scale the replication factor of a Hadoop Distributed File System (HDFS) cluster as the number of nodes increased in order to keep the probability of experiencing any data loss below a certain threshold. My initial reaction to the question was that it would not be affected, I was naively thinking the data loss probability was a product of the replication factor only.

Thankfully, it didn't take me long to realize I was wrong. What is confusing, is that for a constant replication factor as the cluster grows the probability of data loss increases, but the quantity of data lost decreases (if the quantity of data remains constant).

To see why consider a situation in which we have N nodes in a cluster with replication factor K. We let the probability of a single node failing in a given time period be X. This time period needs to be sufficiently small so that we know that the server administrator will not have enough time to replace the machine or drive and recover the data. The probability of experiencing data loss in that time period is the probability of getting K or more nodes failing. The exact value of which is calculated with the following sum:

Although in general a good approximation (a consistent overestimate) is simply:

Clearly as the size of N increases this probability must get bigger.

This got me thinking about how to determine the way the replication factor should scale with the cluster size in order to keep the probability of data loss constant (ignoring the quantity). This problem may have been solved elsewhere, but it was an enjoyable mathematical exercise to go through.

In essence we want to know if the number of nodes in the cluster increases by some value n, then what is the minimum number k such that the probability of data loss remains the same or smaller. Using the approximation from above we can express this as:

Now if we substitute in the formulas for N-choose-K and perform some simplifications we can transform this into:

I optimistically thought that it might be possible to simplify this using Stirling's Approximation, but I am now fairly certain that this is not possible. Ideally we would be able to express k in terms of N,n,K,X, but I do not think that it is possible. If you are reading this and can see that I am wrong please show me how.

In order to get a sense of the relationship between n and k I decided to do some quick numerical simulations in R to have a look at how k scales with n.

I tried various combinations of X, N and K. Interestingly for a constant X the scaling was fairly robust when you varied the initial values of N and K. I have plotted the results for three different values of X so you can see the effect of different probability of machine failure. In all three plots the baseline case was a cluster of 10 nodes with a replication factor of 3.

You can grab the R code used to generate these plots from my GitHub repository.

Using AWK for Data Science

Over the years I have become convinced that one of the essential tools needed by anyone whose job consists of working with data is the unix scripting language AWK. It will save you an awful lot of time when it comes to processing raw text data files.

For example, taking a large delimited file of some sort and pre-processing its columns to pull out just the data you want, perform basic calculations or prepare it for entry into a program that requires a specific format.

AWK has saved me countless hours over the years, so now I am writing a super brief primer that should not only convince you it is worth learning but show you some examples.

The first thing you need to know about AWK is that it is data driven, unlike most other languages for which execution is constrained largely by procedural layout of the instructions. AWK instructions are defined by patterns in the data to which actions should be applied. If you are familiar with the regular expression type control structures available in PERL then this should seem like a comfortable idea.

The programs are also data driven in the sense that the entire program is applied to every line of the file (as long as there are patterns that match) and furthermore the program has inbuilt access to the columns of data inside the file through the $0, $1, 2 ... variables: where $0 contains the entire line and $1 upwards has the data from individuals columns. By default the columns are expected to be TAB separated, but you can follow your AWK script with FS=',' to use a comma or any other field separator.

To run a simple AWK script type:

   >awk 'AWK_SCRIPT_HERE' FILE_TO_PROCESS

The basic syntax of the scripts themselves consists of mutiple pattern action pairs defined like this:
PATTERN {ACTION}

One need not include an PATTERN, in which case the action will be applied to every line inside the file to which the program is applied.

So for example, the following program will out the sum of columns 3 and 4

>awk '{print $3+$4}' FILENAME

If we only wanted this to happen when column 1 contained the value 'COSTS' we have a number of options. We could simply use the pattern equivalent of an IF statement as follows:

>awk '$1=="COSTS" {print $3+$4} FILENAME

Alternatively we could use a PATTERN expression as follows

>awk '/COSTS/ {print $3+$4} FILENAME

The problem with the second solution is that it if for some reason the word COSTS can appear in other fields or places in the file then may not get the results you are looking for. There can be a trade off for using the power and flexibility of the regular expression patterns and their ability to lull us into a false sense of security about what they are doing.

There are several special execution paths that can be included in the program. In place of the pattern you can include the reserved words BEGIN or END in order to execute routine before or after the file processing occurs. This is particularly useful for doing something like calculating a MEAN, shown below:

>awk '{sum+=$1; count+=1} END {print sum/count}' FILENAME

By now you should be seeing the appeal of AWK. You can manipulate your data quickly with small scripts that do not require loading an enormous file into a spreadsheet, or writing a more complicated JAVA or PYTHON program.

Finally here are a few of the kinds of tasks that I do with AWK all the time

1)  Convert some file with 10 or more columns into one with a sum of a few and reformating the others:

>awk '{print toupper($1) "," ($3/100) "," ($2+$4-$5)}' FILENAME


2) Calculate the Mean and Standard Deviation on a column. (The following is fora sample, just change the n-1 to n for a complete population.

> awk 'pass==1 {sum+=$1; n+=1} pass==2 {mean=sum/(n-1); ssd+=($1-mean)*($1-mean)} END {print sqrt(ssd/(n-1))}' pass=1 FILENAME pass=2 FILENAME


3) Calculate the Pearson correlation coefficient between a pair of columns. Again for a sample of the data. Change n-1 to n to do the calculation on the entire population data.

> awk 'pass==1 {sx+=$1; sy+=$2; n+=1} pass==2 {mx=sx/(n-1); my=sy/(n-1); cov+=($1-mx)*($2-my); ssdx+=($1-mx)*($1-mx); ssdy+=($2-my)*($2-my);} END {print cov / ( sqrt(ssdx) * sqrt(ssdy) ) }' pass=1 FILENAME pass=2 FILENAME

If you have any great tips for getting more out of AWK let me know, I am always looking for shortcuts.

Machine Learning for Hackers

I recently read "Machine Learning for Hackers" by Drew Conway and John Myles White.

I'd picked it up because I heard it was a good way to get familiar with the data mining capabilities of R. I also expected the case study based approach to be a good way to see how they approach a broad array of machine learning problems. In these respects I was reasonably well rewarded. You will find a bunch of R code scraps that can be reused with a little effort. Unfortunately the explanation of what the code does (and how) is often absent. In this sense the book is true to its name: you will learn some recipes for tackling certain problems, but you may not understand how the code works, let alone the technique being applied.

The one issue I found unforgivable is that in the instances where the authors talk about machine learning theory, or use its terms, they are often wrong. One example is the application of naive Bayes to spam classification. The scoring function they use is the commonly used likelihood times the prior, leaving off the evidence divisor.

As a method of scoring in Bayesian methods this is appropriate because it is proportional to calculating the full posterior probability, and much more efficient to compute. However, the resulting score is not a probability, yet the authors continuously refer to it as one. This may seem minor, but to me it undermined my confidence in their ability to communicate necessary details about the techniques they are applying.

Another example: in the section on distance metrics the authors state that multiplying a matrix by its transpose computes “the correlation between every pair of columns in the original matrix.” This is also wrong. What they want to say is that it produces a matrix of scores that indicate the correlation between the rows. It is an approximation because the score depends on the length of the columns and whether they have been normalised. These values would not be comparable between matrices. What would be comparable between matrices is a correlation coefficient, but this is not what is being computed.

I am not suggesting that a hacker's guide to machine learning should include a thorough theoretical treatment of the subject. I think only that where terms and theory are introduced they should be used correctly. By this criteria this book is a failure. However, for my purposes (grabbing some code snippets for doing analysis with R) it was moderately successful. My largest disappointment was that given the mistakes I noticed regarding the topics about which I have reasonable knowledge, I have no confidence in their explanation of those areas where I am ignorant.

Maximum Likelihood Estimation

Maximum Likelihood Estimation is widely applicable method for estimating the parameters of a probabilistic model.

Developed by R.A.Fisher in the 1920s the principle behind it is that the ideal parameter settings of a model are the ones that make the observed data most likely.

It is applicable in any situation in which the model can be specified such that the probability of the desired variable y can be expressed as a parameterised function over the vector of observed variables (X).

P(y|X) = f(X,φ)

The parameters φ of the function f(X, φ) are what we want to estimate.

The model is designed to be a function in which the parameters are set and we get back a probability value for a given x. However, we need a process to determine these model parameters. The Likelihood function is defined to be equal to this function, but operating as a function over the parameter space of φ.

L(φ | y,X )= P(y|X,φ)

It is important to recognise that Likelihood is not the probability of the parameters, it is just equal to the probability of y given the parameters. As such it is not a probability distribution over φ.

If we have N observations in our data set, and we let D represent all N of these observations of X and y, then we can express the Likelihood function for this entire data set D as :

L(φ | D ) = ∏^N_i=1 P(y_i|X_i,φ )

Maximum Likelihood is then simply defined as Argmax φ over this function. Finding the value of φ that maximises this function can be done a number of ways.

To find an analytical solution to the Likelihood equation we find the partial derivative of the function with respect to each of the paramters. We then solve this series of equations for the parameter values such the the partial derivatives are equal to zero. This gives us a position that is either a max or min. We then find the second partial derivative with respect to each parameter and make sure it is negative at the points found in the first step. This will give us an analytical peak on the Likelihood surface.

The reality of maximising the Likelihood by searching the parameter space depends a great deal on the problem. Numerous tricks occur to simplify the problem. The natural logarithm of the Likelihood function is often taken because they are monotonically related, so the MLE can be obtained by maximising the log of the Likelihood. In addition, taking the log turns the product into a Sum and can improves the chance of finding an analytical solution, and improve the computational tractability of finding a numerical solution.

In the next post I will summarise the use of the Expectation Maximisation algorithm for situations in which the Likelihood function cannot be solved analytically.