Forem: Kiefer Smith

Clustering in R

Kiefer Smith — Tue, 13 Jun 2017 10:02:26 +0000

Originally published on Wordpress.

Hello, everyone! I've been meaning to get a new blog post out for the past couple of weeks. During that time I've been messing around with clustering. Clustering, or cluster analysis, is a method of data mining that groups similar observations together. Classification and clustering are quite alike, but clustering is more concerned with exploration than actual results.

Note: This post is far from an exhaustive look at all clustering has to offer. Check out this guide for more. I am reading Data Mining by Aggarwal presently, which is very informative.

data("iris")
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

For simplicity, we'll use the iris dataset. We're going to try to use the numeric data to correctly group the observations by species. There are 50 of each species in the dataset, so ideally we would end up with three clusters of 50.

library(ggplot2)
ggplot() +
geom_point(aes(iris$Sepal.Length, iris$Sepal.Width, col = iris$Species))

As you can see, there is already some groupings present. Let's use hierarcical clustering first.

iris2 <- iris[,c(1:4)]  #not going to use the `Species` column.
medians <- apply(iris2, 2, median)
mads <- apply(iris2,2,mad)
iris3 <- scale(iris2, center = medians, scale = mads)

dist <- dist(iris3)
hclust <- hclust(dist, method = 'median')  #there are several methods for hclust
cut <- cutree(hclust, 3)

table(cut)

## cut
##  1  2  3 
## 49 34 67

iris <- cbind(iris, cut)
iris$cut <- factor(iris$cut)
levels(iris$cut) <- c('setosa', 'versicolor', 'virginica')
err <- iris$Species == iris$cut
table(err)

## err
## FALSE  TRUE 
##    38   112

ggplot() +
geom_point(aes(iris2$Sepal.Length, iris2$Sepal.Width, col = iris$cut))

Nice groupings here, but it looks like the algorithm has some trouble determining between versicolor and virginica. If we used this information to classify the observations, we'd get an error rate of about .25. Let's try another clustering technique: DBSCAN.

library(dbscan)

db <- dbscan(iris3, eps = 1, minPts = 20)
table(db$cluster)

## 
##  0  1  2 
##  5 48 97

iris2 <- cbind(iris2, db$cluster)

iris2$cluster <- factor(iris2$`db$cluster`)

ggplot() +
geom_point(aes(iris2$Sepal.Length, iris2$Sepal.Width, col = iris2$cluster))

DBSCAN classifies points into three different categories: core, border, and noise points on the basis of density. Thus, the versicolor/ virginica cluster is treated as one group. Since our data is not structured in such a way that density is meaningful, DBSCAN is probably not a wise choice here.

Let's look at one last algo: the ROCK. No, not that ROCK.

library(cba)

distm <- as.matrix(dist)
rock <- rockCluster(distm, 3, theta = .02)

## Clustering:
## computing distances ...
## computing links ...
## computing clusters ...

iris$rock <- rock$cl
levels(iris$rock) <- c('setosa', 'versicolor', 'virginica')

ggplot() +
geom_point(aes(iris2$Sepal.Length, iris2$Sepal.Width, col = rock$cl))

err <- (iris$Species == iris$rock)
table(err)

## err
## FALSE  TRUE 
##    24   126

While it may not look like it, the ROCK seems to do the best job at determining clusters in this data - the error rate dropped to 16%. Typically this method is reserved for categorical data, but since we used dist it shouldn't cause any problems.

I have been working on a project using some of these (and similar) data mining procedures to explore spatial data and search for distinct groups. While clustering the iris data may not be all that meaningful, I think it is illustrative of the power of clustering. I have yet to try higher-dimension clustering techniques, which might perform even better.

Have any comments? Questions? Suggestions for future posts? I am always happy to hear from readers; please contact me!

Happy clustering,

Kiefer

Which linear model is best?

Kiefer Smith — Wed, 10 May 2017 16:37:16 +0000

Originally published on: Wordpress.

Recently I have been working on a Kaggle competition where participants are tasked with predicting Russian housing prices. In developing a model for the challenge, I came across a few methods for selecting the best regression
model for a given dataset. Let's load up some data and take a look.

library(ISLR)
set.seed(123)
swiss <- data.frame(datasets::swiss)
dim(swiss)

## [1] 47  6

head(swiss)

##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6

This data provides data on Swiss fertility and some socioeconomic
indicators. Suppose we want to predict fertility rate. Each observation
is as a percentage. In order to assess our prediction ability, create
test and training data sets.

index <- sample(nrow(swiss), 5)
train <- swiss[-index,]
test <- swiss[index,]

Next, have a quick look at the data.

par(mfrow=c(2,2))
plot(train$Education, train$Fertility)
plot(train$Catholic, train$Fertility)
plot(train$Infant.Mortality, train$Fertility)
hist(train$Fertility)

Looks like there could be some interesting relationships here. The
following block of code will take a model formula (or model matrix) and
search for the best combination of predictors.

library(leaps)
leap <- regsubsets(Fertility~., data = train, nbest = 10)
leapsum <- summary(leap)
plot(leap, scale = 'adjr2')

According to the adjusted R-squared value (larger is better), the best
two models are: those with all predictors and all predictors less
Examination. Both have adjusted R-squared values around .69 - a decent fit. Fit these models so we can evaluate them further.

swisslm <- lm(Fertility~., data = train)
swisslm2 <- lm(Fertility~.-Examination, data = train)
#use summary() for a more detailed look.

First we'll compute the mean squared error (MSE).

mean((train$Fertility-predict(swisslm, train))[index]^2)

## [1] 44.21879

mean((train$Fertility-predict(swisslm2, train))[index]^2)

## [1] 36.4982

Looks like the model with all predictors is marginally better. We're
looking for a low value of MSE here. We can also look at the Akaike
information criterion (AIC) for more information. Lower is better here
as well.

library(MASS)
step1 <- stepAIC(swisslm, direction = "both")

step1$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## Fertility ~ Agriculture + Examination + Education + Catholic + 
##     Infant.Mortality
## 
## Final Model:
## Fertility ~ Agriculture + Education + Catholic + Infant.Mortality
## 
## 
##            Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                                  36   1746.780 168.5701
## 2 - Examination  1 53.77608        37   1800.556 167.8436

Here, the model without Examination scores lower than the full
model. It seems that both models are evenly matched, though I might be
inclined to use the model without Examination.

Since we sampled our original data to make the train and test datasets,
the difference in these tests is subject to change based on the training
data used. I encourage anyone who wants to test themselves to change the
set.seed at the beginning and evaluate the above results again. Even
better: use your own data!

If you learned something or have questions feel free to leave a comment
or shoot me an email.

Happy modeling,

Kiefer

Using R as a Geographic Information System

Kiefer Smith — Sat, 22 Apr 2017 00:00:00 +0000

Originally published on: realdataweb.wordpress.com

In real estate, spatial data is the name of the game. Countless programs in other domains utilize the power of this data, which is becoming more prevalent by the day.

In this post I will go over a few simple, but powerful tools to get you started using using geographic information in R.

##First, some libraries##
#install.packages('GISTools', dependencies = T)
library(GISTools)

GISTools provides an easy-to-use method for creating shading schemes and choropleth maps. Some of you may have heard of the sp package, which adds numerous spatial classes to the mix. There are also functions for analysis and making things look nice.

Let's get rolling: source the vulgaris dataset, which contains location information for Syringa Vulgaris (the Lilac) observation stations and US states. This code plots the states and vulgaris points.

data("vulgaris")          #load data
par = (mar = c(2,0,0,0))  #set margins of plot area
plot(us_states)
plot(vulgaris, add = T, pch = 20)

One thing to note here is the structure of these objects. us_states is a SpatialPolygonsDataFrame, which stores information for plotting shapes (like a shapefile) within its attributes. vulgaris by contrast is a SpatialPointsDataFrame, which contains data for plotting individual points. Much like a data.frame and $, these objects harbor information that can be accessed via @.

kable(head(vulgaris@data))

	Station	Year	Type	Leaf	Bloom	Station.Name	State.Prov	Lat	Long	Elev
3695	61689	1965	Vulgaris	114	136	COVENTRY	CT	41.8	-72.35	146
3696	61689	1966	Vulgaris	122	146	COVENTRY	CT	41.8	-72.35	146
3697	61689	1967	Vulgaris	104	156	COVENTRY	CT	41.8	-72.35	146
3698	61689	1968	Vulgaris	97	134	COVENTRY	CT	41.8	-72.35	146
3699	61689	1969	Vulgaris	114	138	COVENTRY	CT	41.8	-72.35	146
3700	61689	1970	Vulgaris	111	135	COVENTRY	CT	41.8	-72.35	146

Let's take a look at some functions that use this data.

newVulgaris <- gIntersection(us_states, vulgaris, byid = T)
kable(head(data.frame(newVulgaris)))

	x	y
3 4896	-67.65	44.65
3 4897	-67.65	44.65
3 4898	-67.65	44.65
3 4899	-67.65	44.65
3 4900	-67.65	44.65
3 4901	-67.65	44.65

gIntersection, as you may have guessed from the name, returns the intersection of two spatial objects. In this case, we are given the points from vulgaris that are within us_states. However, the rest of the vulgaris data has been stripped from the resulting object. We've got to jump through a couple of hoops to get that information back.

newVulgaris <- data.frame(newVulgaris)
tmp <- rownames(newVulgaris)
tmp <- strsplit(tmp, " ")
tmp <- (sapply(tmp, "[[", 2))
tmp <- as.numeric(tmp)
vdf <- data.frame(vulgaris)
newVulgaris <- subset(vdf, row.names(vdf) %in% tmp)

	Station	Year	Type	Leaf	Bloom	Station.Name	State.Prov	Lat	Long	Elev	Long.1	Lat.1	optional
3695	61689	1965	Vulgaris	114	136	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE
3696	61689	1966	Vulgaris	122	146	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE
3697	61689	1967	Vulgaris	104	156	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE
3698	61689	1968	Vulgaris	97	134	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE
3699	61689	1969	Vulgaris	114	138	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE
3700	61689	1970	Vulgaris	111	135	COVENTRY	CT	41.8	-72.35	146	-72.35	41.8	TRUE

Look familiar? Now we've got a data frame with the clipped vulgaris values and original data preserved.

vulgarisSpatial <- SpatialPointsDataFrame(data.frame(newVulgaris$Long, newVulgaris$Lat), newVulgaris, proj4string = CRS(proj4string(vulgaris)), bbox = vulgaris@bbox)

After storing our clipped data frame as a SpatialPointsDataFrame, we can again make use of it - in this case we add a shading scheme to the vulgaris points.

shades <- auto.shading(vulgarisSpatial@data$Elev, n = 7, cols = brewer.pal(7, 'YlOrRd'))  #Check cutter arg for more ways to create breaks.
shades$cols <- add.alpha(shades$cols, .5)
plot(us_states)
choropleth(vulgarisSpatial, vulgarisSpatial$Elev,shading = shades, add = T, pch = 20)

Colors are pretty, but what do they mean? Let's add a legend.

us_states@bbox #Get us_states bounding box coordinates.

##           min       max
## r1 -124.73142 -66.96985
## r2   24.95597  49.37173

plot(us_states)
choropleth(vulgarisSpatial, vulgarisSpatial$Elev,shading = shades, add = T, pch = 20)
par(xpd=TRUE)           #Allow plotting outside of plot area.
choro.legend(-124, 30, shades, cex = .75, title = "Elevation in Meters") # Plot legend in bottom left.  Takes standard legend() params.

It looks like there's a lot going on in the Northeastern states. For a closer look, create another clipping (like above) and plot it. Using the structure below, we can create a selection vector. I have hidden the full code since it is repetitive (check GitHub for the full code.)

index <- us_states$STATE_NAME == "Pennsylvania"
         '...'

plot(us_states[index,])
choropleth(vulgarisNE, vulgarisNE$Elev,shading = shades, add = T, pch = 20)
par(xpd = T)
choro.legend(-73, 39.75, shades, cex = .75, title = "Elevation in Meters")

Hopefully this has been a useful introduction (or refresher) on spatial data. I always learn a lot in the process of writing these posts. If you have any ideas or suggestions please leave a comment or feel free to contact me!

Happy mapping,

Kiefer