Detecting Anomolous Events in Washington DC Bike Sharing Dataset

Shin Chin

Posted on Apr 22, 2016

Contributed by Shin Chin. He is one of the remove students in the NYC Data Science Academy 12 week full time Data Science Online Bootcamp program taking place between January 11th to April 1st, 2016. This post is based on his data visualization project.

We want to detect anomolous events in the Washington DC Bike Sharing dataset.

Data Set Information:

Bike sharing systems are new generation of traditional bike rentals where whole process from membership, rental and return back has become automatic. Through these systems, user is able to easily rent a bike from a particular position and return back at another position. Currently, there are about over 500 bike-sharing programs around the world which is composed of over 500 thousands bicycles. Today, there exists great interest in these systems due to their important role in traffic, environmental and health issues.

Apart from interesting real world applications of bike sharing systems, the characteristics of data being generated by these systems make them attractive for the research. Opposed to other transport services such as bus or subway, the duration of travel, departure and arrival position is explicitly recorded in these systems. This feature turns bike sharing system into a virtual sensor network that can be used for sensing mobility in the city. Hence, it is expected that most of important events in the city could be detected via monitoring these data.

Attribute Information:

Both hour.csv and day.csv have the following fields, except hr which is not available in day.csv

instant: record index
dteday : date
season : season (1:springer, 2:summer, 3:fall, 4:winter)
yr : year (0: 2011, 1:2012)
mnth : month ( 1 to 12)
hr : hour (0 to 23)
holiday : weather day is holiday or not (extracted from [Web Link])
weekday : day of the week
workingday : if day is neither weekend nor holiday is 1, otherwise is 0.
weathersit :
1: Clear, Few clouds, Partly cloudy, Partly cloudy
2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist
3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds
4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog
temp : Normalized temperature in Celsius. The values are derived via (t-t_min)/(t_max-t_min), t_min=-8, t_max=+39 (only in hourly scale)
atemp: Normalized feeling temperature in Celsius. The values are derived via (t-t_min)/(t_max-t_min), t_min=-16, t_max=+50 (only in hourly scale)
hum: Normalized humidity. The values are divided to 100 (max)
windspeed: Normalized wind speed. The values are divided to 67 (max)
casual: count of casual users
registered: count of registered users
cnt: count of total rental bikes including both casual and registered

First we load the required libraries as well as read in the “day” data from the bike sharing dataset.

rm(list=ls())

library('ggplot2')

## Warning: package 'ggplot2' was built under R version 3.2.4

library('dplyr')

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

day <- read.csv("~/NYCDataScience_course/Final_Project/Bike-Sharing-Dataset/day.csv")
#View(day)

Next we set season, yr, mnth, holiday, weekday, workingday and weathersit as factors.

day$season = as.factor(day$season)
day$yr = as.factor(day$yr)
day$mnth = as.factor(day$mnth)
day$holiday = as.factor(day$holiday)
day$weekday = as.factor(day$weekday)
day$workingday = as.factor(day$workingday)
day$weathersit = as.factor(day$weathersit)

We now do some visualizations by plotting various variables against each other with points color coded according the magintude of cnt.

scale_color_gradient(low = "#132B43", high = "#56B1F7", space = "Lab", na.value = "grey50", guide = "colourbar")

## <ScaleContinuous>
##  Range:  
##  Limits:    0 --    1

d1 <- qplot(atemp, hum, data=day, color=cnt)
d1 + scale_colour_gradient(limits=c(22, 8714))

d2 <- qplot(temp, hum, data=day, color=cnt)
d2 + scale_colour_gradient(limits=c(22, 8714))

d3 <- qplot(windspeed, temp, data=day, color=cnt)
d3 + scale_colour_gradient(limits=c(22, 8714))

d3 <- qplot(windspeed, hum, data=day, color=cnt)
d3 + scale_colour_gradient(limits=c(22, 8714))

As you can see from the plots, cnt increases with increasing temp as well as decreasing windspeed. Humdity may or may not have some effect on cnt.

Next we plot hitograms of cnt versus various variables.

qplot(cnt, data=day, binwidth=100) + facet_wrap(~weathersit)

by_weather <-group_by(day, weathersit)
summarize(by_weather, total=sum(cnt))

## Source: local data frame [3 x 2]
## 
##   weathersit   total
##       (fctr)   (int)
## 1          1 2257952
## 2          2  996858
## 3          3   37869

qplot(cnt, data=day, binwidth=100) + facet_wrap(~holiday)

by_holiday <-group_by(day, holiday)
summarize(by_holiday, total=sum(cnt))

## Source: local data frame [2 x 2]
## 
##   holiday   total
##    (fctr)   (int)
## 1       0 3214244
## 2       1   78435

qplot(cnt, data=day, binwidth=100) + facet_wrap(~workingday)

by_workingday <-group_by(day, workingday)
summarize(by_workingday, total=sum(cnt))

## Source: local data frame [2 x 2]
## 
##   workingday   total
##       (fctr)   (int)
## 1          0 1000269
## 2          1 2292410

qplot(cnt, data=day, binwidth=100) + facet_wrap(~season)

by_season <-group_by(day, season)
summarize(by_season, total=sum(cnt))

## Source: local data frame [4 x 2]
## 
##   season   total
##   (fctr)   (int)
## 1      1  471348
## 2      2  918589
## 3      3 1061129
## 4      4  841613

qplot(cnt, data=day, binwidth=100) + facet_wrap(~weekday)

by_weekday <-group_by(day, weekday)
summarize(by_weekday, total=sum(cnt))

## Source: local data frame [7 x 2]
## 
##   weekday  total
##    (fctr)  (int)
## 1       0 444027
## 2       1 455503
## 3       2 469109
## 4       3 473048
## 5       4 485395
## 6       5 487790
## 7       6 477807

cnt is greatest for weathersit = 1 and then 2 with low values for 3.
cnt is much larger when holiday = 0 than 1
cnt is much larger for workingday = 1 than 0
cnt is largest for season = 3 followed by 2, 4 and 1.
cnt is roughly similar for all weekday.

We now perform multiple linear regressiom on the “day”" data.

model.saturated = lm(cnt ~ . -instant -dteday -casual -registered, data = day)

summary(model.saturated)

## 
## Call:
## lm(formula = cnt ~ . - instant - dteday - casual - registered, 
##     data = day)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3944.7  -348.2    63.8   457.4  2912.7 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1485.84     239.75   6.198 9.77e-10 ***
## season2       884.71     179.49   4.929 1.03e-06 ***
## season3       832.70     213.13   3.907 0.000102 ***
## season4      1575.35     181.00   8.704  < 2e-16 ***
## yr1          2019.74      58.22  34.691  < 2e-16 ***
## mnth2         131.03     143.78   0.911 0.362443    
## mnth3         542.83     165.43   3.281 0.001085 ** 
## mnth4         451.17     247.57   1.822 0.068820 .  
## mnth5         735.51     267.63   2.748 0.006145 ** 
## mnth6         515.40     282.41   1.825 0.068423 .  
## mnth7          30.80     313.82   0.098 0.921854    
## mnth8         444.95     303.17   1.468 0.142639    
## mnth9        1004.17     265.12   3.788 0.000165 ***
## mnth10        519.67     241.55   2.151 0.031787 *  
## mnth11       -116.69     230.78  -0.506 0.613257    
## mnth12        -89.59     182.21  -0.492 0.623098    
## holiday1     -589.70     180.36  -3.270 0.001130 ** 
## weekday1      212.05     109.49   1.937 0.053187 .  
## weekday2      309.53     107.13   2.889 0.003982 ** 
## weekday3      381.36     107.48   3.548 0.000414 ***
## weekday4      386.34     107.53   3.593 0.000350 ***
## weekday5      436.98     107.44   4.067 5.30e-05 ***
## weekday6      440.46     106.56   4.133 4.01e-05 ***
## workingday1       NA         NA      NA       NA    
## weathersit2  -462.54      77.09  -6.000 3.16e-09 ***
## weathersit3 -1965.09     197.05  -9.972  < 2e-16 ***
## temp         2855.01    1398.16   2.042 0.041526 *  
## atemp        1786.16    1462.12   1.222 0.222261    
## hum         -1535.47     292.45  -5.250 2.01e-07 ***
## windspeed   -2823.30     414.55  -6.810 2.09e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 769.2 on 702 degrees of freedom
## Multiple R-squared:  0.8484, Adjusted R-squared:  0.8423 
## F-statistic: 140.3 on 28 and 702 DF,  p-value: < 2.2e-16

plot(model.saturated) #Assessing the assumptions of the model.

library(car) #Companion to applied regression.
influencePlot(model.saturated)

##       StudRes        Hat     CookD
## 595  2.112057 0.64959535 0.5326909
## 669 -5.336784 0.04056466 0.1998979

alias(model.saturated)

## Model :
## cnt ~ (instant + dteday + season + yr + mnth + holiday + weekday + 
##     workingday + weathersit + temp + atemp + hum + windspeed + 
##     casual + registered) - instant - dteday - casual - registered
## 
## Complete :
##             (Intercept) season2 season3 season4 yr1 mnth2 mnth3 mnth4
## workingday1  0           0       0       0       0   0     0     0   
##             mnth5 mnth6 mnth7 mnth8 mnth9 mnth10 mnth11 mnth12 holiday1
## workingday1  0     0     0     0     0     0      0      0     -1      
##             weekday1 weekday2 weekday3 weekday4 weekday5 weekday6
## workingday1  1        1        1        1        1        0      
##             weathersit2 weathersit3 temp atemp hum windspeed
## workingday1  0           0           0    0     0   0

model.saturated seems like a good candidate for linear regression from the above plots.

As you can see from the alias values, workinday1 is perfectly collinear with holiday1.
Hence remove workingday variable in the next model2.

#vif(model.saturated) #Assessing the variance inflation factors for the variables
#in our model.

model2 = lm(cnt ~ . -instant -dteday -workingday -atemp -casual -registered, data = day)

summary(model2) #R^2 adjusted went up, model still significant, etc.

## 
## Call:
## lm(formula = cnt ~ . - instant - dteday - workingday - atemp - 
##     casual - registered, data = day)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3960.9  -350.9    74.1   456.0  2919.9 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1543.552    235.129   6.565 1.01e-10 ***
## season2       889.302    179.516   4.954 9.12e-07 ***
## season3       832.236    213.204   3.903 0.000104 ***
## season4      1578.947    181.040   8.722  < 2e-16 ***
## yr1          2018.063     58.225  34.660  < 2e-16 ***
## mnth2         136.855    143.747   0.952 0.341396    
## mnth3         545.132    165.480   3.294 0.001036 ** 
## mnth4         456.494    247.615   1.844 0.065667 .  
## mnth5         723.520    267.541   2.704 0.007010 ** 
## mnth6         490.552    281.776   1.741 0.082133 .  
## mnth7           8.404    313.395   0.027 0.978613    
## mnth8         404.912    301.494   1.343 0.179700    
## mnth9         983.948    264.698   3.717 0.000217 ***
## mnth10        520.937    241.636   2.156 0.031432 *  
## mnth11       -111.362    230.816  -0.482 0.629621    
## mnth12        -84.389    182.229  -0.463 0.643439    
## holiday1     -603.605    180.066  -3.352 0.000845 ***
## weekday1      214.877    109.508   1.962 0.050133 .  
## weekday2      309.132    107.171   2.884 0.004041 ** 
## weekday3      377.407    107.467   3.512 0.000473 ***
## weekday4      385.206    107.562   3.581 0.000366 ***
## weekday5      428.604    107.258   3.996 7.12e-05 ***
## weekday6      438.699    106.590   4.116 4.32e-05 ***
## weathersit2  -465.202     77.083  -6.035 2.57e-09 ***
## weathersit3 -1981.357    196.670 -10.075  < 2e-16 ***
## temp         4487.305    411.838  10.896  < 2e-16 ***
## hum         -1518.178    292.208  -5.196 2.68e-07 ***
## windspeed   -2925.438    406.175  -7.202 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 769.5 on 703 degrees of freedom
## Multiple R-squared:  0.848,  Adjusted R-squared:  0.8422 
## F-statistic: 145.3 on 27 and 703 DF,  p-value: < 2.2e-16

plot(model2) #No overt additional violations.

influencePlot(model2) #No overt additional violations; Hawaii actually lowers

##       StudRes        Hat     CookD
## 69  -1.210739 0.16828646 0.1028882
## 668 -4.756498 0.07702696 0.2557740
## 669 -5.356592 0.04026858 0.2033892

#its hat value (leverage).

vif(model2) #VIFs all decrease.

##                  GVIF Df GVIF^(1/(2*Df))
## season     169.713093  3        2.352985
## yr           1.046249  1        1.022863
## mnth       391.715048 11        1.311784
## holiday      1.116829  1        1.056801
## weekday      1.153293  6        1.011956
## weathersit   1.886133  2        1.171907
## temp         7.006223  1        2.646927
## hum          2.135348  1        1.461283
## windspeed    1.221500  1        1.105215

anova(model2, model.saturated)

## Analysis of Variance Table
## 
## Model 1: cnt ~ (instant + dteday + season + yr + mnth + holiday + weekday + 
##     workingday + weathersit + temp + atemp + hum + windspeed + 
##     casual + registered) - instant - dteday - workingday - atemp - 
##     casual - registered
## Model 2: cnt ~ (instant + dteday + season + yr + mnth + holiday + weekday + 
##     workingday + weathersit + temp + atemp + hum + windspeed + 
##     casual + registered) - instant - dteday - casual - registered
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1    703 416284778                           
## 2    702 415401688  1    883090 1.4924 0.2223

model2 seems like a good candidate for linear regression from the above plots.
The VIF values for all variables are low indicating that collinearity is not a problem.
From analysis of variance, seems like model2 and model.saturated are significantly different.
We now have a reduced model where we remove the hum variable.

model.reduced = lm(cnt ~ . -instant -dteday -workingday -atemp -hum -casual -registered, data = day)

anova(model.reduced, model.saturated)

## Analysis of Variance Table
## 
## Model 1: cnt ~ (instant + dteday + season + yr + mnth + holiday + weekday + 
##     workingday + weathersit + temp + atemp + hum + windspeed + 
##     casual + registered) - instant - dteday - workingday - atemp - 
##     hum - casual - registered
## Model 2: cnt ~ (instant + dteday + season + yr + mnth + holiday + weekday + 
##     workingday + weathersit + temp + atemp + hum + windspeed + 
##     casual + registered) - instant - dteday - casual - registered
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1    704 432269190                                  
## 2    702 415401688  2  16867502 14.252 8.563e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.reduced)

## 
## Call:
## lm(formula = cnt ~ . - instant - dteday - workingday - atemp - 
##     hum - casual - registered, data = day)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4189.5  -372.9    67.2   481.8  2893.1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   715.03     175.95   4.064 5.37e-05 ***
## season2       848.08     182.62   4.644 4.08e-06 ***
## season3       841.83     217.10   3.878 0.000115 ***
## season4      1585.86     184.35   8.603  < 2e-16 ***
## yr1          2065.21      58.57  35.263  < 2e-16 ***
## mnth2         180.92     146.12   1.238 0.216071    
## mnth3         638.19     167.52   3.810 0.000151 ***
## mnth4         606.99     250.41   2.424 0.015603 *  
## mnth5         820.13     271.78   3.018 0.002639 ** 
## mnth6         749.76     282.40   2.655 0.008110 ** 
## mnth7         242.98     315.80   0.769 0.441899    
## mnth8         575.45     305.18   1.886 0.059762 .  
## mnth9        1033.99     269.36   3.839 0.000135 ***
## mnth10        541.10     246.02   2.199 0.028175 *  
## mnth11       -109.31     235.04  -0.465 0.642021    
## mnth12       -119.21     185.44  -0.643 0.520540    
## holiday1     -585.60     183.33  -3.194 0.001464 ** 
## weekday1      222.82     111.50   1.998 0.046061 *  
## weekday2      328.50     109.07   3.012 0.002688 ** 
## weekday3      394.44     109.38   3.606 0.000333 ***
## weekday4      440.09     109.00   4.038 5.99e-05 ***
## weekday5      476.69     108.81   4.381 1.36e-05 ***
## weekday6      467.50     108.39   4.313 1.84e-05 ***
## weathersit2  -702.78      63.19 -11.122  < 2e-16 ***
## weathersit3 -2422.60     180.63 -13.412  < 2e-16 ***
## temp         3923.88     404.57   9.699  < 2e-16 ***
## windspeed   -2328.09     396.69  -5.869 6.76e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 783.6 on 704 degrees of freedom
## Multiple R-squared:  0.8422, Adjusted R-squared:  0.8364 
## F-statistic: 144.5 on 26 and 704 DF,  p-value: < 2.2e-16

plot(model.reduced)