Wednesday, 5 October 2011

Modelling with R: part 3

The previous posts, part 1 and part 2, detailed the procedure to successfully import the data and transform the data so that we can extract some useful information from them. Now it's time to get our hands dirty with some predictive modelling.

The dependent variable here is a binary variable taking values "0" and "1", indicating whether the customer defaulted or not. There are specific, but many, modelling techniques that help us with binary dependent variables. Let's start with the one of the simplest, the Logistic Regression model. And this time, let's just dive straight into it and try to make sense as we go along.

The above code is pretty much self-explanatory. We are running a logistic model on our transformed data. We have specified the "logit" link which means that the estimates (or coefficient or weights) will be for the log of odds ratio and not the direct probabilities. Additionally, note that in case we had to regress the dependent variable, "default" on all the independent variables, we could have simply written# m1.logit <- glm(default ~ ., family = binomial(link = "logit"), data = dev)
Here the "." indicates all the variables in the data set except for the dependent variable. We are not using the command here is because we have the original "amount" and "age" variables as well as the transformed "amt.fac" and "age.fac" variables as well and we should use only one form, either the continuous form or the factor form but not both.

We can see that we have five variables which are significant at 1 percent level of significance and one variable significant at 0.1 percent level of significance. Great! Or not so great! How do we judge how good/bad our model is? Have no fear, R is here.

One of the most popular methods to check the performance of a model, which has a binary independent variable, is by plotting the Receiver Operating Characteristic (ROC) curve which is a plot of the sensitivity, or the true positive rate vs. the false positive rate. In simpler terms, it plots the no. of 1s in the model correctly identified as 1 vs. the no. of 0s in the model which are incorrectly identified as 1.

Also, it makes more sense to plot the ROC curve on the validation sample. For this, we would need to score that validation data set. By "scoring" we mean that with the values of the independent variables given in the validation sample, we will run the model equation on each row of the data set and get the probability of the customer defaulting. Simply, use all the Xs and the Y-hats or the predicted Ys.val$m1.yhat <- predict(m1.logit, val, type = "response")# The predict command runs the regression model on the "val" dataset and stores the estimated y-values, i.e, the yhat.

We can then plot the ROC curve by plotting the performance fuction which checks how our model is performing.plot(performance(m1.scores, "tpr", "fpr"), col = "red")abline(0,1, lty = 8, col = "grey")# "tpr" and "fpr" are arguments of the "performance" function indicating that the plot is between the true positive rate and the false positive rate.# "col" is an argument to the plot function and indicates the colour of the line# "abline" plots the diagonal, "lty" is the line type which is used to create the dashed line

The ROC curve's performance is judged by the area under the curve formed by the ROC curve and the 45' line. The idea is that if we randomly try to classify the customers who will default, we will get the 45' line. And if we use a predictive approach, then the approach should try to classify all the 1s and 1, i.e, correctly identify all the people who defaulted. Such a curve would then be a right triangle with vertices at (0,0), (0,1), (1,1). But since predictive modelling is usually accompanied with some error, we get a curve like we see in the graph. And the rule of thumb is, more the area under the curve, better the model.

In addition to the ROC curve, the KS statistic is also a popular metric used to judge model performance, the Kolmogorov-Smirnov statistic. It is the maximum difference between the cumulative true positive rate and the cumulative false positive rate

I am really sorry for the late reply. I had been travelling for the past 3 days and did not have access to an internet connection.

Anyway, coming back to your problem, we don't need to convert the "response" variable so that it lies between 0 and 1. R will do that for us. But we do need to ensure that the "response" variable is a factor variable and not numeric.

I think you are getting the problem because by default R read the "response" variable to be numeric instead of factor and then it is tried to run a logistic regression on a numeric variable instead of a factor variable. Just change it to factor and you should be good to go. You can use

Congratulations on your blog. It is very useful to us 'data miners'. I had two comments about this post:

1. A very basic doubt. The place where you listed down all the variables in the model equation. Is there a way in R to get a list of all the variables present in the data set. I mean do we have to manually enter the variable name separated by a "+" sign. In case where there are too many variables, this might be a tedious task. For instance, I have about 200 variables in my model, out of which I want to use say 170; how do we go about it?

2. This is more of a possible solution/comment motivated by one of the projects I m working on!!So, towards the end, you talk about "performance" of the model. I was just wondering if we can draw upon the concepts from very basic Econometrics. As we remember, a simple transformation of variables (squared/root/logistic/cube etc.) can sometimes improve the performance of the model. When we are faced with huge data, it can be a herculean task to use intuition for each variable and decide the form in which the variable enters the model. What we did in one of our projects is designed an automated process (read: macro) such that we could transform each variable and check the best form in which it can enter the model equation. What I want to ask is whether there is a way to automate this kind of a process in R. I m sure there is but forgive my ignorance if this is a very silly question! I hope I have made myself clear in my doubt.

Of course this is a very basic and elementary kind of tweaking that I can think of (and also limited in applicability since it can't be used on categorical variables!). After reading your next post on decision trees, I m convinced that you have many more things up your sleeve to help us resolve the "performace" issues.

I want to commend for the succint and simple rendering of complicated modelling for dummies like myself.

1. We don't necessarily have to list down all the variables manually. Like mentioned in the post as well, in case you wanted all the two hundred variables to be included you can simply writelm(y ~ .), where "." stands for all the variables except the dependent one.

Other than that if you wanted to use, say, 170 of these variables then there a couple of convenient options. I'll list the ones I know.

A. All 170 variables are arranged in one order, say, from 20 to 190: lm(y ~ data[, c(20:190)]

B. Some are clubbed together: lm(y ~ data[, c(1:50, 60:100, 120:200)]

C. No proper order:This is a little tedious but is still better than writing the variable names. You specify each column that it should take lm(y ~ data[, c(1,2,4,5,11,17,25)]

2. I'll have to look in to your second question. Will update the post as soon as I find some valuable information.

Hi MK,thanks a ton for taking so much efforts and putting this .. love it !this is really really helpful for a beginner like me..

I ran the whole piece of codes you have written and then this struck when I ran the logistic regression part..Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels

I got confused and started searching for an answer, tried googling to understand but it was in vain, and then i stumbled upon thisg.data$amt.fac <- as.factor(ifelse(amount <= 2500, "0-2500", ifelse(amount <= 5000, "2600-5000", "5000+")))

if you look at the amount variable, we have the following valueshead(amount)[1] 5 5 5 13 7 11

So when we are actually binning this variable, we are creating a single bin 0-2500 and thats what causes the error.

I changed the bins to 0-5,5-10,10+ and the logistic regression code runs without any problem..

thanks a ton again for all the effort you have put in.. great help !!!!

I had a look at the data again by reloading it and got the following output with the same command as yours

> head(g.data$amount)[1] 1169 5951 2096 7882 4870 9055

I guess there could be an issue where names were assigned to columns. Could you please rerun and check. If you still get an error, please do let me know. I'd be happy to run your script and mine side by side and see where the difference is.

Variable selection is quite a complex game and usually depends on the problem at hand. For simple problem using a stepwise algorithm is sufficient. In R, it can be implemented using the stepAIC() function from the MASS package.

To use the function, first run a logisitic regression using all the variables. Then pass the created object through the stepAIC() function. It will return an object with the selected variables.

In your example, we have many variables here; however not all are working well. If we decide to build a model of, say, 5 variables, what should we do?

I saw that you already recommended to use Stepwise function in R, however how to treat this for the categorical variables. For some of the dummy variables created, we have good significance (p-value), while for others we do not. What would you recommend to do in R?

We can use stepwise to determine which variables are insignificant and this will work for both continuous and dummy variables.

Typically, for a categorical variable with n categories, we put n-1 dummy variables in the regression model to prevent multi-collinearity. The remaining category is then part of the intercept. There are cases where the coefficients of some of the n-1 dummy variables come out to insignificant. We can remove these variables from the model and they will automatically be part of the intercept as well.

That said, stepwise is just one way of selecting features. There are many other ways like regularized regression, genetic algorithms, importance from decision trees to select the best features. Most of these are easily implementable in R.

1) As per my experience with R, the StepAIC() function eliminates (or adds) variables step by step; however does not specify the categories per each variable. For example, if you have age.fac:0-20, age.fac:20-40 as dummy variables (and 40+ as reference variable), it simply eliminates age.fac in one of the stages if it doesn't fit in the model. Do you think such elimination is done only when none of the dummy variables (categories) work in the model?

2) You mentioned there are other ways for variable selection. It would be great if you can mention how to do these in R (or maybe in another post).

Hello, I used your code on a different data set, and it mostly works. However, when I run the plot(performance(m1.scores, "tpr", "fpr"), col = "red"), I get an error: "Error in approxfun(x.values.1, y.values.1, method = "constant", f = 1, : zero non-NA points"Do you know what this means and how to fix it?