6
$\begingroup$

Sample data

dat <- structure(list(yld.harvest = c(1800, 2400, 2000, 2400, 2160, 
        2400, 2400, 2250, 2400, 2280, 2400, 3120, 3300, 3300,                                     3000, 3000, 2400, 2700, 3000), 
           year = c(1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 
                    2004, 2005, 2006, 2007, 2008, 2009, 2011, 2012, 2013, 
                    2014, 2015)), class = c("tbl_df", "tbl", "data.frame"), 
           row.names = c(NA, -19L))

This data consists of yield of a crop over time. I am interested in identifying any oultiers in the yield data. This is my approach:

sample_size <- nrow(dat)

Construct a model of yield with time

 mod <- lm(dat$yld.harvest ~ dat$year + I(dat$year^2))

Checking if fitting a quadratic term is ok or not

  mod.update <- step(mod, direction = "backward", trace = FALSE)

Calculate cooks distance

  cooksd <- cooks.distance(mod.update)

Find the influential point

  influential <- as.numeric(names(cooksd)[(cooksd > (4/sample_size))])
  # 2, 17  

Remove the influential point

  dat_screen <- dat[-influential, ]

Plot

 plot(dat$year, dat$yld.harvest, col = "red", pch = 1, xlab = "", 
   ylab = "Yield (kg/ha)", type = "b", main = mun, ylim = yrange)
 points(dat_screen$year, dat_screen$yld.harvest, pch = 19, col = "blue", type = "b")

enter image description here

Visually look at the plots, I do not think the point 2 and 17 should be an outlier. Am I doing it correctly?

$\endgroup$
9
$\begingroup$

The method is correct. The problem with outlier detection is: there is not a general solution and real thresholds. It should be used additionally but not primary to assess your data. You find in different literature different threshold. The inventor of this method suggests your threshold with D > 1 1. Your method (4/sample size) is suggested by other authors 2. Another solution is 4*mean(cooksd) like in boxplots (but I did not found a source for this one). If you use this, you don't detect outliers from your code:

cooks distance with two different thresholds

Also the method is for regression models, which is not very accurate with few values, it should be better to use a larger data sample.

So there is not really perfect solution for detecting outliers. And if you think they are not outliers it should be fine not to remove then.

PS: I found a similar post while doing research: How to read Cook's distance plots? There it is already discussed better than by myself.


  1. Cook, R. D. and Weisberg, S. (1984) Residuals and Influence in Regression.
  2. Bollen, Kenneth A.; Jackman, Robert W. (1990). Fox, John; Long, J. Scott, eds. Regression Diagnostics: An Expository Treatment of Outliers and Influential Cases. Modern Methods of Data Analysis.
$\endgroup$
  • 1
    $\begingroup$ Welcome to the site. I hope we'll see more of these in the future. $\endgroup$ – gung Jun 27 '18 at 14:00
3
$\begingroup$

The method is correct IN THIS CASE but not in general . You need a sufficient model to detect an anomaly ( read underlying signal and latent deterministic structure ).

You approach assumes the kind of model i.e. a model with time and versions of time as exogenous predictors premising no pulses/level shifts/seasonal pulses and attempts to find what had been assumed to be not present. Why is my high degree polynomial regression model suddenly unfit for the data? should ne reviewd and especially @whuber comment Does the p-value in the incremental F-test determine how many trials I expect to get correct? reflecting on "fitting polynomials to data can be a deceptively very poor approach "

The problem or opportunity is to simultaneously identify both i.e. a possible hybrid model using both memory and the waiting to be identified latent deterministic structure.

I used my favorite time series package ( specifically authored to do both with my help) and obtained in this trivial case enter image description here and enter image description here . It concluded as you did regarding the pulse and the level shift.

In summary the noise model can be other than white noise i.e. not (0,0,0) in arima notation. The baseline can have break points in trend not simply curvature which is implied by time squares et al.

A solid reference to detecting amomalies is here http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html covering the impact of error variance changes while incorporating memory as it searches for anomalies.

Standard Regression Diagnostics assume uncorrelated errors while attempting to provide clarity to anomalies. If you have an error process that is white noise then simpler procedures often work BUT when you have time series data this is often not the case thus more powerful/correct/robust approaches are needed.

$\endgroup$
  • $\begingroup$ Thank you. This is quite a detailed response and valuable resource. Regards $\endgroup$ – Crop89 Jun 27 '18 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.