Another Data Mining Blog

Useful R Snippets

2012-01-01T19:04:00.000-08:00

Every time I need to do something in R it nearly always means I have to do a Google search or trawl my previous code to see how I did it before. Here I am going to post some snippets of code - mainly for my own use so that I know where to find them. Much of this code will be 'borrowed' and probably not the most efficient (I like to write code the long way so I can follow what is going on) - but it seems to work. If anyone finds it doesn't work or there is a more efficient way then please let me know.

1. Randomly sampling data into a train and test set
totalrecords <- nrow(mydata)
trainfraction = 0.7
trainrecords = as.integer(totalrecords * trainfraction)
allrows <- 1:totalrecords
trainrows <- sample(totalrecords,trainrecords)
testrows <- allrows[-trainrows]
#check
length(trainrows)
length(testrows)
#then build model, something like...
model <- lm(theFormula, data=mydata[trainrows,])

1a. Randomly sampling data into a train and test set
Thanks to Isamoor
trainrows <- runif(nrow(mydata)) > 0.7
testrows <- !trainrows

2. Repeated n fold cross validation

This is to generate a cross-validation set, useful when wanting to know the expected error or for generating a set to use for getting ensemble weightings.

######################################
# the error function
calc_RMSE <- function(act,pred){
    
    aact <- as.matrix(act)
    ppred <- as.matrix(pred)
    
    if(nrow(aact) == nrow(ppred)){ 
    return (sqrt(sum(((ppred) - (aact)) ^ 2) / nrow(aact)))
    } else {
    return (-99)
    }

}
#####################################


###########################
#Load and prepare data
###########################
databuild <- iris
datascore <- iris #put real score set here

#target - what we are predicting
theTarget <- 'Sepal.Length'

#set the formula
theFormula <- as.formula(paste(theTarget," ~ . "))

#find the position of the target
targindex <-  which(names(databuild)==theTarget)

#actuals
build_actuals <- databuild[,targindex]


#######################################
#vectors to score the model outputs
#######################################
buildcases <- nrow(databuild)
scorecases <- nrow(datascore)

pred_train <- vector(length=buildcases)
pred_test <- vector(length=buildcases)
pred_score <- vector(length=scorecases)

pred_trainLoop <- vector(length=buildcases)
pred_testLoop <- vector(length=buildcases)
pred_scoreLoop <- vector(length=scorecases)

#settings
numloops <- 300
numfolds <- 10

test_errors <- vector(length=numloops)
train_errors <- vector(length=numloops)

pred_testLoop <- 0
pred_trainLoop <- 0
pred_scoreLoop <- 0
    
modtype = 'linear regression'

#####################################
# now the work
#####################################        
for(loop in 1:numloops){

    # generate the indicies for each fold    
    id <- sample(rep(seq_len(numfolds), length.out=buildcases))

    # lapply over them:
    indicies <- lapply(seq_len(numfolds), function(a) list(
        test = which(id==a),
        train = which(id!=a)
    ))
    
    #reset the predictions for this loop
    pred_train <- 0
    pred_test <- 0
    pred_score <- 0
    
        for(fold in 1:numfolds){
            
            #set the cases for this fold
            rows_train <- indicies[[fold]]$train
            rows_test  <- indicies[[fold]]$test

            #build the models - use any model
            model <- lm(theFormula, data=databuild[rows_train,])

            #score up the model
            buildPred <- predict(model, databuild, type="response")
            scorepred <- predict(model, datascore, type="response")
            
            #now score the cv and scoring predictions
            z <- buildPred
            z[rows_test] <- 0
            pred_train <- pred_train + z
            pred_test[rows_test] <- buildPred[rows_test]
            pred_score <- pred_score + scorepred
        } #next fold

    
    #average the predictions on the train set
    pred_train <- pred_train / (numfolds - 1)
    pred_score <- pred_score / numfolds
    
    #add to previous loop results
    pred_trainLoop <- pred_trainLoop + pred_train
    pred_testLoop <- pred_testLoop + pred_test
    pred_scoreLoop <- pred_scoreLoop + pred_score
    
    #calculate the errors    
    train_errors[loop] <- calc_RMSE(build_actuals,pred_trainLoop / loop)
    test_errors[loop] <- calc_RMSE(build_actuals,pred_testLoop / loop)
    
    #report
    cat("\nloop = ",loop,"train error = ",train_errors[loop],"cv error = ",test_errors[loop]) 

    #plot a chart as we go
    if(loop>1){
         plot(test_errors[1:loop],col='blue',type='l',main = paste(modtype,numloops,'by',numfolds,'-fold cross validation'), xlab = 'Repetitions', ylab = 'RMSE',ylim = range(rbind(test_errors[1:loop],train_errors[1:loop])))
        abline(h=test_errors[loop],col='blue')
        points(train_errors[1:loop],type='l',col='red')
        abline(h=train_errors[loop],col='red')
        legend('top',c('test','train'),col=c('blue','red'),lty=1)
    }

} #loop
########################


    #the cross validation predictions and scoring set predictions
     #this is what we are after
    cvPredictions <- pred_testLoop / numloops
    scPredictions <- pred_scoreLoop / numloops

    #plot should show decreasing test error with increasing train error
    plot(train_errors,test_errors,type='p')

$3 million - hmmm...

2011-12-26T13:39:00.000-08:00

I've shown this plot a couple of times already, but here it is again, with another that shows the bar we have to limbo under to win $3 million.

Anyone think this is possible?

Any predictions of what the final winning error will be?

I predict 0.453

Place your bets here...

What's Going On Here

2011-12-17T13:54:00.000-08:00

In many of the analytics problems I have been involved in, the problem you end up dealing with is not the one you initially were briefed to solve.

These new problems are always discovered by visualising the data in some way and spotting curious patterns.

Here are a three of examples...

1. Algorithmic Trading Challenge

The Algorithmic Trading Challenge is based on data from the London Stock Exchange and is about things called 'Liquidity Shocks'. I know nothing about these but we had data, so the first thing I did was plot a few graphs to see id I could get a better understanding of things.

The plot below shows the times these 'Liquidity Shocks' occur.

Now it is quite clear there is something going on at 1pm, 2:30pm, after 3:30pm and at 4pm.

Interestingly these spikes are only evident when all commodities are looked at together, they are not as obvious in any individual commodity.

My first question if I was solving a business problem would be to return to the business to get more insight in what was going on here. My initial thoughts were lunch breaks and the opening times of other Stock Exchanges around the world - as 3:30pm London time could be around opening time in New York.

Understanding the cause of these peaks is important as you would expect the reaction to them (the problem to solve) to be a function of the cause.

If we did discover it was the opening times of other exchanges, then I would ask for extra information like the specific dates, so I could calculate when these peaks would occur in the future when the clocks changed. We do not have this information at the current time, or even the day of the week (it can be inferred but not accurately as there will be public holidays when the exchanges are closed)

As it stands any models built could potentially fail on the leaderboard (or real life) data as our model might think 2:30pm is a special time, wheras really it is when another exchange opens, or when people come back from lunch. We need this causal information rather than just dealing with the effect - time differences change - lunch breaks may change.

The current competition data is potentially lacking the full information required to build a model that is as robust as possible over time.

2. Interesting Distributions

One of the first things I do when receiving a data set is to scan the distributions of all variables to sanity check them for anything that looks out of place - but still things can sneak past you.

The following is exam mark data in the range 0-100. If we bin it in 20 bins then things look reasonable, but if we zoom in then we get the 'what is going on here' question again. It is quite clear what is going on, but if exam marks is the thing we are trying to predict, how do we deal with this phenomenon and how would our algorithm cope looking at it blindly? And what if the pass mark changed or rules changes - the algorithm would fail. Again, we need to be aware of the underlying root cause and not just the effect.

3. Don't Get Kicked

This is another Kaggle Competition...

Kicked cars often result when there are tampered odometers, mechanical issues the dealer is not able to address, issues with getting the vehicle title from the seller, or some other unforeseen problem. Kick cars can be very costly to dealers after transportation cost, throw-away repair work, and market losses in reselling the vehicle.

Modelers who can figure out which cars have a higher risk of being kick can provide real value to dealerships trying to provide the best inventory selection possible to their customers.

The challenge of this competition is to predict if the car purchased at the Auction is a Kick (bad buy)

This is a binary classification task and a quick way to spot data issues with this type of problem is to throw it in a decision tree in order to spot what are called 'gimmees'. These are cases that are easily perfectly predictable and are more than often a result of giving prediction data that just shouldn't be there as it is not known at the time (future information) - an extraction issue that would result in a useless model (It is common that people think they have built really good predictive models using future information without really questioning why their models are so good!).

Another reason 'gimmees' occur are poorly defined target variables, that is not excluding certain cases (and example in target marketing would be not excluding dead people from your mailing list and then predicting they won't respond to your offer!)

After a bit of data prep I threw the Don't Get Kicked Data into a Tiberius Decision Tree - the visual below immediately tells me there are clear cut cases of cars that will be kicked - it is almost black and white.

These 'gimmees' can be described by the rules...

[WheelTypeID] = 'NULL' AND [Auction] <> 'MANHEIM'

MANHEIM is an auctioneering company where cars are auctioned - there are 2 main auctioneers in the data set plus 'other'.

Having worked extensively with car auction data before I know that there are certain auctions where only 'write off' cars are sold, that is those that are sold for scrap because they have been in accidents. I also know that different auction houses will record data differently.

The above simple rule easily identifies cars that are more than likely going to be 'knocked' - but this is probably because they are 'knocked' in the first place (are we saying that someone in a coma is more likely to die). Is this useful? Is this a poorly defined definition of what is 'knocked'? Why does a missing value for WheelTypeID make such a big difference between auction houses?

A bit more digging reveals location and the specific buyer drills down on these gimmees even more...

[WheelTypeID] = 'NULL' AND [Auction] <> 'MANHEIM' AND [VNST] in ('NC','AZ') AND [BYRNO] NOT IN (99750,99761)

and after excluding these 'gimmees' it becomes clear there are certain buyers that just don't but knocked cars, especially 99750 and 99761...

byrno = 99750 and VNST in ('SC','NC','UT','ID','PA','WV','MO','WA')
byrno = 99761 and Auction = 'MANHEIM'
byrno = 99761 and MAKE = 'SUZUKI'
byrno = 99761 and SIZE = 'VAN'
byrno = 99761 and VNST IN ('FL','VA')

Now is this actually useful?

The challenge of this competition is to predict if the car purchased at the Auction is a Kick (bad buy)

The model is going to focus on who bought the car rather than the characteristics of the car itself. What happens if buyers suddenly change their policy? Wouldn't we rather just go and speak to these buyers to understand what their policy is and hence get some business understanding? Why is specific auction house location so important? Is it because of the specific auction house itself or that specific cars are actually routed to specific places (this does happen).

Basically if this was a real client engagement I would be going back to them with a lot of questions to help me understand the data better so it can be used in a way that is going to be useful to them.

In Summary

When doing predictive modelling, you can throw the latest hot algorithm at a problem such as a GBM, Neural Net or Random Forest and get impressive results, but unless you thoroughly understand and account for the real dynamics of what is going on then the models could disastrously fail when these dynamics change. I find visualisation the key to spotting and interpreting these dynamics - which is why I would rather have a good data miner who knows what he is doing using free software over a poor data miner with the most expensive software - see http://analystfirst.com/analyst-first-101/

Two Become One

2011-12-16T01:58:00.000-08:00

In the previous post I looked at the HHP leaderboard and discovered some interesting patterns regarding certain teams.

It looks as the evidence proved out to be true, with SD_John and Lily now all of a sudden merging into a single team.

Interestingly they have also been in other competitions with very similar results.

This was the final standing in the Give Me Some Credit competiton,

What is actually more interesting here is the demonstration of overfitting to the leaderboard. Opera Solutions & JYL are more than likely working together and we know Lily & SD_John are working together. If you look at the leaderboard just before the competition ended (on the 30%) you will see Opera near the top but the final position on the 70% was much worse. Similarly a few others found that relying on the leaderboard as an indication of the final position can be misplaced trust.

If you followed the competition forum, you will see team VSU also had multiple accounts for the same person, and they seem to have also fallen into the same trap of overfitting to the leaderboard - they ended up 9th on the 70% when they were first on the 30%.

The data mining lesson here is that you need to take all necessary steps to avoid overfitting, rather than just relying on the leaderboard feedback.

Congratulations to Nathaniel, Eu Jin (small world - I used to work with Nathaniel at the National Australia Bank and regularly see Eu Jin at the Melbourne R user group) and Alec, who clearly did not overfit. A Perfect Storm!

Phantom of the Opera

2011-12-14T13:01:00.000-08:00

There have been some recent announcements on Kaggle reminding competitors about the rules regarding teams and that a single person can't have muliple accounts in order to get around the daily submission limit.

I used the HHP leaderboard as an interesting data source to educate myself on the data manipulation capabilities in R and it became very evident that there was some curious behaviour going on.

From a data scientist viewpoint, this demonstrates the power of the human eye in picking up things that will give you the insight that an algorithm won't. In most (probably all) of my professional projects the important data issues and findings have been a result of looking at visualisations of the data and asking the question "what's going on here!".

The first curiosity on the leaderboard was by trying to discover if the competition was attracting new entrants by looking at the dates of the first submissions of entrants. The two plots below show different ways of looking at the same data. What is obvious is that the 29th Nov had an unusual number of new entrants.

and looked at in another way...

What's going on here?

If you look at the team name of the entrants it is clear that all these accounts are somewhat connected - so no real mystery as to the cause of the blip for this date.

"accnt002" "accnt003" "accnt004" "accnt005" "accnt006"
"accnt007" "accnt008" "accnt009" "cyclops" "Faber"
"Farbe" "Fortis" "glad5" "glad55" "gladiator"
"gladiator1" "gladiator2" "gladiator3" "jackie" "Kaggleacctk"
"KaggleK2" "sashik"

The next two plots show the scores of the first submission of teams.

What's going on here?

The common scores where the steps are seen are the all zeros benchmark, optimised constant benchmark and the code we posted in our writeup - so this is explained. There is another common first score which is another very simple model that many teams independently thought of.

What does raise an eyebrow from the cumulative plot is one team stands out as having a very impressive first score. This is team YARONG who posted a very impressive model of 0.457698 on the first attempt and it still remains their best score 22 attempts later. This is possible (you don't need to submit models to blend them if you have your own holdout set - see the IBM writeup in the KDD Cup Orange Challenge) but somewhat unlikely as we know from the writeups that an individual model will get you no where near this score.

If you look at the dates teams submit and look at some sort of correlation of entry dates, one team appears twice towards the top - SD_John, and they are also at the top of the leaderboard.

td.row td.col pairs correl
UCI-CS273A-RegAll Alex_Tot 27 0.9979902
rutgers HappyAcura 29 0.9978254
SD_John lily 34 0.9974190
Roger99 Krakozjabra 21 0.9956643
SD_John JYL 24 0.9950884
The_Cuckoo's_Nest NumberNinja 23 0.9931073
NumberNinja Chris_R 29 0.9924864

What's is going on here?

If you plot the submissions and scores you will see SD_John and Lily seem to perfectly track each other in both the days they submit, the times they submit and the scores they get.

And on one particular day they get exactly the same score within 5 minutes of each other...

SD_John and JYL seem to also track each other in submission dates. Interestingly JYL has a very similar profile to a member of Opera, and a little digging would suggest this is one and the same person.

So here we can hypothesize that SD_john, lily, JYL and Opera (and evidence also suggests many more teams) are collaborating in some way.

Interesting - all from following your nose, which is what good data mining is all about.

In conclusion, the top of the leaderboard is not really what it appears to be - which I hope will encourage others to keep trying.

The main reason for this investigation was to help me discover what R can do to manipulate data - and the answer is basically anything you want it to do. You first have to know what you want to achieve then do some Googling and you will find some code to help you somewhere.

The Pack is Catching Up

2011-11-17T04:18:00.000-08:00

If you have been keeping an eye on the leaderboard, then you will notice there has been apparently little activity since the milestone 1 deadline. On some occasions there has been no change in the top 40 positions for over one week. This is quite an eerie silence and I suspect there may be a lack of submissions so teams can merge (the total submissions of teams merging has to be less than the number of days the comp has been running).

There have also been some very interesting movements if you look closer - more to come on this in a later post.

If you look further down though, the pack is catching up. In the past six weeks more teams are heading towards the 0.461 mark, which is the point the early leaders got to straight away and was the score to beat. Now it is only a top 50 place.

This score is a good single model. To improve dramatically from there though it is probably necessary to ensemble various models. What is pretty clear though is that the benchmark of 0.40 for the 3 million is impossible (hopefully this might be adjusted?).

I tried to put some nice colours in the chart below, which is generated in R, but could not find any up to date listing of colour codes in R. This is one of the disadvantages of the open source movement - documentation is very low on the contributors list of priorities (and what documentation there is leads a lot to be desired if R is to be used by 'regular' types of people).

I did find the following link though, which is where I got the colours for the plot,

http://colorbrewer2.org/

click on the image to enlarge

If you can't measure it, you can't manage it!

2011-10-07T15:46:00.000-07:00

"If you can't measure it, you can't manage it!" - this is often quoted in the marketing and analytics world.

Apparently

"It comes from a remark attributed to GALILEO, who said 'Count what is
countable, measure what is measurable. What is not measurable, make
measurable'."

see here for more details on the origins of the quote.

Anyway, the point is that if you don't measure the before and after, you have no idea if the actions made a difference.

Last week in my previous blog post, I made available some code that would result in a score of around 0.4635 on the HHP leaderboard. The hope was that others would be able to get this code to work and re-ignite interest in the HHP and solving the worlds health problems.

When I looked at the leaderboard today, it looked on the face of it that this was happening, with a few high risers around the 0.4635 mark.

This was the after, but the before is not shown, so on this evidence there is a bit of a hole in the claim that the code was responsible for this change - we are only showing the part of the leaderboard that we want to - such movement could be going on all over the place.

Good old Jeff Moser at Kaggle has been hard at work making the leaderboard dynamic, so we can actually go back in time. Knowing this I modified a previous R script I posted to plot the leaderboard at 3 points in time, the day I released the code, and a week before and after (the new script can be found at the end of this post).

This resulted in the following 2 plots...

These are now a lot more convincing that the code did make a difference. If you follow Eu Jin Lok on the first plot, you will see he has improved his score but actually gone backwards in rank over the two week period. The second plot shows that the number of submissions to score 0.4635 has dropped.

Anyway, the point of this post is that data scientist have the power to make massive differences to the bottom lines of companies. I have been there and done it, but failed to get too excited as I am unassuming and knew exactly what the result would be because my models told me in advance. If analytics - and us analysts - are to be taken more seriously and given the recognition we deserve, we need to start banging our own drum more than we do. Those who get the recognition and rewards are not always those who do the best work, but those who shout the loudest!

And here is the R script to generate those plots...

##########################################
# Some R function to plot the Kaggle 
# leaderboard at different points in time
##########################################

plotKaggleLeaderboard <- function(theURL
,theDates
,myTeam
,topX=100
,title
,plottype=1){
    
#this library does all the hard work
#for windows users to install this packages see
#http://cran.r-project.org/bin/windows/contrib/r-release/ReadMe
#and http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.13/
library(XML)

theColours <- c(552,254,26)

for (i in 1:length(compDates))
{

#go and read the tables from the web page
thisDate <- paste(theURL,"?asOf=",compDates[i],sep="")
tables <- readHTMLTable(thisDate)

#get the table of interest
#names(tables)
lb <- tables[['leaderboard-table']]

#see what the columns are
#colnames(lb)
#nrow(lb)

#convert to numeric - see ?readHTMLTable
#numeric_cols <- c('AUC','Entries')
#numeric_cols <- c('RMSLE','Entries')
numeric_cols <- c(4,5)
lb[numeric_cols] = lapply(lb[numeric_cols], function(x) as.numeric(gsub(".* ", "", as.character(x))))


#extract the team name from a messy field
team_col <- c('Team Name')
#lb[team_col]

#split the field by "\r\n" than denotes the break between fields within a field 
lb[,team_col] <- sapply(strsplit(as.character(lb[,team_col]), "[\r\n]"), function (x) x[1]) 


#####################
#      plot
#####################
myRank <- which(lb[team_col] == myTeam)
myText = paste("following team",myTeam)

error_data <- lb[,4]
entries_data <- lb[,5]

if(plottype==1){
if(i==1)  plot(error_data[1:topX],col = theColours[i],type="l",xlab='Rank',ylab='Error',main = title, sub=myText)
if(i>1) points(error_data[1:topX],col = theColours[i],type="l")
          
#mark position
points(myRank,error_data[myRank],col=theColours[i],pch=19)
}


if(plottype==2){
if(i==1){ plot(error_data[1:topX],entries_data[1:topX],type = 'p'
            ,main = title
            , xlab = 'Error'
            , ylab = 'Entries'
            , col= theColours[i]
    )}

if(i>1) points(error_data[1:topX],entries_data[1:topX],type = 'p', col= theColours[i])

lines(lowess(error_data[1:topX],entries_data[1:topX]), col=theColours[i], lwd=2.5) # lowess line (x,y)

#mark my position
points(error_data[myRank],entries_data[myRank],col=theColours[i],pch=19)
}

}

legend("topleft",legend=theDates,text.col=theColours)

}
###########################
# End of Function
# plotKaggleLeaderboard
###########################



compURL <- "http://www.heritagehealthprize.com/c/hhp/Leaderboard"
compDates <- c('2011-09-24','2011-10-01','2011-10-07')
compTeam <- 'Eu Jin Lok'
compTopX <- 150
compTitle <- 'HHP Leaderboard Through Time'

plotKaggleLeaderboard(theURL = compURL
            ,theDates = compDates
                      ,myTeam = compTeam
                      ,topX = compTopX
                      ,title = compTitle
                       ,plottype=1
                      )

Code for a Respectable HHP Model

2011-10-01T06:43:00.000-07:00

Our team Market Makers did well in the HHP first milestone prize, taking out top spot just above Willem Mestrom. Congratulations to all who took part.

For those of you who haven't seen the presentation, there is a recording at http://www.youtube.com/watch?v=xziMY3nDB_E

Reports were written on our methods and are available here...
http://www.heritagehealthprize.com/c/hhp/Leaderboard/milestone1

Other competitors now have 30 days to read the reports and come back with any concerns.

In our report we included some source code, which is a bit of a pain to copy and paste, so I have reproduced it here. There are two bits of code to run, and at the end a file ready for submission should appear.

The first bit is the data preparation which I did in SQL. I used SQL Server, but it should be easy to run it in any SQL database without much messing around. This took about 8 minutes to run, and at the end is a table ready for use by the second part, which is the modelling performed in R, which takes about 15 minutes.

There are two modifications to the code to make - in the SQL you need to set the correct path to the source files - and in R you need to put in the correct server name so the data can be read directly from the database.

So hopefully in two mouse clicks there will be a file sitting there ready for submission.

If you are not sure about using SQL Server, then the earlier blog posts should point you in the right direction...

http://anotherdataminingblog.blogspot.com/2011/05/progress-loading-hhp-data.html

Note that the SQL provided does not generate all the variables we used - others were created, all using SQL, but the code was a little more complicated. The extra variables are described in our write up. Also the parameters used in the GBM modelling are not the most efficient. You will need more trees and a lower shrinkage factor to tease out a better model (see the R documentaion) - but this one runs in a reasonable time.

I am keen to see if this code is useful to anybody - if you can't get it to work then please post a comment hear. If you do use it then I would also like to hear from you!

SQL Code

/**************************************************************
* SQL Code to create an example data set for the HHP
*
* Edit the path in the 'bulk insert' commands to locate
* the source data
* The end result is a table called 'modelling_set' which can
* then be used to build predictive models
*
* created in SQL server express
* http://www.microsoft.com/sqlserver/en/us/editions/express.aspx
*****************************************************************/

/**************************
create a new database
**************************/
CREATE DATABASE HHP_comp
GO
USE HHP_comp

/**************************
load in the raw data
**************************/

--claims
CREATE TABLE Claims
(
MemberID VARCHAR(8) --integers starting with 0, could be text!
, ProviderID VARCHAR(7) --integers starting with 0, could be text!
, Vendor VARCHAR(6) --integers starting with 0, could be text!
, PCP VARCHAR(5) --integers starting with 0, could be text!
, Year VARCHAR(2)
, Specialty VARCHAR(25)
, PlaceSvc VARCHAR(19)
, PayDelay VARCHAR(4)
, LengthOfStay VARCHAR(10)
, DSFS VARCHAR(12)
, PrimaryConditionGroup VARCHAR(8)
, CharlsonIndex VARCHAR(3)
, ProcedureGroup VARCHAR(4)
, SupLOS TINYINT
)

BULK INSERT Claims
FROM 'F:\comps\hhp\raw data\HHP_release2\Claims.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

--members
CREATE TABLE Members
(
MemberID_M VARCHAR(8) --integers starting with 0, could be text!
, AgeAtFirstClaim VARCHAR(5)
, Sex VARCHAR(1)
)

BULK INSERT Members
FROM 'F:\comps\hhp\raw data\HHP_release2\Members.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

-- drug count
CREATE TABLE DrugCount
(
MemberID INT
, Year VARCHAR(2)
, DSFS VARCHAR(12)
, DrugCount VARCHAR(2)
)

BULK INSERT DrugCount
FROM 'F:\comps\hhp\raw data\HHP_release3\DrugCount.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

-- Lab Count
CREATE TABLE LabCount
(
MemberID INT
, Year VARCHAR(2)
, DSFS VARCHAR(12)
, LabCount VARCHAR(3)
)

BULK INSERT LabCount
FROM 'F:\comps\hhp\raw data\HHP_release3\LabCount.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

--DaysInHospital_Y2
CREATE TABLE DaysInHospital_Y2
(
MemberID INT
, ClaimsTruncated TINYINT
, DaysInHospital TINYINT
)

BULK INSERT DaysInHospital_Y2
FROM 'F:\comps\hhp\raw data\HHP_release2\DaysInHospital_Y2.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

-- DaysInHospital_Y3
CREATE TABLE DaysInHospital_Y3
(
MemberID INT
, ClaimsTruncated TINYINT
, DaysInHospital TINYINT
)

BULK INSERT DaysInHospital_Y3
FROM 'F:\comps\hhp\raw data\HHP_release2\DaysInHospital_Y3.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

-- Target
CREATE TABLE Target
(
MemberID INT
, ClaimsTruncated TINYINT
, DaysInHospital TINYINT
)

BULK INSERT Target
FROM 'F:\comps\hhp\raw data\HHP_release2\Target.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

/*************************
adjust the claims data to
convert text to integers
**************************/

-- PayDelay
ALTER TABLE Claims
ADD PayDelayI integer
GO

UPDATE Claims
SET PayDelayI = CASE WHEN PayDelay = '162+' THEN 162 ELSE CAST(PayDelay AS integer) END

--dsfs
ALTER TABLE Claims
ADD dsfsI integer
GO

UPDATE Claims
SET dsfsI =
CASE
WHEN dsfs = '0- 1 month' THEN 1
WHEN dsfs = '1- 2 months' THEN 2
WHEN dsfs = '2- 3 months' THEN 3
WHEN dsfs = '3- 4 months' THEN 4
WHEN dsfs = '4- 5 months' THEN 5
WHEN dsfs = '5- 6 months' THEN 6
WHEN dsfs = '6- 7 months' THEN 7
WHEN dsfs = '7- 8 months' THEN 8
WHEN dsfs = '8- 9 months' THEN 9
WHEN dsfs = '9-10 months' THEN 10
WHEN dsfs = '10-11 months' THEN 11
WHEN dsfs = '11-12 months' THEN 12
WHEN dsfs IS NULL THEN NULL
END

-- CharlsonIndex
ALTER TABLE Claims
ADD CharlsonIndexI INTEGER
GO

UPDATE Claims
SET CharlsonIndexI =
CASE
WHEN CharlsonIndex = '0' THEN 0
WHEN CharlsonIndex = '1-2' THEN 2
WHEN CharlsonIndex = '3-4' THEN 4
WHEN CharlsonIndex = '5+' THEN 6
END

-- LengthOfStay
ALTER TABLE Claims
ADD LengthOfStayI INTEGER
GO

UPDATE Claims
SET LengthOfStayI =
CASE
WHEN LengthOfStay = '1 day' THEN 1
WHEN LengthOfStay = '2 days' THEN 2
WHEN LengthOfStay = '3 days' THEN 3
WHEN LengthOfStay = '4 days' THEN 4
WHEN LengthOfStay = '5 days' THEN 5
WHEN LengthOfStay = '6 days' THEN 6
WHEN LengthOfStay = '1- 2 weeks' THEN 11
WHEN LengthOfStay = '2- 4 weeks' THEN 21
WHEN LengthOfStay = '4- 8 weeks' THEN 42
WHEN LengthOfStay = '26+ weeks' THEN 180
WHEN LengthOfStay IS NULL THEN null
END

/**************************
create a summary table
at the member/year level
***************************/
SELECT
year
,Memberid

,COUNT(*) AS no_Claims
,COUNT(DISTINCT ProviderID) AS no_Providers
,COUNT(DISTINCT Vendor) AS no_Vendors
,COUNT(DISTINCT PCP) AS no_PCPs
,COUNT(DISTINCT PlaceSvc) AS no_PlaceSvcs
,COUNT(DISTINCT Specialty) AS no_Specialities
,COUNT(DISTINCT PrimaryConditionGroup) AS no_PrimaryConditionGroups
,COUNT(DISTINCT ProcedureGroup) AS no_ProcedureGroups

,MAX(PayDelayI) AS PayDelay_max
,MIN(PayDelayI) AS PayDelay_min
,AVG(PayDelayI) AS PayDelay_ave
,(CASE WHEN COUNT(*) = 1 THEN 0 ELSE STDEV(PayDelayI) END) AS PayDelay_stdev

,MAX(LengthOfStayI) AS LOS_max
,MIN(LengthOfStayI) AS LOS_min
,AVG(LengthOfStayI) AS LOS_ave
,(CASE WHEN COUNT(*) = 1 THEN 0 ELSE STDEV(LengthOfStayI) END) AS LOS_stdev

,SUM(CASE WHEN LENGTHOFSTAY IS NULL AND SUPLOS = 0 THEN 1 ELSE 0 END) AS LOS_TOT_UNKNOWN
,SUM(CASE WHEN LENGTHOFSTAY IS NULL AND SUPLOS = 1 THEN 1 ELSE 0 END) AS LOS_TOT_SUPRESSED
,SUM(CASE WHEN LENGTHOFSTAY IS NOT NULL THEN 1 ELSE 0 END) AS LOS_TOT_KNOWN

,MAX(dsfsI) AS dsfs_max
,MIN(dsfsI) AS dsfs_min
,MAX(dsfsI) - MIN(dsfsI) AS dsfs_range
,AVG(dsfsI) AS dsfs_ave
,(CASE WHEN COUNT(*) = 1 THEN 0 ELSE STDEV(dsfsI) END) AS dsfs_stdev

,MAX(CharlsonIndexI) AS CharlsonIndexI_max
,MIN(CharlsonIndexI) AS CharlsonIndexI_min
,AVG(CharlsonIndexI) AS CharlsonIndexI_ave
,MAX(CharlsonIndexI) - MIN(CharlsonIndexI) AS CharlsonIndexI_range
,(CASE WHEN COUNT(*) = 1 THEN 0 ELSE STDEV(CharlsonIndexI) END) AS CharlsonIndexI_stdev

,SUM(CASE WHEN PrimaryConditionGroup = 'MSC2a3' THEN 1 ELSE 0 END) AS pcg1
,SUM(CASE WHEN PrimaryConditionGroup = 'METAB3' THEN 1 ELSE 0 END) AS pcg2
,SUM(CASE WHEN PrimaryConditionGroup = 'ARTHSPIN' THEN 1 ELSE 0 END) AS pcg3
,SUM(CASE WHEN PrimaryConditionGroup = 'NEUMENT' THEN 1 ELSE 0 END) AS pcg4
,SUM(CASE WHEN PrimaryConditionGroup = 'RESPR4' THEN 1 ELSE 0 END) AS pcg5
,SUM(CASE WHEN PrimaryConditionGroup = 'MISCHRT' THEN 1 ELSE 0 END) AS pcg6
,SUM(CASE WHEN PrimaryConditionGroup = 'SKNAUT' THEN 1 ELSE 0 END) AS pcg7
,SUM(CASE WHEN PrimaryConditionGroup = 'GIBLEED' THEN 1 ELSE 0 END) AS pcg8
,SUM(CASE WHEN PrimaryConditionGroup = 'INFEC4' THEN 1 ELSE 0 END) AS pcg9
,SUM(CASE WHEN PrimaryConditionGroup = 'TRAUMA' THEN 1 ELSE 0 END) AS pcg10
,SUM(CASE WHEN PrimaryConditionGroup = 'HEART2' THEN 1 ELSE 0 END) AS pcg11
,SUM(CASE WHEN PrimaryConditionGroup = 'RENAL3' THEN 1 ELSE 0 END) AS pcg12
,SUM(CASE WHEN PrimaryConditionGroup = 'ROAMI' THEN 1 ELSE 0 END) AS pcg13
,SUM(CASE WHEN PrimaryConditionGroup = 'MISCL5' THEN 1 ELSE 0 END) AS pcg14
,SUM(CASE WHEN PrimaryConditionGroup = 'ODaBNCA' THEN 1 ELSE 0 END) AS pcg15
,SUM(CASE WHEN PrimaryConditionGroup = 'UTI' THEN 1 ELSE 0 END) AS pcg16
,SUM(CASE WHEN PrimaryConditionGroup = 'COPD' THEN 1 ELSE 0 END) AS pcg17
,SUM(CASE WHEN PrimaryConditionGroup = 'GYNEC1' THEN 1 ELSE 0 END) AS pcg18
,SUM(CASE WHEN PrimaryConditionGroup = 'CANCRB' THEN 1 ELSE 0 END) AS pcg19
,SUM(CASE WHEN PrimaryConditionGroup = 'FXDISLC' THEN 1 ELSE 0 END) AS pcg20
,SUM(CASE WHEN PrimaryConditionGroup = 'AMI' THEN 1 ELSE 0 END) AS pcg21
,SUM(CASE WHEN PrimaryConditionGroup = 'PRGNCY' THEN 1 ELSE 0 END) AS pcg22
,SUM(CASE WHEN PrimaryConditionGroup = 'HEMTOL' THEN 1 ELSE 0 END) AS pcg23
,SUM(CASE WHEN PrimaryConditionGroup = 'HEART4' THEN 1 ELSE 0 END) AS pcg24
,SUM(CASE WHEN PrimaryConditionGroup = 'SEIZURE' THEN 1 ELSE 0 END) AS pcg25
,SUM(CASE WHEN PrimaryConditionGroup = 'APPCHOL' THEN 1 ELSE 0 END) AS pcg26
,SUM(CASE WHEN PrimaryConditionGroup = 'CHF' THEN 1 ELSE 0 END) AS pcg27
,SUM(CASE WHEN PrimaryConditionGroup = 'GYNECA' THEN 1 ELSE 0 END) AS pcg28
,SUM(CASE WHEN PrimaryConditionGroup IS NULL THEN 1 ELSE 0 END) AS pcg29
,SUM(CASE WHEN PrimaryConditionGroup = 'PNEUM' THEN 1 ELSE 0 END) AS pcg30
,SUM(CASE WHEN PrimaryConditionGroup = 'RENAL2' THEN 1 ELSE 0 END) AS pcg31
,SUM(CASE WHEN PrimaryConditionGroup = 'GIOBSENT' THEN 1 ELSE 0 END) AS pcg32
,SUM(CASE WHEN PrimaryConditionGroup = 'STROKE' THEN 1 ELSE 0 END) AS pcg33
,SUM(CASE WHEN PrimaryConditionGroup = 'CANCRA' THEN 1 ELSE 0 END) AS pcg34
,SUM(CASE WHEN PrimaryConditionGroup = 'FLaELEC' THEN 1 ELSE 0 END) AS pcg35
,SUM(CASE WHEN PrimaryConditionGroup = 'MISCL1' THEN 1 ELSE 0 END) AS pcg36
,SUM(CASE WHEN PrimaryConditionGroup = 'HIPFX' THEN 1 ELSE 0 END) AS pcg37
,SUM(CASE WHEN PrimaryConditionGroup = 'METAB1' THEN 1 ELSE 0 END) AS pcg38
,SUM(CASE WHEN PrimaryConditionGroup = 'PERVALV' THEN 1 ELSE 0 END) AS pcg39
,SUM(CASE WHEN PrimaryConditionGroup = 'LIVERDZ' THEN 1 ELSE 0 END) AS pcg40
,SUM(CASE WHEN PrimaryConditionGroup = 'CATAST' THEN 1 ELSE 0 END) AS pcg41
,SUM(CASE WHEN PrimaryConditionGroup = 'CANCRM' THEN 1 ELSE 0 END) AS pcg42
,SUM(CASE WHEN PrimaryConditionGroup = 'PERINTL' THEN 1 ELSE 0 END) AS pcg43
,SUM(CASE WHEN PrimaryConditionGroup = 'PNCRDZ' THEN 1 ELSE 0 END) AS pcg44
,SUM(CASE WHEN PrimaryConditionGroup = 'RENAL1' THEN 1 ELSE 0 END) AS pcg45
,SUM(CASE WHEN PrimaryConditionGroup = 'SEPSIS' THEN 1 ELSE 0 END) AS pcg46

,SUM(CASE WHEN Specialty = 'Internal' THEN 1 ELSE 0 END) AS sp1
,SUM(CASE WHEN Specialty = 'Laboratory' THEN 1 ELSE 0 END) AS sp2
,SUM(CASE WHEN Specialty = 'General Practice' THEN 1 ELSE 0 END) AS sp3
,SUM(CASE WHEN Specialty = 'Surgery' THEN 1 ELSE 0 END) AS sp4
,SUM(CASE WHEN Specialty = 'Diagnostic Imaging' THEN 1 ELSE 0 END) AS sp5
,SUM(CASE WHEN Specialty = 'Emergency' THEN 1 ELSE 0 END) AS sp6
,SUM(CASE WHEN Specialty = 'Other' THEN 1 ELSE 0 END) AS sp7
,SUM(CASE WHEN Specialty = 'Pediatrics' THEN 1 ELSE 0 END) AS sp8
,SUM(CASE WHEN Specialty = 'Rehabilitation' THEN 1 ELSE 0 END) AS sp9
,SUM(CASE WHEN Specialty = 'Obstetrics and Gynecology' THEN 1 ELSE 0 END) AS sp10
,SUM(CASE WHEN Specialty = 'Anesthesiology' THEN 1 ELSE 0 END) AS sp11
,SUM(CASE WHEN Specialty = 'Pathology' THEN 1 ELSE 0 END) AS sp12
,SUM(CASE WHEN Specialty IS NULL THEN 1 ELSE 0 END) AS sp13

,SUM(CASE WHEN ProcedureGroup = 'EM' THEN 1 ELSE 0 END ) AS pg1
,SUM(CASE WHEN ProcedureGroup = 'PL' THEN 1 ELSE 0 END ) AS pg2
,SUM(CASE WHEN ProcedureGroup = 'MED' THEN 1 ELSE 0 END ) AS pg3
,SUM(CASE WHEN ProcedureGroup = 'SCS' THEN 1 ELSE 0 END ) AS pg4
,SUM(CASE WHEN ProcedureGroup = 'RAD' THEN 1 ELSE 0 END ) AS pg5
,SUM(CASE WHEN ProcedureGroup = 'SDS' THEN 1 ELSE 0 END ) AS pg6
,SUM(CASE WHEN ProcedureGroup = 'SIS' THEN 1 ELSE 0 END ) AS pg7
,SUM(CASE WHEN ProcedureGroup = 'SMS' THEN 1 ELSE 0 END ) AS pg8
,SUM(CASE WHEN ProcedureGroup = 'ANES' THEN 1 ELSE 0 END ) AS pg9
,SUM(CASE WHEN ProcedureGroup = 'SGS' THEN 1 ELSE 0 END ) AS pg10
,SUM(CASE WHEN ProcedureGroup = 'SEOA' THEN 1 ELSE 0 END ) AS pg11
,SUM(CASE WHEN ProcedureGroup = 'SRS' THEN 1 ELSE 0 END ) AS pg12
,SUM(CASE WHEN ProcedureGroup = 'SNS' THEN 1 ELSE 0 END ) AS pg13
,SUM(CASE WHEN ProcedureGroup = 'SAS' THEN 1 ELSE 0 END ) AS pg14
,SUM(CASE WHEN ProcedureGroup = 'SUS' THEN 1 ELSE 0 END ) AS pg15
,SUM(CASE WHEN ProcedureGroup IS NULL THEN 1 ELSE 0 END ) AS pg16
,SUM(CASE WHEN ProcedureGroup = 'SMCD' THEN 1 ELSE 0 END ) AS pg17
,SUM(CASE WHEN ProcedureGroup = 'SO' THEN 1 ELSE 0 END ) AS pg18

,SUM(CASE WHEN PlaceSvc = 'Office' THEN 1 ELSE 0 END) AS ps1
,SUM(CASE WHEN PlaceSvc = 'Independent Lab' THEN 1 ELSE 0 END) AS ps2
,SUM(CASE WHEN PlaceSvc = 'Urgent Care' THEN 1 ELSE 0 END) AS ps3
,SUM(CASE WHEN PlaceSvc = 'Outpatient Hospital' THEN 1 ELSE 0 END) AS ps4
,SUM(CASE WHEN PlaceSvc = 'Inpatient Hospital' THEN 1 ELSE 0 END) AS ps5
,SUM(CASE WHEN PlaceSvc = 'Ambulance' THEN 1 ELSE 0 END) AS ps6
,SUM(CASE WHEN PlaceSvc = 'Other' THEN 1 ELSE 0 END) AS ps7
,SUM(CASE WHEN PlaceSvc = 'Home' THEN 1 ELSE 0 END) AS ps8
,SUM(CASE WHEN PlaceSvc IS NULL THEN 1 ELSE 0 END) AS ps9

INTO claims_per_member
FROM Claims
GROUP BY year,Memberid

-- remove some nulls
UPDATE claims_per_member
SET LOS_max = 0 WHERE LOS_max IS NULL

UPDATE claims_per_member
SET LOS_min = 0 WHERE LOS_min IS NULL

UPDATE claims_per_member
SET LOS_ave = 0 WHERE LOS_ave IS NULL

UPDATE claims_per_member
SET LOS_stdev = -1 WHERE LOS_stdev IS NULL

UPDATE claims_per_member
SET dsfs_max = 0 WHERE dsfs_max IS NULL

UPDATE claims_per_member
SET dsfs_min = 0 WHERE dsfs_min IS NULL

UPDATE claims_per_member
SET dsfs_ave = 0 WHERE dsfs_ave IS NULL

UPDATE claims_per_member
SET dsfs_stdev = -1 WHERE dsfs_stdev IS NULL

UPDATE claims_per_member
SET dsfs_range = -1 WHERE dsfs_range IS NULL

UPDATE claims_per_member
SET CharlsonIndexI_range = -1 WHERE CharlsonIndexI_range IS NULL

/***********************************
Members
***********************************/

-- create binary flags for age
ALTER TABLE Members ADD age_05 INT
ALTER TABLE Members ADD age_15 INT
ALTER TABLE Members ADD age_25 INT
ALTER TABLE Members ADD age_35 INT
ALTER TABLE Members ADD age_45 INT
ALTER TABLE Members ADD age_55 INT
ALTER TABLE Members ADD age_65 INT
ALTER TABLE Members ADD age_75 INT
ALTER TABLE Members ADD age_85 INT
ALTER TABLE Members ADD age_MISS INT

GO

UPDATE Members SET age_05 = CASE WHEN ageATfirstclaim = '0-9' THEN 1 ELSE 0 END
UPDATE Members SET age_15 = CASE WHEN ageATfirstclaim = '10-19' THEN 1 ELSE 0 END
UPDATE Members SET age_25 = CASE WHEN ageATfirstclaim = '20-29' THEN 1 ELSE 0 END
UPDATE Members SET age_35 = CASE WHEN ageATfirstclaim = '30-39' THEN 1 ELSE 0 END
UPDATE Members SET age_45 = CASE WHEN ageATfirstclaim = '40-49' THEN 1 ELSE 0 END
UPDATE Members SET age_55 = CASE WHEN ageATfirstclaim = '50-59' THEN 1 ELSE 0 END
UPDATE Members SET age_65 = CASE WHEN ageATfirstclaim = '60-69' THEN 1 ELSE 0 END
UPDATE Members SET age_75 = CASE WHEN ageATfirstclaim = '70-79' THEN 1 ELSE 0 END
UPDATE Members SET age_85 = CASE WHEN ageATfirstclaim = '80+' THEN 1 ELSE 0 END
UPDATE Members SET age_MISS = CASE WHEN ageATfirstclaim IS NULL THEN 1 ELSE 0 END

--create binary flags for sex
ALTER TABLE Members
ADD sexMALE INT
GO

UPDATE Members
SET SexMALE =
CASE
WHEN Sex = 'M' THEN 1 ELSE 0
END

ALTER TABLE Members
ADD sexFEMALE INT
GO

UPDATE Members
SET SexFEMALE =
CASE
WHEN Sex = 'F' THEN 1 ELSE 0
END

ALTER TABLE Members
ADD sexMISS INT
GO

UPDATE Members
SET SexMISS =
CASE
WHEN Sex IS NULL THEN 1 ELSE 0
END

/******************
DRUG COUNTS
******************/

-- convert to integers
ALTER TABLE drugcount ADD DrugCountI INT
GO
UPDATE DRUGCOUNT
SET DrugCountI =
CASE WHEN DrugCount = '7+' THEN 7 ELSE DrugCount END

SELECT
memberID AS memberID_dc
,Year AS YEAR_dc
,MAX(drugcountI) AS drugCount_max
,MIN(drugcountI) AS drugCount_min
,AVG(drugcountI * 1.0) AS drugCount_ave
,COUNT(*) AS drugcount_months
INTO DRUGCOUNT_SUMMARY
FROM
drugcount
GROUP BY
memberID
,Year

/******************
LAB COUNTS
******************/

-- convert to integers
ALTER TABLE LabCount ADD LabCountI INT
GO
UPDATE LabCount
SET LabCountI =
CASE WHEN LabCount = '10+' THEN 10 ELSE LabCount END

SELECT
memberID AS memberID_lc
,Year AS YEAR_lc
,MAX(labcountI) AS labCount_max
,MIN(labcountI) AS labCount_min
,AVG(labcountI * 1.0) AS labCount_ave
,COUNT(*) AS labcount_months
INTO LABCOUNT_SUMMARY
FROM
labcount
GROUP BY
memberID
,Year

/********************************
Targets
********************************/

SELECT *
INTO DIH
FROM
(
SELECT
MemberID AS MemberID_t
,'Y1' AS YEAR_t
,ClaimsTruncated
,DaysInHospital
,1 AS trainset
FROM DaysInHospital_Y2

UNION ALL

SELECT
MemberID AS MemberID_t
,'Y2' AS YEAR_t
,ClaimsTruncated
,DaysInHospital
,1 AS trainset
FROM DaysInHospital_Y3

UNION ALL

SELECT
MemberID AS MemberID_t
,'Y3' AS YEAR_t
,ClaimsTruncated
,null AS DaysInHospital
,0 AS trainset
FROM Target
) a

/*****************************
Now merge them all together to
create the modeling data SET
******************************/
SELECT a.*,b.*
INTO #temp1
FROM
DIH a
LEFT OUTER JOIN
members b
on a.MemberID_t = B.Memberid_M

ALTER TABLE #temp1 DROP COLUMN Memberid_M
ALTER TABLE #temp1 DROP COLUMN AgeAtFirstClaim
ALTER TABLE #temp1 DROP COLUMN Sex
GO

SELECT a.*,b.*
INTO #temp2
FROM
#temp1 a
LEFT OUTER JOIN
claims_per_member b
on a.MemberID_t = B.Memberid
AND a.YEAR_t = b.year

ALTER TABLE #temp2 DROP COLUMN Memberid
ALTER TABLE #temp2 DROP COLUMN year
GO

SELECT a.*,b.*
INTO #temp3
FROM
#temp2 a
LEFT OUTER JOIN
DRUGCOUNT_SUMMARY b
on a.MemberID_t = B.Memberid_dc
AND a.YEAR_t = b.YEAR_dc

ALTER TABLE #temp3 DROP COLUMN Memberid_dc
ALTER TABLE #temp3 DROP COLUMN YEAR_dc
GO

SELECT a.*,b.*
INTO #temp4
FROM
#temp3 a
LEFT OUTER JOIN
LABCOUNT_SUMMARY b
on a.MemberID_t = B.Memberid_lc
AND a.YEAR_t = b.YEAR_lc

ALTER TABLE #temp4 DROP COLUMN Memberid_lc
ALTER TABLE #temp4 DROP COLUMN YEAR_lc
GO

-- removel nulls for those who had
-- no lab or drug information
ALTER TABLE #temp4 ADD labNull INT
ALTER TABLE #temp4 ADD drugNull INT
GO

UPDATE #temp4 SET labNull = 0
UPDATE #temp4 SET labNull = 1 WHERE labCount_max IS NULL

UPDATE #temp4 SET drugNull = 0
UPDATE #temp4 SET drugNull = 1 WHERE drugCount_max IS NULL

UPDATE #temp4 SET labCount_max = 0 WHERE labCount_max IS NULL
UPDATE #temp4 SET labCount_min = 0 WHERE labCount_min IS NULL
UPDATE #temp4 SET labCount_ave = 0 WHERE labCount_ave IS NULL
UPDATE #temp4 SET labcount_months = 0 WHERE labcount_months IS NULL

UPDATE #temp4 SET drugCount_max = 0 WHERE drugCount_max IS NULL
UPDATE #temp4 SET drugCount_min = 0 WHERE drugCount_min IS NULL
UPDATE #temp4 SET drugCount_ave = 0 WHERE drugCount_ave IS NULL
UPDATE #temp4 SET drugcount_months = 0 WHERE drugcount_months IS NULL

SELECT *
INTO modelling_set
FROM #temp4

R code

########################################
# Example GBM model for HHP
# scores ~ 0.4635 on leaderboard
# which would be 55th position of 510
# as at 9th Sept 2011
#
# Requires the data having been prepared
# using the SQL supplied
#
########################################

starttime <- proc.time()

########################################
#load the data
########################################
library(RODBC)

#set a connection to the database 
conn <- odbcDriverConnect("driver=SQL Server;database=HHP_comp;server=servername\\SQLEXPRESS;")

#or this method involves setting up a DSN (Data Source Name) called HHP_comp
#conn <- odbcConnect("HHP_comp")

alldata <- sqlQuery(conn,"select * from modelling_set")


########################################
# arrange the data
########################################

#identify train and leaderboard data
trainrows <- which(alldata$trainset == 1)
scorerows <- which(alldata$trainset == 0)

#sanity check the size of each set
length(trainrows)
length(scorerows)

#display the column names
colnames(alldata)

#memberid is required as key for submission set
memberid <- alldata[scorerows,'MemberID_t']

#remove redundant fields
alldata$MemberID_t <- NULL
alldata$YEAR_t <- NULL
alldata$trainset <- NULL

#target - what we are predicting
theTarget <- 'DaysInHospital'

#put the target on the log scale
alldata[trainrows,theTarget] <- log1p(alldata[trainrows,theTarget]) 

#find the position of the target
targindex <-  which(names(alldata)==theTarget)


########################################
# build the model
########################################

#GBM model settings, these can be varied
GBM_NTREES = 500
GBM_SHRINKAGE = 0.05
GBM_DEPTH = 4
GBM_MINOBS = 50

#build the GBM model
library(gbm)
GBM_model <- gbm.fit(
             x = alldata[trainrows,-targindex]
            ,y = alldata[trainrows,targindex]
            ,distribution = "gaussian"
            ,n.trees = GBM_NTREES
            ,shrinkage = GBM_SHRINKAGE
            ,interaction.depth = GBM_DEPTH
            ,n.minobsinnode = GBM_MINOBS
            ,verbose = TRUE) 

#list variable importance
summary(GBM_model,GBM_NTREES)

#predict for the leaderboard data
prediction <- predict.gbm(object = GBM_model
              ,newdata = alldata[scorerows,-targindex]
              ,GBM_NTREES)

#put on correct scale and cap
prediction <- expm1(prediction)
prediction <- pmin(15,prediction)
prediction <- pmax(0,prediction)

#plot the submission distribution
hist(prediction, breaks=500)


########################################
#write the submission to file
########################################
submission <- cbind(memberid,prediction)
colnames(submission) <- c("MemberID","DaysInHospital")
fnname <- "C:\\GBM_demo1.csv"
write.csv(submission, file=fnname, row.names = FALSE)

elapsedtime <- proc.time() - starttime
cat("\nFinished\n",elapsedtime)

First Milestone

2011-08-31T23:30:00.000-07:00

The time has now been and gone for submitting an entry for the first milestone, and the chance to win either $30,000 or $20,000 for the top two places.

The images below are the state of play (on the 30% visible data!) just after the first milestone ( I didn't capture this image in time so some teams have made an extra submission since).

There was a bit of movement in the last 2 days.

Tom & Jerry changed their name to reveal themselves (Petterson & Caetano). This was interesting for me as they are from a relatively new machine learning lab here in Australia that I had not been aware of. They jumped up to 3rd position before the last submission but were overhauled at the last minute
Willem Mestrom made a big leap on his last submission to get into 3rd position. He was obviously keeping his best (or was he?) until the very last moment. This might be a good tactical play if he looking towards the $3 million rather than the short term smaller prizes.
Planet Melbourne, after making consistent daily submissions throughout the competition did not take the opportunity to make one last submission.
It is good to see Chris R still in the top 10. Chris has been up there since the very start.Interestingly Chris was cut off from submitting for a while due to Hurricane Irene taking his power away.

Interestingly, positions 2,3 & 4 are now taken by members of The Ensemble from Netflix fame. Position 1 contains Dave, who was in an original team called Ensemble Team, and also entered Netflix as a member of Ensemble Experts , and I (also a member of Market Makers) ran a competition on Ensembling the Netflix Data and have looked at ensembling other competitions.
There seems to be a pattern here!

Note you can click on the images to enlarge them.

The plot below shows the Lowess line fit of the error as a function of the number of entries submitted.

First Milestone Just Around the Corner

2011-08-27T22:29:00.000-07:00

There is about 4 days to go until the deadline for the first milestone, so I thought it would be interesting to capture the leaderboard now and then in 4 days, to see if things change wildly. At this stage it can be like a game of poker.

All the top 10 teams are still in the hunt and making daily submissions. Tom & Jerry were late starters but have made astonishing progress to 3rd place. They currently split 3 others who were part of 'The Ensemble', Edward, Willem and B.Yang.

So many teams at the top all have a proven history in data mining competitions. I will be very happy if those ensemble dudes repeat their performance in the Netflix Prize ;-).

1 minute of fame

2011-08-19T03:34:00.000-07:00

The HHP competition and Kaggle are catching the attention of the media over here in Australia.

I was on TV recently (and yes, it does add a few Kg to you),

http://www.abc.net.au/catalyst/stories/3296837.htm

and in the paper (2nd last paragraph interesting!),

http://www.smh.com.au/business/miners-of-a-different-kind-in-line-for-riches-of-their-own-20110731-1i6bi.html

Blog posts may be limited in the next 2 weeks as its nose down to try and catch Edward.

Bath Mawr Coch

2011-08-03T09:02:00.000-07:00

In one of my first posts you will see how my eagerness to gloat of my success unleashed Dave into the HHP. He has subsequently been at the top of the leaderboard from the beginning, whereas my progress was more gradual. Eventually we seemed to converge to a similar score and even held positions 1 & 2 for a while, which was very surprising to me given the caliber of contestants that this competition has attracted.

We have now decided to join forces, hence the appearance of team Bath Mawr Coch. Hopefully there will be some synergy in our methods.

Dave and I are both competition junkies. We first cam into contact as a result of the KDDCup 2004 contest, which Dave won - and you will see my effort in a lowly 11th place. The leaderboard was left open after the conclusion of the contest, and I eventually got to top spot, which is obviously not that impressive when the winners had already revealed what they did - but I had learnt things and been able to implement them. I then contacted Dave to gloat - and he immediately overtook me.

We have since worked on a couple of other competitions together, and won the NCDM Analytic challenge in 2008.

Competitions like the HHP are a great way to learn and to also meet others with similar interests. I'm sure many collaborations and friendships will be made as a result, to the benefit of everyone.

Gone Shopping

2011-08-02T04:50:00.000-07:00

Another New Comp

Kaggle have just posted a new competition Dunnhumby's Shopper Challenge that is a very interesting one. The data is historical records of when customers visited a store and how much they spent. The goal is to predict when customers will next visit the store and how much they will spend.

This is probably a novel real world data set and I expect the interest in this competition to be high. The reasons I think this are

The data is simple
There is not a massive amount of data, so processing power will not be an issue
It is a novel problem where creativity is needed and a new algorithm will probably have to be developed. It is not going to be a case of who can pre-process the data best and build the best ensemble using existing algorithms

So, where to start?

The first thing I always do is run the data through my Nifty Tool (see here for example) to check it and generate some SQL so I can load it into a database.

Here is what comes out...

CREATE DATABASE dunnhumby
USE dunnhumby

CREATE TABLE training
(
    customer_id    int
,    visit_date    date
,    visit_spend    float
)


BULK INSERT training
FROM 'E:\comps\dunnhumby\training.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

--(12146637 row(s) affected)


CREATE TABLE test
(
    customer_id    int
,    visit_date    date
,    visit_spend    float
)


BULK INSERT test
FROM 'E:\comps\dunnhumby\test.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

The data set is relatively nice, just 3 coulumns;

customer ID
date
amount spent

the dates are from 1st April 2010 to 18th July 2011.

A Naive Guess

I always like to get in a submission without actually really looking at the data, just to make sure I am reading the data correctly and that the submission file is in the correct format.

In forecasting, a naive prediction is not going to be far wrong in the long run. If you want to forecast the weather, then saying tomorrow is going to be the same as today is going to serve you well in the long run if you have no other information to go off. So for this problem, we will say the next spend will be the same as the last spend, and the gap in days will also be the same.

To achieve this in SQL is a bit tricky, but can be done. The code below does the trick and when submitted gets 9.5% correct on the leaderboard, the same as the simple baseline benchmark.

/*************************************************
SQL to generate a naive prediction
- spend will be the same as the last visit
- the next visit will be in the same number
  of days as the gap between the previous visit
**************************************************/
 
-- append a visit number to the data - visit 1 = most recent visit
select *, Rank() over (Partition BY customer_id order by visit_date desc) as visit_number
into #temp1
from test

-- create field which is days since previous visit
select a.*
,b.visit_date as previous_visit_date
,DATEDIFF(DD,b.visit_date,a.visit_date) as days_since_previous_visit
into #temp2
from #temp1 a inner join #temp1 b
on a.customer_id = b.customer_id
and a.visit_number = b.visit_number - 1
where a.visit_number = 1

-- generate the submission file, makink sure 1st April is earliest data
select

customer_id
,
    (case 
        when dateadd(dd,days_since_previous_visit,visit_date) < '2011-04-01' 
        then '2011-04-01' 
        else dateadd(dd,days_since_previous_visit,visit_date)
    end) 
    as visit_date
,visit_spend
from #temp2
order by customer_id

Now Lets Look at the Data

In order to get a feel for things, we will aggregate the data to a daily level and look at how the population as a whole goes shopping. We will create a count of the total number of visits per day and generate some time based fields that will help us understand the daily and seasonal patterns in shopping. The SQL below will generate that data, which is small enough to save in Excel or as a text file and then load in your favourite analysis package, or we could just create a table in the database and load directly from there (don't leave as a view as it takes a while to generate).

-- add a field which is the days since 1st April (1st April 2010)
alter table training add daysSinceStart int
update training set daysSinceStart = DATEDIFF(dd,'2010-04-01',visit_date)

-- create a daily summary dataset
select visit_date

-- fields of interest
, COUNT(*) as visits
, avg(visit_spend) as avg_visit_spend
, stdev(visit_spend) as stdv_visit_spend

-- time based predictor variables
, min(daysSinceStart) as daysSinceStart
, sin(2 * pi() * (DATEPART(DAYOFYEAR, visit_date) * 1.0 / 365.0)) as doySin 
, cos(2 * pi() * (DATEPART(DAYOFYEAR, visit_date) * 1.0 / 365.0)) as doyCos
, (case when (DATEPART(WEEKDAY, visit_date)) = 1 then 1 else 0 end) as dowSun
, (case when (DATEPART(WEEKDAY, visit_date)) = 2 then 1 else 0 end) as dowMon
, (case when (DATEPART(WEEKDAY, visit_date)) = 3 then 1 else 0 end) as dowTue
, (case when (DATEPART(WEEKDAY, visit_date)) = 4 then 1 else 0 end) as dowWed
, (case when (DATEPART(WEEKDAY, visit_date)) = 5 then 1 else 0 end) as dowThu
, (case when (DATEPART(WEEKDAY, visit_date)) = 6 then 1 else 0 end) as dowFri
, (case when (DATEPART(WEEKDAY, visit_date)) = 7 then 1 else 0 end) as dowSat

from dbo.training
group by visit_date
order by visit_date asc

Tiberius is an ideal tool for analysing this type of data, as I initially wrote it to look at a similar type of data set - electricity consumption - and the patterns I want to discover will be similar.

I have built a model to predict visits as a function of time. The model predictors are weekly, seasonal and trend components, so we are trying to model the number of visits based solely on time. The plot below shows the errors in the model, the points are in date order.

(Click on the plots to enlarge)

What we see are that there are certain dates where the errors stand out. These are easily identifed in Tiberius by rolling your mouse over the image and the dates are displayed.

The two big anomalies are 25th Dec 2010 and 1st Jan 2011 - Christmas Day and New Years Day. So if we are trying to deduce where this data set is from, China is probably crossed off the list.

Other dates with errors are

4th April 2010
3rd May 2010
31st May 2010
30th August 2010

A bit of research will identify that these are all England & Wales public holidays, so I think we have figured out where the data set originates from (hello Tesco!).

The model of visits can be decomposed into weekly, annual and a trend components.

The weekly component shows Sunday is the quietest day, and Friday the busiest day. Monday, Tuesday and Wednesday are all similar.

The Seasonal component shows early April is the busiest time of year.

And unfortunately the trend looks to be going downhill.

Finally here is an interesting plot of the standard deviation of the spend,

So we are on our way to understanding this problem.

Good luck everyone, hope the challenge keeps you busy!

No Change at the Top

2011-07-30T17:24:00.000-07:00

I just spotted something on the HHP leaderboard that I have not seen before. There has been no change in the top 12 positions for 1 week.

Does this mean we are all running out of ideas? Teams are submitting but the incremental improvement seems to be relatively minuscule.

Dave (He is smarter than me!) and Willem seem to have a big gap at the top, and Edward is close on my tail for 3rd place. The team I would look out for though are some Old Dogs With New Tricks.

A1 Melbourne and MelbURN: R vs SAS

2011-07-29T01:52:00.000-07:00

Last night was another double header in the Melbourne analytic scene. Someone else has already described the events I attended, so pointing to their post will save me some typing.

A1 Melbourne and MelbURN: R vs SAS

Damn Computers II

2011-07-20T17:01:00.000-07:00

In an earlier post I grumbled a bit about how computers can be very frustrating...

http://anotherdataminingblog.blogspot.com/2011/05/damn-computers.html

here is another one for the collection!

Chrome v Internet Explorer 9

As described in the previous post I have just upgraded my PC. Now generally my experience with Windows is that when new, computers are quite quick but with time they seem to slow down and will eventually take 30 mins just to turn on. There is a school of thought that this is an inbuilt Microsoft conspiracy to ensure you keep upgrading your operating system when they release a new one - which seems to be every second week at the moment. The other theory is that in time you install more software that want to do things for you without asking (why does Adobe think I need it to check for updates every 5 mins).

Wanting to not clutter my new machine I am resisting installing anything that is not essential at the moment. I declined the 'trial' version of the latest Microsoft Office suite during the install of the o/s. This caused a problem when I quickly wanted to jot some numbers down in a spreadsheet. Rather than install my old licenced version of Office 2000, I downloaded the free Open Office and was very impressed - seamless integration with Microsoft Office so no issues you won't be able to read files others send you and it basically looks the same - I think I will be sticking with that from now on.

I also resisted downloading the Firefox browser and stuck with IE, until now that is.

I wanted to reply to a blog post. Having spent a while composing the reply, I tried to submit it and got stuck in an endless loop of being prompted to log into my Google account. Having experienced things like this before, my first option these days is to try a different browser. I went for Chrome this time, just on the hunch that seeing Chrome was a Google offering, as is Blogger then it should work. And it did - first time.

Having developed a few web sites, these differences between browsers are a bit of a pain. You really have to test everything out for all scenarios these days. On the Kaggle leaderboard I'm getting different views depending on what browser I am using. In IE my team name is just in bold, in Firefox and Chrome the whole cell is red making it easier to spot me, which I think is the look that it intended.

Internet Explorer - no colouring of the cells

Chrome & Firefox - I think this is the intended look

Keyboard Issues

With this new PC, I am having trouble typing a single quote ( ' ), which is required quite a lot in programming. If I press the key once, nothing happens. I then press it again an two of them appear and I have to delete one to get the single quote required.

I'm sure there is a setting somewhere to prevent this as I've never had this issue before. Rather than me waste all day trying to figure it out, I'm hoping someone can point me in the right direction.

Thanks in advance....

Pump Up The Volume

2011-07-20T05:24:00.000-07:00

The HHP competition convinced me that it was time to upgrade my tired old computer. As data sets get bigger, having more RAM and a 64 bit operating system can make things a little quicker.

For those who aren't aware, if you have a 32-bit operating system then there is no point having more than 4GB of RAM, as it can't be used. So if you need more RAM, then you will also need a 64-bit operating system.

I had a chat with the guy in the computer shop and he convinced me that it was going to be worth my while just to buy a new PC rather than upgrade my old one. My specifications were 16GB RAM, Windows 7 64-bit and the latest i7 processors. Size was also a consideration - due to the impending take over of my office by my daughter and me being relegated to the cupboard under the stairs.

I considered a laptop, but apparently the i7 processors that go in them are not as fast as the desktop versions. Instead I opted for a really small case, which increased the price as the case was more expensive and apparently you need a smaller motherboard to fit in them. A bog standard graphics card was added and Bobs your uncle.

The next day the shop rang me - with some bad news. The power supply that is required for the i7 processors won't fit into the nice little case I had opted for, so they recommended another for me that is kind of normal size (but still half the size and weight of my incumbent that comes with 5 fans and neon lights!).

The downgrade in case saved me some money, and the total cost was AU$1,300 - no monitor, keyboard or mouse included though. Incidentally, it would have still cost about $800 to upgrade my old PC (or so the salesman told me), and I now have 2 (actually about 8 as I never seem to get around to getting rid of my old ones).

Keen to see what it could do, I ran my old algorithms through the new beast, and it seemed to accomplish the same tasks as my Win 7 64-bit i5 laptop in 66% of the time. I then looked at the task manager to see what was going on under the hood.

The first thing I noticed is that with no applications running at all, just the operating system, it seems to be using 1.71GB, which seems a little excessive.

I then kicked off my algorithms and was a bit dismayed to see that only 14% of the available CPU was being used - so essentially only 1 of the 8 cores.

This is because the code I was running was doing things serially. It was doing a set of tasks that needed the prior task to be completed before the next one could begin. The extra processors would enable me to work on different tasks at the same time (such as browsing the Internet) - but the whole point of the investment was to speed up the number crunching.

One solution to utilise the extra processing power would have been to run several algorithms at the same time - but I wanted one algorithm to complete 8 times quicker.

A particularly popular algorithm is Random Forests. Now this is essentially an ensemble of lots of individual decision trees that are all built independently. There is no need to wait for one tree being finished before we can start building the next - an ideal candidate for parallel processing.

It took me a while to hack together a solution, but eventually I got 100% CPU usage and my model built in 1/8th of the time.

Below is the R code that will multi thread the building of Random Forests. The multithreading piece came from a generous Kaggle post in the Don't Overfit competition by Chris Pardy, and with that bit of info, it was trivial to adapt for my purposes.

Incidentally, Chris posted this code completely unexpectedly in response to some other code I had already shared. This was a nice surprise and is why the Kaggle concept is such a great learning environment.

There is now a dilemma though to sharing - a slight matter of $3million.

From what I understand, in order to win one of the milestone cash prizes you essentially have to tell everyone what you did in a way that they can replicate your results. If you do this, you are essentially putting yourself at the bottom of the leaderboard, as everyone else has your solution plus their own.

The milestone prizes increase in value every 6 months, so would it be a good strategy for the leader to not take the $30,000 in the first one and hold out for $50,000 6 months later or $60,000 12 months later, given that you will potentially be improving your model by learning from those who do take the money.

Any thoughts?

R code for parallel random forests

######################################
# some example code to multithread
# the building of random forests
######################################


#load the required libraries
library(randomForest)
library(snowfall)    
library(rlecuyer)

#load the data
theData <- iris

#set the formula
theFormula <- as.formula("Species ~ . ")

#specify number of threads
threads <- 8

#total number of trees in the random forest
totalTrees <- 1000

#number of trees to build per thread
trees <- as.integer(totalTrees / threads)


###########################
# now the funky stuff
###########################

#the function each thread calls
    parallelRF <- function(i){
    return(randomForest(theFormula, data=theData, ntree=trees))
}    
 
 # Initialise "cluster"
 sfInit(parallel = TRUE, cpus = threads, type = "SOCK")

 # Make data available to other R instances / nodes
 sfExport(list = c("theData", "theFormula","trees"))

 # To load a library on each R instance / node
 sfClusterEval(library(randomForest))
 
 # Use a parallel RNG to avoid correlated random numbers
 # Requires library(rlecuyer) installed on all nodes
 sfClusterSetupRNG()

# build the random forests
 allForests <- sfClusterApplyLB(1:threads, parallelRF)

 sfStop()
 
#everything finished, so
#merge all the forests into one
 myRandomForest <- allForests[[1]]
 for(i in 2:threads){
 myRandomForest <- combine(myRandomForest,allForests[[i]])
 }

#convince oursleves we have them all combined
myRandomForest$ntree

#what is important
myRandomForest$importance

Small World in Melbourne

2011-07-14T14:28:00.000-07:00

Having two young kids and no family in Melbourne (Australia) to babysit, my evening social activities are very minimal at the moment. They normally revolve around analytic type user groups that start after work and last for a couple of hours - then home for bath time.

Last night was a double header.

The Melbourne Tableau User Group had a meeting to demonstrate the new features of Tableau that are due to be released later this month. Tableau is a package that will draw you graphs and plot geographic data on maps. It is one of my favourite pieces of software for quickly drilling down on issues - you get to what you want very quickly. It costs ~$1,000 for a single licence or $2,000 for an advanced licence that lets you connect to ODBC data sources.

They gave away an IPad at the presentation - as one of the new features is dashboards for mobile devices. With only 27 people there I was feeling lucky, but #5 won (who incidentally already had one and was using it in the meeting to take notes), I was #9.

I spoke to the new Melbourne based Tableau employee, and asked him what he did, as I thought the biggest sales people for Tableau are the users themselves. He agreed - it is a great product and I've not heard a bad word about it.

The meeting was concluded with free beer and finger food that was nice. Here I met one of only two previous Australian KDDCup winner Bhavani Raskutti, who I know very well (the other was the University of Melbourne).

Immediately as this event drew to a close at 6:30, it was 100m down the road to a Data Melbourne talk. This one was even more social as a beer was thrust in my hand before the talk even began.

The reason I attended this event (permission for a leave pass had to be booked in advance) was that the speaker was non other than Matt Kwan from the HHP team Planet Melbourne who are currently 6th on the leaderborad. Matt talked about his PhD reasearch on using mobile phone data to analyse daytime population movements, which was really quite interesting. He then talked about the HHP, just basically introducing the concept of the competition.

He didn't give a way any secrets, but some interesting outtakes were that he had never heard of Kaggle before and was drawn to this by the lure of big money. He was an accomplished programmer but had done no analytics and had to learn linear algebra from scratch and is coding all his own software and wishes he had a pc with more grunt. All in all he is doing very well.

One interesting comment was that wouldn't it be more profitable to put all the effort into creating your own startup than spending all these hours for a potential but not very likely $500k. It is worth noting that several companies started up as a result of the Netflix Prize, including Commendo.

Melbourne is a small place. It is also the spiritual home of Kaggle. In fact I saw Anthony Goldblooms first ever presentation of the Kaggle concept at AusDM 2009, well before it was actually launched (incidentally I was running a data mining comp for this conference that was about blending solutions from the Netflix prize to get the best mega ensemble - may be handy for the HHP?).

At an R user group about 12 months ago, Anthony and I were sat at the same table as one Jeremy Howard, who became intrigued with the Kaggle concept and went on to become Kaggles most successful competitor and now part of the Kaggle team.

So, small world in Melbourne - my adopted home and a great place to live and work.

The Power of Twitter and Analyst First

2011-07-08T15:25:00.000-07:00

I don't use twitter, but it must be important as I have just bought a new Panasonic DVD player that can connect to the internet - but the only things you can seem to do are look at Twitter, YouTube, Picasso and the Weather. I also watched a presentation from a colleague recently and every 10 seconds these little messages kept popping up on his screen saying who had tweeted what. Not being Twitter user I had not seen this before, but it makes you realise that information exchange is instantaneous these days.

I am saying this as I had another big spike in views of this blog, in a specific half hour. Being a data miner and generally curious chap I obviously wanted to know what was the cause of this, or more precisely, if I had the skills to find out what the cause was.

The plot below shows the spike...

and it occurred in a specific half hour exactly 19 hours ago...

and a look at the web stats for the day pointed me to the source...

and then a search of Twitter pointed me to the guilty party...

Note this was tweeted exactly 19 hours ago, the exact time the spike occurred.

Now this is all very interesting you may say, but so what?

My point is that you need a curious mind to be a good data scientist - not expensive software.

In business today, those who make the decisions think that buying the most expensive analytic software is going to solve all the problems know to man, and you can pluck any monkey out to use it. To me this is putting the cart before the horse. If you want the best results, this could not be further from the truth. You need to spend money on the best analysts first and then think about the software.

I am currently 2nd in the HHP leaderboard, and have not used any software that was not freely available or written by myself to generate my predictions (except for some commercial compilers bought years ago and maybe Windows if you count that).

There is a new movement called Analyast First that has started up here in Australia. If you are an analytic manager I recommend watching the video below.

Analyst First 101

If you build it, they will come

2011-06-30T14:14:00.000-07:00

Or will they?

(Apparently, the actual quote is "If you build it, he will come" from Field of Dreams.)

I figured if I am to make money out of this blog by writing about the HHP, then it would be a good advert if I was actually doing well in the competiton. I recently have been putting effort in on this front, and team Sali Mali eventually got to the top of the leaderboard.

Being a data geek, I looked at the subsequent stats of this blog and saw a huge spike for a particular half hour around the time I got to the top. 'Bingo' I thought - thats the way to generate traffic.

I then looked at my linkedin stats, as the only real way to get to this blog is via my Kaggle profile, which will take you to my linkedin page and then to this blog. Surprisingly there was no such spike there - so what was going on?

Blog Views:

LinkedIn views:

As I don't know what exact time zones everything is in, my blog stats pointed me to this post on the kaggle blog, which I am assuming is the cause for the spike.

http://www.heritagehealthprize.com/c/hhp/forums/t/664/cross-validation-discrepancies/4381#post4381

Anyway, the data mining point of all this is that sometimes people are quick to jump to conclusions that are completely wrong. The real answers are always in the data - which is why I think the HHP will be won by a data scientist and prior expert medical knowledge will pay no part at all.

The more you practice, the luckier you become

2011-06-18T14:27:00.000-07:00

Chris R has now opened up a gap at the top of the leaderboard. It may look like nothing, but when you plot it, it is quite significant.

If you plot the number of entries per team, we see that progress seems to come with more effort. Chris R and a few other teams are almost doing daily submissions.

I am now progressing steadily as I investigate more of the data and more algorithms. I have gone from 40th to 19th by just refining the algorithm, rather than adding any different data.

The technique I am using to determine if things are improving are quite simple. Rather than have an independent set to test the errors against, I randomly split the set in half, build a model on one set and then calculate the error on the rest. The error on the building set is the 'in bag error' and the error on the unseen set is 'out of bag' error. This random splitting is repeated many times with the predictions just averaged. What I am looking for is the OOB error to go down with improvements in the algorithm. This seems to work, and the leaderboard scores improve as by OOB scores improve, although the leaderboard error is about 0.01 higher than my OOB error.

For a sanity check, you can plot the distributions of the predictions by in bag, out of bag and leaderboard set. The plot below shows a remarkable similarity, meaning the data distributions in the sets are very similar.

Data Anomalies

2011-06-08T14:38:00.001-07:00

This is a place where I will put some quirks of the HHP data that may or may not be relevant in wining the $3million

1) Be aware of the claims truncated flag for the claims in Y1.

2) There are some who have claimed for pediatrics who are clearly not children. This might indicate family accounts or incorrect age?

Scraping Up the Leaderboard

2011-06-05T18:50:00.000-07:00

Leaderboards in analytic competitions can be very revealing.

From a competitor viewpoint, you can see if you are in the right ballpark or completely missing something that others have obviously spotted.

If you are the competitor who is in the lead, then revealing your hand can be very detrimental as it encourages others to try harder and you tend to take your own foot off the gas if you are ahead.

I have been on both ends...

The right end a)

The pakdd09 Data Mining Competition was about predicting if credit card customers would pay their bills. Similar to the HHP, you had to predict into the future. You can see the leader board here.

Now the team at the top of the training data leaderboard, 'Tiberius Data Mining Solutions' is my team, but we weren't always there. You will notice the team in 2nd place, Suncorp, is way ahead of the field and they had obviously spotted something no one else had. Curiosity to figure out exactly what it was eventually lead us to discover new insights and overtake them.

Curiously in this competition we ended up 4th and the winners came from nowhere on the training data leaderboard to win on the out of time dataset.

The right end b)

The plot above shows the final standings for the Tourism2 Kaggle competiton. On the training leaderboard, Kaggles own Jeremy Howard (in his pre Kaggle days) was in front by a country mile after only a handful of submissions. This set a benchmark to aspire to, and pushed others to try to figure out what Jeremy had discovered. In the end I got to 2nd place on the leaderboard but was still way off Jeremy. Eventually though on the holdout sample, my submission did prove to be the winner.

The Wrong End

The Kaggle Informs contest was about predicting stock prices. This was one where I was surprised to make a submission that shot to the top of the leaderboard by a wide margin and stayed there for a while until the last few hours when I was overtaken by 3 teams. I had taken my eye off the ball and but would like to think had encouraged the eventual winners to squeeze every last bit of juice out of the data.

Leaderboard Plots

this blog post on the recent Chess comp is very interesting and studies how the leaderboard evolved through time and how big leaps in scores are made.

I ran a recent Kaggle competition called 'Don't Overfit' and posted some leaderbaord plots to the forum during the course of the comp. The way I did this seemed quick the first time around - just copy and paste into Excel and plot a graph. Doing this as a one off got me to see what I wanted fast, but it became tedious in time, doing it every week.

I then figured this task could be automated, and a bit of Googling soon provided me with the clues I needed to proceed. If I run the R code below, it will plot me the current leaderboard for any Kaggle comp and show me where I am (and how much more effort I need to put in).

What is interesting about the current scores is that maybe the top 5 teams are a bit ahead, but everyone else is on a nice trajectory, unlike some of the other leaderboards we have seen where there are marked step increases.

My own score is going slowly in the right direction as I explore more of the data and tighten up the techniques, but I am going backwards as others are doing this at a quicker pace!

The R code below will plot the live leaderboard positions as shown above...

##########################################
# Some R function to plot your position
# on a Kaggle leaderboard
##########################################

plotKaggleLeaderboard <- function(theURL,myTeam,topX=100,title){
    
#this library does all the hard work
#for windows users to install this packages see
#http://cran.r-project.org/bin/windows/contrib/r-release/ReadMe
#and http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.13/
library(XML)

#go and read the tables from the web page
tables <- readHTMLTable(theURL)

#get the table of interest
#names(tables)
lb <- tables[['leaderboard-table']]

#see what the columns are
#colnames(lb)
#nrow(lb)

#convert to numeric - see ?readHTMLTable
#numeric_cols <- c('AUC','Entries')
#numeric_cols <- c('RMSLE','Entries')
numeric_cols <- c(4,5)
lb[numeric_cols] = lapply(lb[numeric_cols], function(x) as.numeric(gsub(".* ", "", as.character(x))))


#extract the team name from a messy field
team_col <- c('Team Name')
#lb[team_col]

#split the field by "\r\n" than denotes the break between fields within a field 
lb[,team_col] <- sapply(strsplit(as.character(lb[,team_col]), "[\r\n]"), function (x) x[1]) 


#####################
#      plot
#####################
myRank <- which(lb[team_col] == myTeam)
mytext = paste("Team",myTeam,"is in position",myRank,"out of",nrow(lb))

error_data <- lb[,4]

plot(error_data[1:topX]
    ,col = 'blue'
    ,main = paste(title,"- leaderboard at",Sys.Date())
    , sub = mytext
    , xlab = 'Rank'
    , ylab = 'Error'
    )

#mark my position
points(myRank,error_data[myRank],col="red",pch=19)

}
###########################
# End of Function
# plotKaggleLeaderboard
###########################




compURL <- "http://www.heritagehealthprize.com/c/hhp/Leaderboard"
compTeam <- 'Sali Mali'
compTopX <- 125
compTitle <- 'HHP'



plotKaggleLeaderboard(theURL = compURL
                      ,myTeam = compTeam
                      ,topX = compTopX
                      ,title = compTitle 
                      )

There’s more than one way to skin a cat

2011-05-30T22:54:00.000-07:00

In a previous post we examined at the strange looking log function that we are supposed to minimise to win the $3 million HHP and showed that all errors are not equal with this metric.

For example, a difference of 0.1 between your prediction and the actual value is more severely punished if the actual is 0 than if the actual is 15. The point of this is because it is more important to know if the customer is going to go to hospital or not than knowing exactly how many days they will be in hospital.

Now this 'non-standard' way the model is being judged might cause issues if we are using standard algorithms that seek to minimise different error functions. For example, in linear regression, the algorithm seeks to find a solution that minimises the sum of the squared errors, with no regard for what the actual value is;

sqrt( ( pred - act ) ^ 2 )

So, does this mean that we have to come up with a new algorithm that seeks to minimise this specific log function, or is there anything else that can be done so that we don't have to reinvent the wheel?

Now the error function we are asked to minimise is very similar,

sqrt( ( log(pred+1) - log(act+1) ) ^ 2 )

Note that the fact we have taken the log(act+1) is what introduces the issue that the magnitude of the actual value plays a part. The plot below shows x v log(x+1). A small change in x will give different changes in log(x+1) depending on what x was in the first place (ie the gradient of the line).

So, what to do?

If we look at the two equations above we can see that they are the same if we substitute

pred1 = log(pred + 1) and act1 = log(act + 1)

So, if rather than predict DIH, we predict log(DIH + 1), then we can use least squares minimisation, and use some standard algorithms that have already been developed for us.

What we need to remember though is that when we make a submission we have put the prediction back on the right scale.

DIH_tran = log(DIH + 1)

DIH + 1 = exp(DIH_tran)

DIH = exp(DIH_tran) - 1

If we plot the data with these adjustments we can see how the error surface shows parallel lines, meaning all errors are equal, and least squares minimisation will work OK. When we put the errors back on the original scale, all errors are not equal.

#############################################
# NOW ADJUSTING THE TARGET SO WE CAN
# MINIMISE THE RMSE
#############################################


# generate the data
dih <- seq(from=0, to=15, by = 0.05)

#adjusted days in hospital
dih1 <- log(dih + 1)
dat <- expand.grid(act1 = dih1, pred1 = dih1)

#adjusted error function
dat$err1 <- sqrt((dat$pred1 - dat$act1) ^ 2)

#'un-adjust' the days in hospital
dat$act <- exp(dat$act1) - 1 
dat$pred <- exp(dat$pred1) - 1

#plot the error on the adjusted scale
contourplot(err1 ~ act1 * pred1, data = dat
          ,region = TRUE
          ,cuts = 10
          ,col.regions = terrain.colors
          )
    
#plot the real scale
contourplot(err1 ~ act * pred, data = dat
          ,region = TRUE
          ,cuts = 10
          ,col.regions = terrain.colors
          )

Hunting the Hidden Dimensions

2011-05-29T00:43:00.000-07:00

Last week I came across a documentary on TV called 'Hunting the Hidden Dimensions'. Now generally there is nothing worth watching on Australian TV (unless you are into cooking or talent contests) but this one actually got me interested (although I must admit to falling asleep half way through due to the late for me time slot).

http://www.sbs.com.au/documentary/program/huntingthehiddendimensions/index

The program was about fractals and the hidden order that appears in everything in nature. I didn't realise that fractal geometry was involved in film animation to generate mountain ranges and volcanic lava flow!

One thing that did strike an accord was how Mandelbrot, when he first published his work, got a bit of a thumbs down from the academic mainstream because the ideas didn't conform to the nice orderly nature of the world people were used to (ie the type of mathematics the Egyptians used to build the pyramids). The thing is, fractals are perfectly orderly, but looked at in a way that was different to the norm.

This is in a data mining context is similar to neural networks, with many people who should know better saying they are uninterpretable black boxes. When neural networks give improved predictions, they have found 'hidden dimensions' that are very easy to find if you are willing to look. Like Mandelbrot, you will need to look at data visually rather than mathematically in order to figure things out.

You can watch the documentary online.

Minimise the what?

2011-05-21T22:08:00.000-07:00

The other week I was on Amazon to purchase Data Mining with R, and Amazon did some cross-selling, telling me I might be interested in ggplot2. Now I remember listening to a talk by Jeremy Howard, and seemed to recall he mentioned this book was worth a read, so I ticked the approproate box, and now the books are here.

I prefer looking at pictures to interpret information, so I wanted to see if I could use my new book to figure out what this log function was all about that the HHP prize people want us to minimize.

The following plots are all variations of the same data, plotted with different packages and functions. The generated data is basically what the errors would be for all possible combinations of Actual Days In Hospital v Predicted Days in Hospital.

Now generally, if you predict 15 and the real answer is 14, that is considered just as good as if you predict 2 and the answer is 1 - both have the same error - 1. Many methods of determining what is a good prediction, such as the mean squared error (MSE), don't care about the magnitude of the prediction, just the magnitude of the error.

Other metrics such as mean absolute percentage error (MAPE) do consider the magnitude. In electric load forecasting, if the actual load is 2MegaWatts (MW) and you have forecast 1MW, it can be far more serious than if the actual load is 900MW and you have forecast 899MW. In both cases the error is only 1MW but the consequences can be quite different.

So what does the error function of the HHP tell us?

Basically it says that we are penalised less for the same absolute error the higher the actual Days in Hospital is. So rather than concentrate on improving the prediction for someone who was 15 days in hospital and we only predicted 14 days, we would get far more reward by putting the effort into improving the prediction for someone who was 2 days in hosptal and we only predicted 1 day.

or...

a small improvement in accuracy for a small Days in Hospital is worth as much as a large improvement in accuracy for a large Days in Hospital.

Now all this insight was gleaned from just creating a graph. No Maths!

Now if we zoom in to what is the business end (remember 0.2 is a good guess) we see things are a bit more linear...

And here is what a we get for a normal type of error (just strip out the log and +1 from the data generation function). Note the lines are parallel.

and for percentage errors, note the lines are converging, similar to the log error, but the worse position to be in is making high predictions when the actual is low.

Here is the R code to generate the plots...

#############################################
# generate the data
#############################################

dih <- seq(from=0, to=15, by = 0.2)
dat <- expand.grid(act = dih, pred = dih)
dat$err <- sqrt((log(dat$pred+1) - log(dat$act+1)) ^ 2)
dat1 <- matrix(sqrt((log(dat$act+1) - log(dat$pred+1)) ^ 2),length(dih))



###########################################
# now plot the error surface of act v pred
# using several methods
###########################################
library(ggplot2)
library(lattice)
library(vcd)
require(grDevices) # for colours


### using ggplot2 ###
df <- data.frame(dat)
names(df) <- c("actual","predicted","error")

#version1
errplot <- ggplot(df, aes(actual,predicted, fill=error))
errplot <- errplot + geom_tile()
errplot <- errplot + scale_x_continuous(expand = c(0,0)) 
errplot <- errplot + scale_y_continuous(expand = c(0,0))
errplot

#version2
fill_gradn <- function(pal) {
  scale_fill_gradientn(colours = pal(7),limits = c(0,.5))
}
errplot <- errplot + fill_gradn(rainbow_hcl)
errplot


### using wireframe ###
wireframe(err ~ act * pred, data = dat
          ,scales = list(arrows = FALSE)
          ,drape = TRUE
          ,colorkey = TRUE,
          screen = list(z = 30, x = -60))
    
contourplot(err ~ act * pred, data = dat
          ,region = TRUE
          ,cuts = 10
          ,col.regions = terrain.colors
          )
    
    
## filled contour plot
filled.contour(x = dih, y = dih, z = dat1
                ,nlevels=10
                ,color = terrain.colors
                ,plot.title = title(main = "HHP Error Function (the funny log one!)",
                xlab ='actual DIH', ylab = "predicted DIH")
                ,key.title = title(main="Error"),
                )

State of the Union

2011-05-18T13:13:00.000-07:00

80%+ of a dataminers time is normally spent manipulating the data so they can get to a point where they can start extracting insights. This data preparation is actually the most important part of the process where you should be most diligent.

Here is another example (to add to the one in the previous post) that demonstrates you don't always get what you think you should get.

We have been supplied with two data sets of names and want to combine them into a single table. The UNION command will oblige us, but we have got to be very careful how we pose the question...

(This is MS SQLServer)

The result of this is not what we might think. Union seems to ignore the actual name of the field, so our firstnames and surnames are now all muddled up because they were in different column positions in each source table.

So to get what we really wanted, we can get the columns aligned first...

Don't assume the source data is consistent and the fields are in the same order in the database. Fortunately if you have lots of fields there is a way to check this...

If the field orders are not consistent, then it can be a pain explicitly typing out all the fields in your SQL statement to get the orders to match. Now there is also a time saver for this...

We can then just copy the result and paste it into our SQL to save us typing.

Beware the Integers of SQL

2011-05-18T03:52:00.000-07:00

Most data miners improve by learning from their mistakes. Many software programs have their own little quirks and sometimes deliver answers that are not quite the answers to the question you thought you were asking...

select AVG(DaysInHospital) from DaysInHospital_Y2
-- 0

select AVG(DaysInHospital * 1.0) from DaysInHospital_Y2
-- 0.467101

select SUM(DaysInHospital) / COUNT(*) from DaysInHospital_Y2
-- 0

select SUM(DaysInHospital) * 1.0 / COUNT(*) from DaysInHospital_Y2
-- 0.467101542670

Learning from the Leaderboard - Part 1

2011-05-16T07:40:00.000-07:00

It is now day 12 of the HHP. I know this as I have been trying to make a submission every day, and the leaderboard says Sali Mali has made 12 submissions.

Dave is keeping his position near the top but Sali Mali is going downhill, as others overtake me.

I've now started to look at the problem in a little more detail, and decided I needed to a way of replicating the error function that is being used to judge this competiton. This will help in figuring out what sort of leaderboard score I expect my submissions to achieve.

I've also not been totally wasting my submissions during the time when I have had no useful model to submit. In previous time series competitons, it has been possible to learn a lot from the leaderboard (real examples will be given part 2 of this post), so I figured pinging a few constant valued submissions might help me extract some insight about the future (although I've not really thought about what that insight might be at the moment).

The following R code demonstrates

how the function to calculate the error was put together
how we can easily generate the errors for a sequence of constant valued predictions on Y2 and Y3
how the Y4 leaderboard results were read from an Excel file and compared to the Y2 and Y3 results.

Doing this was a sanity check that my error function gave similar results to the leaderboard calculations. It shows that a constant of 0.2 is a good guess and also that the Y4 curve looks more like Y2 than Y3. This is interesting and needs further thought on exactly what this means, if anything!

###########################################
#function to calculate the model error
###########################################
calc_error <- function(act,pred)
{
    aact <- as.matrix(act)
    ppred <- as.matrix(pred)
    
    if(nrow(aact) == nrow(ppred)){
    return (sqrt(colSums((log(ppred+1) - log(aact+1)) ^ 2) / nrow(aact)))
    } else {
    return (-99)
    }

}
#### EOF to calcualte model error ####



###########################################
# function to calculate the model errors
# given a sequence of constant values for
# the predictions and the known outcome
###########################################
err_seq <- function(sequence,act)
{
    
    actual <- as.matrix(act)

    #vector of errors
    errors <- vector(length=length(sequence))

    #vector for the predictions (will be constants)
    predictions <- vector(length=nrow(actual))

        ind <- 0
        for(i in sequence){
            
            #report progress
            ind = ind + 1
            cat("\n",i)
            flush.console()
            
            #set the prediction to the constant     
            predictions[] <- i    
            #calculate the error
            errors[ind] <- calc_error(act=actual,pred=predictions)
        }

    return (errors)
}
###EOF to calculate sequence of errors ####




#connect to data
library(RODBC) #for data connection
conn <- odbcConnect("sql_server_HHP")

#load the actual days in hospital
actualY2 <- sqlQuery(conn,"select DaysInHospital from DaysInHospital_Y2")
actualY3 <- sqlQuery(conn,"select DaysInHospital from DaysInHospital_Y3")

#set up a sequence of constants to be used as the predictions
const_preds <- seq(from=0, to=1, by=0.01)

#calculate the error sequences for Y2 & Y3
Y2Err <- err_seq(sequence=const_preds,act=actualY2)
Y3Err <- err_seq(sequence=const_preds,act=actualY3)

#read in the errors for Y4 from an excel file 
#where they were recorded
library(xlsReadWrite)
#this is also required to be run
#see http://www.swissr.org/software/xlsreadwrite
#xls.getshlib() 
Y4 <- read.xls("E:\\comps\\hhp\\constants.xls")



#############################
#plot prediction v error
#############################
ymin <- min(Y2Err,Y3Err)
ymax <- max(Y2Err,Y3Err)
plot(const_preds,Y2Err
    ,type='l'
    ,main='Learning from the Leaderboard'
    ,xlab='Predicted Constant'
    ,ylab='Error'
    ,ylim=c(ymin, ymax)
    ,col='blue'
    )
lines(const_preds,Y3Err,col='green')

#add the line for Y4, discovered via leaderboard
lines(Y4[,1],Y4[,2],col='red',lwd=2)

legend('bottomright', c('Y2 - known','Y3 - known','Y4 - via leaderboard'),lty=1, col=c("blue","green","red"))

Connecting R to our Database and Plotting Histograms

2011-05-13T16:06:00.000-07:00

So far with the HHP I have

loaded the data into a database
written an R script to check the submission file is correct, with the added bonus of plotting a distribution for a sanity check

As yet I haven't even really looked at the data and don't really know what the problem is we are being asked to solve. But not to worry, all that will hopefully fall into place once we have our tools all set up and talking to each other.

Data lives in a database, but it is no good to anyone if we can't get at it to analyse. It is possible to access this data, as the database developers have very genourosly provided what are called database drivers that allow these other applications to get to the data, and even change it if so desired. Many analytic applications have provided the means to utilise these drivers so the data can be sucked into their tools for analysis.

The task of connecting your applications to databases can be very painful if you are not really a database person and don't understand the lingo of what you are being asked to do. Some applications make it very easy (I have found Tableau remarkably painless and my own software Tiberius has instructions on how to do this) and some you have to Google around a lot find the answers to your questions. It is amazing how many products in their documentation just say 'contact your system administrator' - which is no help if it is you.

In the end though, once you have connected that is the biggest hurdle over, and then the world is your oyster.

So, we now want to connect SQL Server to R. First the code, then the detail.

Note: If you are viewing this in Internet Explorer or Chrome, you might get text wrapping in the code below and no horizontal scroll bar. If you are using FireFox then you get the scroll bar and no text wrapping. Either way you can still copy and paste the code OK. Another Damn Computer!

# load the required libraries
library(RODBC) #for data connection
library(lattice) #for histograms

#######################################
# set up a connection to the database
#######################################

#method 1 - using a connection string
conn <- odbcDriverConnect("driver=SQL Server;database=HHP;server=PHIL-PC\\SQLEXPRESS;")

#method 2 - involves setting up a DSN (Data Source Name)
#conn <- odbcConnect("sql_server_HHP")


########################################################
# extract the data using this connection and some SQL
########################################################
mydata <- sqlQuery(conn,"select  * from DaysInHospital_Y2")
#mydata <- sqlQuery(conn,"select  * from Claims")


############################
# take a look at the data
############################
colnames(mydata)
head(mydata)
summary(mydata)


############################
# plot the distributions
############################

#set up the plot layout
maxplots <- 9 #upper limit on plots
plots <- NCOL(mydata)
plots <- min(plots,maxplots)
sideA <- ceiling(sqrt(plots))
sideB <- ceiling(plots/sideA)
cells <- sideA * sideB
#par(mfrow=c(sideA,sideB)) #use this for standard graphics
plotposition <- matrix(1:cells,sideA,sideB,byrow = T)

#draw all the plots
for(i in 1:plots){
myplot <- histogram(mydata[,i],main=colnames(mydata[i]),xlab='')
myposition <- rev(which(plotposition == i, arr.ind=TRUE))
print(myplot, split=c(myposition,sideA,sideB), more=TRUE)
}

The code above is R code and is run from R. It assumes you have the library RODBC installed. If you don't, then just type

install.packages('RODBC')

So, how do you get it to work?

There are 2 options to connect to the database, methods 1 & 2.

In method 1, you need to know 3 bits of information (you also might need user IDs and passwords depending on how access to the database is set up),

the driver name (SQL Server in our case)
the database server name. This can be tricky to remember. For SQL Server the easiest way to find it is when you start up the SQL Server Management Studio. In the picture below you will see the server name is PHIL-PC\SQLEXPRESS.

There is one gotcha here. You will notice in the R script we have an extra '\'

server=PHIL-PC\\SQLEXPRESS

This is because '\' is a special charater in R, so you need to put two of them together to tell it that you really do mean '\'. This is useful to know when you have paths of files etc.

Another thing to note above, we are in Windows Authentication mode, which means we won't have to worry with passwords.
The name of the database in the database server that you want to connect to. You can discover your options from within the Management Studio. Here we have HHP and test

So when we string all this together we get,

conn <- odbcDriverConnect("driver=SQL Server;database=HHP;server=PHIL-PC\\SQLEXPRESS;")

Method 2 involves creating a Data Source Name (DSN), which can be thought of as the above bits of information but stored away somewhere in your computer, just wrapped under a name. With many applications, you just then have to specify this name to connect to the data.

On Windows machines, a DSN is created via a wizard that takes you through the steps. Here are instructions for XP.

The one Gotcha that I commonly encounter with SQL Server is that it tries to populate the list of available servers so you can just select it from a list, sometimes you can be waiting forever and you are not guaranteed to get your server in the list.

You can see above that it lists PHIL-PC, but this will fail, it needs to be PHIL-PC\SQLEXPRESS.

Also, if you are connected to a network, the list of servers could be pretty big and picking the one you need would be difficult to remember. If you do this often it can be worthwhile just having your common servers in a list and just pasting them in.

With the wizard, it can test the connection for you, so you know if your settings are correct.

The end result of all this is that we can connect R to our database and plot a graph. We can now get a sense of what the data is all about. Progress!

Interesting observations - Days in Hospital is mainly 0, then 1-5, with another hardly visible blip at 15.

hist(mydata$DaysInHospital[mydata$DaysInHospital > 1])

Interesting, obviously an anomaly in the data? Or maybe we should now start to read about the problem we are trying to solve and how the data has been put together.

We can just change the query to look at a different table in the database. The Claims data looks like this...

This gives us an idea of what we are dealing with. Now to 'shake rattle and roll!'

Relevant Adds are Here!

2011-05-11T06:13:00.000-07:00

I posted a piece on loading in the data, and now I see the adds appearing related to this. This is good, I even wanted to click on them myself to see what they were, but I resisted as I want the kitty to be non funded by me.

So click away boys and girls if you like the adds. I am $2 down today already as I went to the laundrette to dry my clothes (rained in Melbourne today) and put my $2 in the machine to get some change and all that came out was a flattened hearing aid battery.

Beware the HHP submission sort order

2011-05-11T05:56:00.000-07:00

When a solution is submitted for the HHP, the file has to contain the memberID and the predicted Days In Hospital for year 4.

It was noticed early on by a competitor that the memberIDs in the example file provided were not in numerical order, in fact they were quite random. Kaggle confirmed that the submissions had to be in the same memberID order as the supplied example file.

This now causes a concern, being as the data is now in a database. There is no real concept of row order in a database, just because you enter data in a particular order, there is no guarantee it will come out in that order. If you need a specific order, then the order by statement should be used in the SQL to specify the order required.

So how can we get the specified order?

We need to add a new column to the Target.csv file that contains the row order. A simple way to do this in Excel is just to add an extras column going from 1 to the number of cases, in ascending order.

Another way to create this extra column is to use R. When R reads data it conserves the order the file is read. The code below will create a new file with an extra column.

#########################################################
# the customer ID is not in any particular order in
# the submission file. This code creates a sort order
# that we can use to ensure things are sorted correctly.
#########################################################
setwd('E:\\comps\\hhp\\raw data\\HHP_release2')

#read in the original
Target <- read.csv('Target.csv')

#add an extra column called orderbyme
#with the row number as the value
Target$orderbyme <- 1:nrow(Target)

#take a look at the result
head(Target)

#write out a new file
write.csv(Target, file='Target_with_order.csv',row.names=FALSE,na = "")

(the syntax highlighting is not perfect for R, it is for Perl, but lets you see the comments easily)

Now if we load this into SQL Server with the rest of the data, we have a field we can always sort by.

Once a submission file is generated, we can also double check it in R and plot the distribution to give it a sanity check.

#########################################################
# this R code reads in a submission file, checks it and
# plots the distribution of predicted days in hospital
#########################################################

#change this to the name of the submission file to check
file_name <- 'submission6.csv'

#the submission files and
#Target.csv should be here
setwd('E:\\comps\\hhp\\submissions')

#read in the files
sub_file <- read.csv(file_name, header=TRUE)
compare_file <- read.csv('Target.csv',header=TRUE)
 
#check the file has the correct number of rows and columns
if(nrow(sub_file) != 70942){cat('\n\n***invalid number of rows***\n\n')}
if(ncol(sub_file) != 3){cat('\n\n***invalid number of columns***\n\n')}

#check the example file
if(nrow(compare_file) != 70942){cat('\n\n***Target.csv is corrupted***\n\n')}
if(ncol(compare_file) != 3){cat('\n\n***Target.csv is corrupted***\n\n')}

#check the ids match
differences <- abs(sub_file[,1] - compare_file[,1])

if(max(differences) != 0){
cat('\n\n***IDs in wrong order***\n\n')
}else{
cat('\n\n***IDs are in the correct order***\n\n')
}

#check the range
if(max(sub_file[,3]) > 15){
cat('\n\n***you have predictions > 15 days!***\n\n')}

if(min(sub_file[,3]) < 0){
cat('\n\n***you have predictions < 0 days!***\n\n')}


#####################################
#plot the submission distribution
#####################################
num_bins <- 50
theplot <- hist(sub_file[,3]
    ,breaks=num_bins
    ,xlim = c(0, max(sub_file[,3]) + 0.1)
    ,main = paste("Distribution of ",file_name)
    ,xlab = "Predicted Days in Hospital"
    )

#the mean prediction
mean_pred <- colMeans(sub_file[3])

xpos = max(sub_file[,3]) * 0.5
ypos = max(theplot$counts) * 0.75
mytext = paste("Mean submitted DIH= ",sprintf("%.4f", mean_pred))
text(xpos,ypos,mytext,pos=4)

and this is what I get for one of my submissions to date...

and interestingly if you change num_bins to 5,000, you get...

which will lead to some investigation why we get this funny distribution.

Some are saying this guy will win the #drflix prize

2011-05-10T15:08:00.000-07:00

I suddenly started receiving emails saying I had a first visitor from Belguim, UK etc. What was all this about - I'd only told 2 people about this blog and was not ready to go public just yet.

A bit of data mining via the blog web stats pointed to twitter, and then I saw this on the HHP page next to the leaderboard. See the top post on the right hand side - I clicked on the link and was taken aback when it lead to this blog.

I'm not into twitter, so can someone explain what #drflix is?

The number of visitors is now 85 from 16 countries (about the same as the number of teams now entered) but the kitty is still empty, so I guess I need to start talking about data mining more to get some relevant adds. I'm also getting some error messages because of that web counter gadget, page loading slow because of video type things - it may have to go, when I get passed 100 visitors maybe.

As to who will win, Dave is looking pretty strong at the moment. He remains at the top - I knew I shouldn't have bragged to him about my early triumph. Meanwhile I am sitting at no. 23, tucked in the middle of the peleton. Hope to tell you how I got their soon.

Damn Computers

2011-05-10T05:45:00.000-07:00

MICROSOFT

Been having a few issues with my Windows 7 64 bit laptop. Basically the other day it decided to tell me not to switch it off as it was going to install 125 updates. After a few hours waiting everything finished, but some very useful software I had wrote using Microsofts very own compiler no longer worked. It did work once but it looked as if the font was Chinese, and in subsequent attempts it seemed to want to run as a service. This was odd as the very same executable worked fine on my 32 bit machine.

To try to solve this issue I recompiled on the 64 bit machine and everything seemed OK . To my delight it seems now that a single executable is all that is need to run on 32 and 64 bit machines, which is good. The new issue was that I compiled with a setting that says you need dot net 4, and my 32 bit machine didn't have that. Rather than make people have to install dot net 4, I tried to recompile with a setting that says you only need dot net 2, but when I changed this setting the compiler got stuck in some kind of neverending loop and didn't do anything.

The reason for saying this is that I am going to mention this software in future posts - but if you want to use it Microsoft decrees you will need to install the latest version of dot net, unless anyone can tell me how to solve this problem.

BLOGGER & FIREFOX

Another issue that I've had is with this blog. I tried to edit posts and again, the little circle wizzed around indicating that it was thinking, but and hour later it was still thinking. I also couldn't seem to get any love from the template settings options.

I was using Firefox, which had also recently updated itself, so thought I'd try to see if it was the browser. Being as the blog was a Google thing I thought I'd download Chrome. Now when I tried to goto the Chrome download page, I got the message from Google saying wait for the download and the same wizzing circles, which also wanted to be there for hours. I suspect this is a Firefox anti Google thing. I then tried IE and everything works fine in that, so Microsoft have redeemed themselves a little.

RATTLE

The next issue I had was with trying to get Rattle to run in 64 bit mode. There are a few documented issues of things you need to check, but I tried them all but couldn't get it to work. I have it working in 32 bit R, but 64 bit would be nice. Rattle is a nice useful bit of data mining software that I will hopefully be mentioning at some point in this blog.

PGP Encription

I have recently had some Pretty Good Privacy put on my hard drive, which at first glance seems a pretty good idea if your computer contains valuable information.

Now for me, the valuable information is large datasets that can contain sensitive details you don't want to find make their way into the newspaper should you leave your laptop on the train. The software I mentioned earlier reads these data files, but what I am finding is that it will read the first 30,000 or so rows in a few seconds, but then grind to a really slow pace after that, which is bad when you have a few million rows.

In order to eliminate the fact that it might be a result of the new compilation of the software, I moved on of the datasets to a different drive that was not encrypted, and voila - it worked as expected.

Java & Eclipse

I was using SPSS (a tool that can do things with data) and was quite merrily putting together a model when the window quite rudely decided it was not playing.

Initially I suspected it was an issue with SPSS rather than something I had corrupted on my PC , especially when a colleague who had exactly the same PC and OS got the same problem. This was an issue for me, particualry as the blackout was unannounced and made you lose any unsaved work.

SPSS was redeemed when I also then discovered it happedned in another application, unfortunately for me the one that helped me get my expenses reimbursed. The common denominator was that the applications were both built on Eclipse. All evidence pointed to the video drivers and java not being a happy couple.

I found there were some new drivers available for download and installed them. Hopefully this has fixed the issue, but as this issue happens unannanounced half way through some work, I have learnt to save things quite frequently, just in case it happens again.

Progress - Loading the HHP Data

2011-05-07T14:09:00.000-07:00

A little progress...

My daughter made it through her first swimming lesson without needing a swimming nappy. This was mainly my fault as we had ran out and the lesson was too early to have time to go to the shop. Anyway, that should save me some cash, so maybe I will no longer need the full 3 million.
I managed to find a way of posting code with syntax highlighting into this blog
I managed to load the HHP data into a database
I picked up some bugs and improved some of my software while doing 3).

So, now the details....

1) Well done Princess no. 3
(but please tell us where you hid the TV remote)

2) Highlighting Code in a blog

When you do data mining, you would generally at some point have to write some computer code that does stuff to your data. You will notice that when you write the code it generally turns nice colours to highlight things that make the code easier to interpret.

I will hopefully be doing lots of this and posting code snippets to this blog so others can get ideas or use the exact same code to try out (not sure how this fits in with the rules of the contest, but lets go for it anyway). It would be nice if the code also coloured itself in the blog.

A Google search came up with this site http://tohtml.com/. You just paste your code in a window, select what computer language the code is written in, press a button and presto, it generates html code for you. You then just edit the html of your blog and paste it in, and should get something like this...

Sub IsThisRight()

    MsgBox (IsNumeric("50+"))

End Sub

and if you tweak the background colour by editing the html...

Sub IsThisRight()

    MsgBox (IsNumeric("50+"))

End Sub

and if the code is big you can put it in a window by wrapping a div tag around in.

Sub IsThisRight()

    MsgBox (IsNumeric("50+"))

End Sub

Now this is probably not the most elegant solution and the background colour of the code is lost outside the code window - but I'm sure that can be easily fixed - but it works for me at the moment, but open to other suggestions.

One Gotcha is that when you edit the html in Blogger it is too smart and wants to add extra tags and rearrange things for you. When I pasted the div code above, height and width were next to each other, but Blogger decided overflow and auto should go in between them.

A call for help - how do I set the default font in Blogger? I can't find the option anywhere.

3) Loading the HHP data into a database

Well, at last something related to data mining...

When you are doing analytics as a consultant, you get given data files in typically 3 ways.

an Excel spreadsheet
some other proprietary data file format such as SAS or SPSS
a text file that is an extract from a database

the Excel spreadsheet can be because of 2 reasons:

Excel is the database used and where the raw data lives (worrying but very common)
The person doing the extract extracted it from a database for you and then decided they were doing you a favour by putting in Excel to pretty it up (also worrying)

I've had the following explanations given to me...

"There were a few million rows in the database but I can only get 64,000 on each sheet of excel so it will take me a while to get it ready for you."

"I extracted it but had to put it in excel because I new some of the data was wrong so I manually edited it for you."

When you hear this, just take a deep breath and do a bit of hand holding.

You often get proprietary files sent, which basically means it can be a pain to open them if you don't have the proprietary software. These are generally binary files and are the ones that look like hieroglyphics when you try to open them in a text editor. Binary files are a more efficient format for storing and retrieving data.

The bog standard ascii text file extract is the most painless way to go, which is the way it works in Kaggle competitions. There are two ways you could get these files:

fixed width
character separated variables (*.csv files)

The data represents fields/variables/columns/features and rows/cases/examples. Text files are clever and know when there is a new row of data, but specifying which bit of data should be in which column is another matter.

In fixed width files, each column starts at the same place in each row. This is not necessarily best as you may have to do a lot of padding with blank spaces to get the columns aligned, which makes the file bigger. If all the data in each column has the same width, then you will get a smaller file as it saves having to put a delimiter character in there. They are also not particularly simple to load - these are the ones in Access where you have to click to tell Access where the columns need to go. The field names are also generally not on the first row, they are tucked away in some other field definition file.

CSV (C Separated Values) files put a particular character in each line to indicate where the column boundaries are. There are a few rules on how the files should be generated (note: there used to be a wikipedia entry that said csv stood for character, not comma, but I can't find this anymore - but there is a new entry DSV.)

The common character is a comma (","), but using this is fraught with danger which can lead to your file not loading properly.

If you have a text field in your data where a comma is actually contained in the text, such as address ("1 Smith St, Melbourne") the text field need to be "qu,oted" in order to determine that the comma is part of the data rather than a delimiter. If these quotes are not included your data will get stuffed up and will have to be extracted again. This happens very often.

Also, some countries use the "," as a decimal point, which can also complicate things (see here).

When receiving data, I always ask for pipe delimited (this is the pipe character | ) files with no quotes around string fields. The pipe character is unlikely to be present in any text fields (as not many people actually know what it is) so is a fair choice for a delimiter. Not including quotes around text, although this is defined in the standard, is going to be more friendly to most software you will use. Tab delimited files are also another option.

Well, enough of the side track, back to the HHP data.

When you receive a delimited text file, I recommend you check it before you do anything - don't assume it has been extracted correctly. The check is simply counting that each row has the same number of delimiters.

Now you could load it into Excel to eyeball and it might 'look right', but Excel won't deal with bad files, and your columns may become miss-aligned. I also don't trust Excel as it also seems to have a mind of its own sometimes and does things to data that you rather it wouldn't, and might not even realise it has done (try typing (22) into Excel, including the brackets. Dates are also a concern for me).

You can also load the data into numerous databases using the import wizards they provide. Now some will tell you there are issues, but not all.

Microsoft SQL Server is one example of a database that used to be good at checking things. I often got messages such as 'error at row 2,758,946' and the would go and check the file at this row and discover there was an issue. You would then have to correct that row, re-load, get the next error message, correct, re-load etc.

Now, in order to show you what the error message was, I tried running an example in Microsoft SQL Server 2008 R2. I used the data import wizard to load the following CSV file:

and this is what it did:

This is not what I was expecting. You will see column c, row 2 has the value 3,4. This is probably not what was intended. Now the last time I used the import wizard, in version 2000, I'm sure it would have thrown an error, so beware version changes.

Anyway, data lives in a database, and there are many free ones available. Most of the major vendors offer free versions that have limitations on the size of the database, and there are also open source ones, see this list.

I am going to use MS SQL Server 2008 R2 which is free. It has a 10 gig limit per database, but you can create multiple databases and communicate between them, so unless you have a really massive file it should accommodate most needs. Installation from what I remember had a few tricky questions you probably will think 'what the **** is that on about' if you are not used to database lingo, but once installed it is a pretty user friendly tool.

In order to load the data, you can use the data import wizard, but this is not something I personally like. I like to run scripts, which are sets of instructions from a file, so you can readily reproduce your steps by running the script, without having to remember any manual things you did. Scripts also make you work completely portable and shareable with others so you can work in teams. As long as you have access to the raw data, you can email your bunch of scripts to your web-mail account, fly to the other side of the world and reproduce the exact steps without having to take your computer with you. This is also why it is nice to use software that is free and downloadable from the internet.

As I needed to frequently check files and automatically generate scripts to load the data, I wrote a tool (niftytools.exe, windows only) that has been a massive time saver. This checks the file, corrects dogy rows or removes them and generates the script to load data into SQL Server.

Here it is in use on the claims data,

and this is the script it generated for loading the data,

CREATE DATABASE HHP
USE HHP

CREATE TABLE Claims
(
    MemberID    int  -- varchar(8)  --integers starting with 0, could be text!
,    ProviderID    int  -- varchar(7)  --integers starting with 0, could be text!
,    Vendor    int  -- varchar(6)  --integers starting with 0, could be text!
,    PCP    int  -- varchar(5)  --integers starting with 0, could be text!
,    Year    varchar(2)
,    Specialty    varchar(25)
,    PlaceSvc    varchar(19)
,    PayDelay    varchar(4)
,    LengthOfStay    varchar(10)
,    DSFS    varchar(12)
,    PrimaryConditionGroup    varchar(8)
,    CharlsonIndex    varchar(3)
,    ProcedureGroup    varchar(4)
,    SupLOS    tinyint
)


BULK INSERT Claims
FROM 'E:\comps\hhp\raw data\HHP_release2\Claims.csv'
WITH
(
MAXERRORS = 0,
FIRSTROW = 2,
FIELDTERMINATOR = ',',
ROWTERMINATOR = '\n'
)

Note it has determined the maximum size the varchar (variable length character) fields need to be for optimal size, which makes the database as small as possible. It also determines what type of integer is most efficient. You will see that the tinyint is not coloured, which is because this is a relatively new datatype and the web page mentioned earlier that did this syntax highlighting obviously hasn't got it in its list of special words. I emailed the developer, so we shall see if it changes. Also the field name Year is highlighted because it as also a reserved word in SQL Server. When this happens it can be sometimes worthwhile renaming the field (just change the word Year to something else).

One other thing worth noting, there are comments starting with

-- integers starting with 0

It was noted on the HHP forum that the IDs in the submission entry example did not match the ids in the data. This is because they started with zeros, so whatever software was used to generate the submission example had treated them as integers and lost the leading zero. Thanks to this I was able to modify my nifty tool to generate warnings if this is the case, and if gives you an alternative if you want to treat the field as a character field, which will preserve the zero.

Leading zeros are sometimes important, as 003 and 03 are not the same, but would become the same if you treated them as integers. I learned this when I used to deal with address files that had postcodes beginning with 0.

Another 'bug' I discovered in my nifty tool was that "55+" passes the isNumeric test (try in an Excel macro - see the sample code at the top of this post) , when it is not really, and SQL Server won't accept it. This improved the tool again, and highlighted the fact that there was some data like this that would need to be closely examined later on.

Running the generated script will load the data and it can then be viewed (click on the image below to see how)

So we are on our way! All the data is in and we can take a peak at it to see what it looks like.

If you want to use niftytools to check and help load data then it is freely available. You will first need to download and install my Tiberius package, and then look in the program files folder to find it. Note you will need dot net 4 - see this blog post find out why.

Who needs the HERITAGE HEALTH PRIZE anyway?

2011-05-06T14:32:00.000-07:00

Just been playing around with this blogging software and discovered that it wants to pay me. Not sure how this works (I presume writing HERITAGE HEALTH PRIZE everywhere will help?) and whether posting this post is breaking the agreement I just agreed to (something to do with not tricking people to click on the adds). Anyway, lets see if I can top $3 million in 2 years! Click away guys and girls.

HERITAGE HEALTH PRIZE
DATA MINING
KAGGLE
HERITAGE HEALTH PRIZE ;-)

The number of visitors will obviously have to increase...

UPDATE

I had a visitor (the one person I told), as informed by the email below which came from the web counter you will see on the top right (which I discovered by clicking the 'next blog' button and saw on another blog).

Note in the email is a target add. Very clever - the person is from Australia and they want me to visit them!

Google just informed me my adds are up, with an email that says,

I guess google has some smart algorithm that determines relevance based on what I am gibbering on about. I think I will have to start talking about Data Mining and hopefully we will see some change. The first add I got was a bit worrying...

UPDATE 2.

While I type another visitor. Hi Dave! And another email. Gee, their algorithm must have cost millions to develop. Can't wait for my first 'visit the South Pole' email.

UPDATE 3 - relevant adds?

This add probably getting a bit more relevant, it's a competition and only 3 zeros missing, otherwise spot on!

And this double add has picked where 2 of my 3 visitors have come from (although probably more likely it knows where I am), although I've no idea what the first is is advertising, the second has picked up on the word heritage. Pretty impressive so far - but still the kitty is empty!

The day I led the HHP (briefly)

2011-05-06T14:06:00.000-07:00

The $3 million HHP competition started on Kaggle the other day. I was briefly leading (hope the leaderboard looks similar in 2 years time). Ecstatic in my efforts, I checked out who was online to brag to and discovered Dave.

'Hey Dave - winning $3 million and haven't even downloaded the data yet!'

Lesson 1: Dave is smarter then me!