Kaggle and R

Following up on last week’s post on doing a Kaggle competition, I then decided to see if I could explore the data more in R on my local desktop.  The competition is about analyzing a large group of house claims to give them a risk score.

I started the R studio to take a look at the initial data:

1 train <- read.csv("../Data/train.csv")
2 head(train)
3 summary(train)
4 
5 plot(train$Hazard)

A couple of things popped out.  All of the X variables look to be categorical.  Even the result “Hazard” is an integer with most of the values falling between 1 and 9.

With that in mind, I decided to split the dataset into two sections: the majority and the minority.

1 train.low <- subset(train, Hazard < 9)
2 train.high <- subset(train, Hazard >= 9)
3 
4 plot(train.low$Hazard)
5 plot(train.high$Hazard)

With the under as:

And the over 9 is like this

But I want to look at the Hazard score from a distribution point of view:

1 hazard.frame <- as.data.frame(table(train$Hazard))
2 colnames(hazard.frame) <- c("hazard","freq")
3 hist(hazard.frame$freq)
4 plot(x=hazard.frame$hazard, y=hazard.frame$freq)
5 plot(x=hazard.frame$hazard, log(y=hazard.frame$freq))
6 

The hist shows the left skew

and the log plot really shows the distribution

So there is clearly a diminishing return going on.   As of this writing, the leader is at 40%, which is about 20,400 of the 51,000 entries.   So if you could identify all of the ones correctly, you should get 37% of the way there.  To test it out, I submitted to Kaggle only ones:

LOL, so they must take away for incorrect answers as it is same as “all 0” benchmark.  So going back, I know that if I can predict the ones correctly and make a reasonable guess at the rest, I might be OK.   I went back and tuned my model some to get me out of the bottom 25% and then let it be.  I assume that there is something obvious/industry standard that I am missing because there are so many people between my position and the top 25%.

The Counted Part 2: Analysis Using R

Following up on this post last week on analyzing The Counted using F# and R, I decided to look a bit closer at the data.  In last week’s post, I had a data frame of all of the people killed by law enforcement collected by The Guardian for the 1st half of 2015.  Although interesting to looks at, I am not sure that the map tells us anything.  The data frame looks like this:

My first thought is to look at killing by population by US State.  Step #1 was to sum up the number of rows by state code:

1 the.counted.state <- data.frame(table(the.counted$state))
2 colnames(the.counted.state ) <- c("StateCode","NumberKilled")
3 summary(the.counted.state)
4 

I then brought in the latest population figures by state from the US Census:

1 state.population <- read.csv("http://www.census.gov/popest/data/state/asrh/2014/files/SCPRC-EST2014-18+POP-RES.csv")
2 state.population
3 

And finally I bought in a cross walk table of US State Codes (what The.Counted data is in and the US State Names, which is what US Census data is in)

1 state.crosswalk <- read.csv("http://www.fonz.net/blog/wp-content/uploads/2008/04/states.csv")
2 state.crosswalk
3 

I then merged all three data frames together using the state name and state code as the common key:

 1 state.population.2 <- state.population[c(5,6)]
 2 state.population.3 <- merge(x=state.population.2,
 3                             y=state.crosswalk,
 4                             by.x="NAME",
 5                             by.y="State")
 6 #The Counted With Population
 7 the.counted.state <- merge(x=the.counted.state,
 8                            y=state.population.3,
 9                             by.x="StateCode",
10                             by.y="Abbreviation")
11 

I then tried to add in a column that took the total number of killed individuals by the number of people in the state.

1 the.counted.state.2 <- the.counted.state
2 the.counted.state.2$KilledRatio <- the.counted.state.2$NumberKilled/the.counted.state.2$POPESTIMATE2014
3 

The problem became quickly obvious: there were not enough people in the numerator to make a meaningful straight division.  To compensate for this issue, I divided the number of people in each state by 10,000.  I also increased the kill ratio by a factor of 10 so that we have a scale between of 0 and 1 of .1 which is easily digestible.  Finally, I renamed the variable “Name” to “StateName” because my OCD couldn’t let such an affront to the naming gods go unpunished.

1 the.counted.state.3 <- the.counted.state
2 the.counted.state.3$AdjustedPopulation <- the.counted.state.2$POPESTIMATE2014/10000
3 the.counted.state.3$KilledRatio <- the.counted.state.3$NumberKilled/the.counted.state.3$AdjustedPopulation
4 the.counted.state.3$AdjKilledRatio <- the.counted.state.3$KilledRatio * 10
5 
6 names(the.counted.state.3)[names(the.counted.state.3)=="NAME"] <- "StateName"
7 the.counted.state.3$StateName <- tolower(the.counted.state.3$StateName)

With the data prepped, I created a choropleth to show the kill ratio by state using a gradient scale:

 1 choropleth <- merge(x=all.states, 
 2                y=the.counted.state.3, 
 3                sort = FALSE, 
 4                by.x = "region", 
 5                by.y = "StateName",
 6                all.x=TRUE)
 7 choropleth <- choropleth[order(choropleth$order), ]
 8 summary(choropleth)
 9 
10 qplot(long, lat, data = choropleth, group = group, fill = AdjKilledRatio,
11         geom = "polygon")

Note that I had to use the all.x=TRUE to account for the fact that South Dakota and Vermont did not have any killings so far in 2015.  This is equiv to a left outer join for you Sql folks.  On a side note, what’s up with Oklahoma? 

I then decided to bin the data into high,medium, and low categories.  Looking at the detail of the adjustedKillRatio, there seems to be some natural breaks around 10% and 20%:

1 the.counted.state.4$AdjKilledRatio
2 summary(the.counted.state.4$AdjKilledRatio)
3 

So I binned like that:

1 the.counted.state.4$KilledBin <- cut(the.counted.state.4$AdjKilledRatio,
2                                      breaks=seq(0,1,.1))
3 summary(the.counted.state.4$KilledBin)
4 

The problem with my code is that this gives me 10 bins and I only really need 3.  Fortunately, this stack overflow post helped me re-write the bin into 3 factors.  Note the Inf on the high side and the labels.

1 the.counted.state.4$KilledBin <- cut(the.counted.state.4$AdjKilledRatio,
2                                      breaks=c(seq(0,.2,.1),Inf),
3                                      labels=c("low","med","high"))
4 

And this gives me a pretty good distribution of bins:

With things binned up, I added another chiropleth and map:

 1 choropleth.2 <- merge(x=all.states, 
 2                     y=the.counted.state.4, 
 3                     sort = FALSE, 
 4                     by.x = "region", 
 5                     by.y = "StateName",
 6                     all.x=TRUE)
 7 choropleth.2 <- choropleth.2[order(choropleth.2$order), ]
 8 summary(choropleth.2)
 9 
10 qplot(long, 
11       lat, 
12       data = choropleth.2,
13       group = group,
14       fill = KilledBin,
15       geom = "polygon")

If your squint, it almost looks like a map of the civil war, no?

WCPSS Scores and Property Tax Valuations Using R

With all of the data gathered and organized I was ready to do some analytics using R.  The first thing I did was to load the four major datasets into R.

• NCScores is the original dataset that has the school score.  I had already done did an analysis on it here.
• SchoolValuation is the aggrgrate property values for each school as determined by scraping the Wake County Tax Website and Wake County School Assignment websites.  You can read how it was created here and here.
• SchoolNameMatch is a crosswalk table between the school name as found in the NCScores dataframe and the School Valuation dataframe.  You can read how it was created here
• WakeCountySchoolInfo is an export from WCPSS that was tossed around at open data day.

Step one was to reduce the North Carolina Scores data to only Wake County

1 #Create Wake County Scores From NC State Scores
2 WakeCountyScores <- NCScores[NCScores$District == 'Wake County Schools',]
3 

The next step was to add in the SchoolNameMatch so that we have the Tax Valuation School Name

1 #Join SchoolNameMatch to Wake County Scores
2 WakeCountyScores <- merge(x=WakeCountyScores, y=SchoolNameMatch, by.x="School", by.y="WCPSS")
3 

Interestingly, R is smart enough that the common field not duplicated, just the additional field(s) are added

The next step was to add in the Wake County Property Values, remove the Property field as it is no longer needed, and convert the TaxBase field from string to numeric

1 #Join Property Values
2 WakeCountyScores <- merge(x=WakeCountyScores, y=SchoolValuation, by.x="Property", by.y="SchooName")
3         
4 #Remove Property column
5 WakeCountyScores$Property = NULL
6         
7 #Turn tax base to numeric
8 WakeCountyScores$TaxBase <- as.numeric(WakeCountyScores$TaxBase)
9 

Eager to do an analysis, I pumped the data into a correlation

1 #Do a Correlation
2 cor(WakeCountyScores$TaxBase,WakeCountyScores$SchoolScore,use="complete")
3 

So clearly my expectations that property values track with FreeAndReducedLunch (.85 correlation) were not met.  I decided to use Practical Data Science with R Chapter 3 (Exploring Data)  as a guide to better understand the dataset.

1 #Practical Data Science With R, Chapter3
2 summary(WakeCountyScores)
3 summary(WakeCountyScores$TaxBase)
4 

So there is quite a range in tax base!  The next task was to use some graphs to explore the data.  I added in ggplot2

and followed the books example for a histogram.  I started with score and it comes out as expected.  I then tried a historgram on TaxBase and had to tinker with the binwidth to make a meaningful chart:

1 #Historgrams
2 ggplot(WakeCountyScores) + geom_histogram(aes(x=SchoolScore),binwidth=5,fill="gray")      
3 ggplot(WakeCountyScores) + geom_histogram(aes(x=TaxBase),binwidth=10000,fill="gray")
4 #Ooops
5 ggplot(WakeCountyScores) + geom_histogram(aes(x=TaxBase),binwidth=5000000,fill="gray")
6 

The book then moves to an example studying income, which is directly analogous to the TaxBase so I followed it very closely.  The next graph were some density graphs.  Note the second one is a logarithmic one:

1 #Density
2 library(scales)
3 ggplot(WakeCountyScores) + geom_density(aes(x=TaxBase)) + scale_x_continuous(labels=dollar)
4 ggplot(WakeCountyScores) + geom_density(aes(x=TaxBase)) + scale_x_log10(labels=dollar) + annotation_logticks(sides="bt")
5 

So kinda interesting that most schools cluster in terms of their tax base, but because there is such a wide range with a majority clustered to the low end, the logarithmic curve is much more revealing.

The book then moved into showing the relationship between two variables.  In this case, SchoolScore as the Y variable and TaxBase as the X variable:

1 #Relationship between TaxBase and Scores
2 ggplot(WakeCountyScores, aes(x=TaxBase, y=SchoolScore)) + geom_point()
3 ggplot(WakeCountyScores, aes(x=TaxBase, y=SchoolScore)) + geom_point() + stat_smooth(method="lm")
4 ggplot(WakeCountyScores, aes(x=TaxBase, y=SchoolScore)) + geom_point() + geom_smooth()
5 

So what is interesting is that there does not seem to be a strong relationship between scores and tax base.  There looks like an equal number of schools both below the score curve than above it.  Note that using a smoothing curve is much better than the linear fit curve in showing the relationship of scores to tax base.  You can see the dip in the lower quartile and the increase at the tail.  It makes sense that the higher tax base shows an increase in scores, but what’s up with that dip?

Finally, the same data is shown using a hax chart

1 library(hexbin)
2 ggplot(WakeCountyScores, aes(x=TaxBase, y=SchoolScore)) + geom_hex(binwidth=c(100000000,5)) + geom_smooth(color="white",se=F)
3 

So taking a step back, it is clear that there is a weakness in this analysis.  Some schools have thousands of students, some schools have a couple hundred.  (high schools versus elementary students). Using the absolute dollars from the tax valuation is misleading.  What we really need is revenue per student.  Going back to the SchoolInfo dataframe, I added it in and pulled a student count column.

1 WakeCountyScores <- merge(x=WakeCountyScores, y=WakeCountySchoolInfo, by.x="School", by.y="School.Name")      
2 names(WakeCountyScores)[names(WakeCountyScores)=="School.Membership.2013.14..ADM..Mo2."] <- "StudentCount"
3 WakeCountyScores$StudentCount <- as.numeric(WakeCountyScores$StudentCount)
4         
5 WakeCountyScores["TaxBasePerStudent"] <- WakeCountyScores$TaxBase/WakeCountyScores$StudentCount
6 summary(WakeCountyScores$TaxBasePerStudent)
7 

Interestingly, the number of records in the base frame dropped from 166 to 152, which means that perhaps we need a second mapping table.  In any event, you can see that the average tax base is $6.5 million with a max of $114 Million.  Quite a range!

Going back to the point and hex graphs

1 ggplot(WakeCountyScores, aes(x=TaxBasePerStudent, y=SchoolScore)) + geom_point() + geom_smooth()
2 ggplot(WakeCountyScores, aes(x=TaxBasePerStudent, y=SchoolScore)) + geom_hex(binwidth=c(25000000,5)) + geom_smooth(color="white",se=F)
3 

There is some interesting going on.  First, the initial conclusion that a higher tax base leads to a gradual increase in scores is wrong once you move from total tax to tax per student.

Also, note the significant drop in school scores once you move away from the lowest tax base schools, the recovery, and then the drop again.  From a real estate perspective, these charts suggest that the marginal value of a really expensive or really inexpensive house in Wake County is not worth it (at least in terms of where you send you kids), and there is a sweet spot of value above a certain price point. 

You can find the gist here and the repo is here.

Some lessons I learned in doing this exercise:

• Some records got dropped between the scores dataframe and the info dataframe -> so there needs to be another mapping table
• Make the tax base in millions
• What’s up with that school with 114 million per student?
• An interesting question is location of dollars to school compared to tax base.  I wonder if that is on WCPSS somewhere.  Hummm…
• You can’t use the tick(‘) notation, which means you do a lot of overwriting of dataframes.  This can be a costly and expensive feature of the language.  It is much better to assume immutably, even if you clutter up your data window.

As a final note, I was using the Console window b/c that is what the intro books do.  This is a huge mistake in R Studio.  It is much better to create a script and send the results to the console

So that you can make changes and run things again.  It is a cheap way of avoid head-scratching bugs…