Week 7 Project Anomaly Detection
John Neville
Prompt:
The following data set is the number of earthquakes per year magnitude 7.0 or greater. 1900-1998:
Use outlier detections techniques including box plot, or other techniques such as R open source from twitter (AnomalyDection package), tsoutliers package, h2o package to analyze outliers. Discuss your findings and suggested actionable decision.
The first piece of this exercise begins by downloading the data set that is provided. I have gone ahead and done that, which is read in below, and where we begin our redition of this exercise.
setwd('C:\\Users\\JP\\Downloads')Error in names(frame)[names(frame) == "x"] <- name :
names() applied to a non-vector
The working directory was changed to C:/Users/JP/Downloads inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.data<-read.csv(file='number-of-earthquakes-per-year-m.csv',header=T)Before we dive in, let’s take a look at the data qaulities and see what we imported.
str(data)'data.frame': 100 obs. of 2 variables:
$ Year : Factor w/ 100 levels "1900","1901",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Number.of.earthquakes.per.year.magnitude.7.0.or.greater..1900.1998: int 13 14 8 10 16 26 32 27 18 32 ...summary(data) Year Number.of.earthquakes.per.year.magnitude.7.0.or.greater..1900.1998
1900 : 1 Min. : 6.00
1901 : 1 1st Qu.:15.00
1902 : 1 Median :20.00
1903 : 1 Mean :20.02
1904 : 1 3rd Qu.:24.00
1905 : 1 Max. :41.00
(Other):94 NA's :1 plot(data)
The data appears to be very straight forward data over time, one record for each year, like a time series. In this case we are only looking for anomalies, and not trying to predict the future just yet. The most basic way to detect anomalies is to use a simple boxplot technique. This will help us to identify outliers. We do need to be careful here though. Our data is sequential and outliers are not always an anomaly.
names(data)[2]<-c('Earthquakes')
boxplot(data$Earthquakes)
There appear to be two obvious outliers in our dataset, but the boxplot does not consider the time series aspect, which should not be ignored.
The next approach we can take is using a package designed to identify anomalies for time series data. This is the package mentioned in the requirements above. Below we load the package and transform our data to meet the needs of the package. The Year field need to be converted to a date. We will also use ggplot2 to show plots of our data.
library(devtools)
library(AnomalyDetection)
library(ggplot2)
Year<-NULL
Year$year<-as.data.frame(as.POSIXlt(data$Year,format = '%Y'))
quakes<-as.data.frame(cbind(Year$year,data$Earthquakes))
names(quakes)<-c('Year','NumberOfQuakes')
ggplot(quakes, aes(Year, NumberOfQuakes)) + geom_line() + scale_x_datetime() + xlab("") + ylab("Number of Earthquakes")
It is good to see the data plotted as such, next we will use the package that we downloaded to identify the anomalies. To use the package it is helpful to go over all of the arguments that it takes. max.anoms is the argument that determines how many maximum anomalies there could be. alpha determines a level of significance for qaulifying anomalies. We can try a few of these out in different scenarios. Something that this package requires is a period of more than 1. We will just look at this in increments of 5 years for the purposes of getting results using this method.
m1<-AnomalyDetectionTs(quakes,max_anoms=.02,direction='both', plot=T,na.rm = T)
m1$plot
We found one this way, which is less than our boxplot. We can change a few things though. We can lower the alpha and see if any other anomalies show up.
m2<-AnomalyDetectionTs(quakes,max_anoms=.1,alpha =.25,direction='both', plot=T,na.rm = T)
m2$plot
At .25 we got a second, but its likely not significant at that level. This package appears to work alright, but not seeing the dates at the bottom is bothersome. The arguments are also very fragile. If we specified everything that matched our data set, we would get errors. Even the guide to this package on R-Bloggers had to change what they were doing to present their data in a way that satisfied their needs.
So far we know that the year with ~40 earthquakes is most likely an anomaly. Using the boxplot, and when we lower the significance, the year with just under 40 earthquakes also appears to be an anomaly. These two points in time are close together and therefore may not be anomalies, but outliers.
To dig further, we can try a few other methods. We will look at the package “tsoutliers” next, the second package mentioned in the prompt above. This one requires time series objects:
library(tsoutliers)
library(forecast)
quakesTS<-ts(quakes$NumberOfQuakes,start = c(1900))
quakesTSL<-log(quakesTS)
plot(quakesTS)
plot(quakesTSL)
We created a time series object to work on here and ran the auto arima function. We also took the log values of our original data to help with distances.
?tso()there is no package called <U+393C><U+3E31>anomalyDetection<U+393C><U+3E32>ts1<-tso(quakesTS,types = c("AO","LS","TC","IO",'SLS'),tsmethod = 'auto.arima')
ts1Series: quakesTS
Regression with ARIMA(1,0,0) errors
Coefficients:
ar1 intercept IO58
0.6138 19.9440 15.8789
s.e. 0.0789 1.4161 3.4839
sigma^2 estimated as 31.14: log likelihood=-309.9
AIC=627.8 AICc=628.22 BIC=638.22
Outliers:
type ind time coefhat tstat
1 IO 58 1957 15.88 4.558Interesting, here we have a third point just before 1960, the spike that can be seen in the plots above. Here it is idenfitied as an innovative outlier, which means that it does skew our results quite as much.
ts2<-tso(quakesTS,types = c("AO","LS","TC","IO",'SLS'),remove.method = 'en-masse',tsmethod = 'auto.arima')
ts2Series: quakesTS
Regression with ARIMA(1,0,0) errors
Coefficients:
ar1 intercept IO58
0.6138 19.9440 15.8789
s.e. 0.0789 1.4161 3.4839
sigma^2 estimated as 31.14: log likelihood=-309.9
AIC=627.8 AICc=628.22 BIC=638.22
Outliers:
type ind time coefhat tstat
1 IO 58 1957 15.88 4.558ts3<-tso(quakesTS,types = c("AO","LS","TC","IO",'SLS'),remove.method = 'bottom-up',tsmethod = 'auto.arima')
ts3Series: quakesTS
Regression with ARIMA(1,0,0) errors
Coefficients:
ar1 intercept IO58
0.6138 19.9440 15.8789
s.e. 0.0789 1.4161 3.4839
sigma^2 estimated as 31.14: log likelihood=-309.9
AIC=627.8 AICc=628.22 BIC=638.22
Outliers:
type ind time coefhat tstat
1 IO 58 1957 15.88 4.558ts4<-tso(quakesTSL,types = c("AO","LS","TC","IO",'SLS'),remove.method = 'bottom-up',args.tsmethod = list(1,0,1))
ts4Series: 1
ARIMA(1,0,1) with non-zero mean
Coefficients:
ar1 ma1 mean
0.8381 -0.4290 2.9087
s.e. 0.0774 0.1197 0.1057
sigma^2 estimated as 0.09805: log likelihood=-24.77
AIC=57.55 AICc=57.97 BIC=67.97
No outliers were detected.We try a few different options here, but this method still only shows the one outlier.
The final package that is recommended in the prompt is the ‘h20’ package. h20 is a deep learning algorithm and needs to have data presented in factor. We will add the
We were able to isolate points with the greatest error. Interestingly, we do not have the 1960 year in here, but rather, our points that the boxplot and the open source AnomalyDetectionTS dectected earlier, 1950 and 1943. To conclude this exercise, it is most likely that 1943 was a year that is likely an anomaly in terms of the number of earthquakes. 1950 is also likley an anomaly. 1957 could also be considered an anomaly, but has little impact on the data, as it spikes rather than misleads a trend.