Merge pull request #42 from jread-usgs/master

rolling window for MAD(x,w)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: sensorQC
 Type: Package
 Title: sensorQC
-Version: 0.3.4
+Version: 0.4.0
 Date: 2015-10-23
 Author: Jordan S Read, Brad Garner, Brian Pellerin, Luke Loken
 Maintainer: Jordan S Read <jread@usgs.gov>
@@ -15,9 +15,11 @@ Copyright: This software is in the public domain because it contains materials
 Imports:
     yaml,
     dplyr,
+    RcppRoll,
     lazyeval,
     readr
 Suggests:
     testthat
 LazyLoad: yes
 LazyData: yes
+RoxygenNote: 5.0.0

NAMESPACE

@@ -1,4 +1,4 @@
-# Generated by roxygen2 (4.1.1): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 S3method("[",sensor)
 S3method(calc_flags,sensor)
@@ -24,6 +24,7 @@ export(persist)
 export(read)
 export(read.default)
 export(sensor)
+importFrom(RcppRoll,roll_median)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,group_by_)
 importFrom(dplyr,mutate_)

R/custom-functions.R

@@ -14,13 +14,44 @@ MAD.values <- function(vals, b = 1.4826){
   return(MAD.normalized)
 }
 
+
+MAD.roller <- function(vals, window){
+  b = 1.4826
+  warning('MAD.roller function has not been robustly tested w/ NAs')
+  u.i <- is.finite(vals)
+  left.fill <- median(head(vals[u.i], ceiling(window/2)))
+  right.fill <- median(tail(vals[u.i], ceiling(window/2)))
+  medians <- roll_median(vals[u.i], n=window, fill=c(left.fill, 0, right.fill))
+  abs.med.diff <- abs(vals[u.i]-medians)
+  left.fill <- median(head(abs.med.diff, ceiling(window/2)))
+  right.fill <- median(tail(abs.med.diff, ceiling(window/2)))
+  abs.med <- roll_median(abs.med.diff, n=window, fill=c(left.fill, 0, right.fill))
+  MAD <- abs.med*b
+  MAD.normalized = rep(NA,length(vals))
+  MAD.normalized[u.i] <- abs.med.diff/MAD # division by zero
+  MAD.normalized[is.na(MAD.normalized)] <- 0
+  return(MAD.normalized)
+}
+
 MAD.windowed <- function(vals, windows){
+
   stopifnot(length(vals) == length(windows))
-  . <- '_dplyr_var'
-  mad <- group_by_(data.frame(x=vals,w=windows), 'w') %>% mutate_(mad='sensorQC:::MAD.values(x)') %>% .$mad
-  return(mad)
+  if (length(unique(windows)) == 1){
+    w = unique(windows)
+    x = vals
+    return(MAD.roller(x, w))
+  } else {
+    . <- '_dplyr_var'
+    mad <- group_by_(data.frame(x=vals,w=windows), 'w') %>% mutate_(mad='sensorQC:::MAD.values(x)') %>% .$mad
+    return(mad)
+  }
+
+
 }
 #'@title median absolute deviation outlier test
+#'
+#' @description Median Absolute Deviation test
+#'
 #'@name MAD
 #'@aliases MAD
 #'@aliases median.absolute.deviation
@@ -29,6 +60,7 @@ MAD.windowed <- function(vals, windows){
 #'@return a vector of MAD normalized values relative to an undefined rejection criteria (usually 2.5 or 3).
 #'@keywords MAD
 #'@importFrom dplyr group_by_ mutate_ %>%
+#' @importFrom RcppRoll roll_median
 #'@author
 #'Jordan S. Read
 #'@export

R/window.R

@@ -22,13 +22,20 @@ window.sensor<- function(x, ..., type){
 
   if (type=='auto'){
     windowed.data <- auto.chunk.time(x$sensor)
+  } else if (type=='rolling') {
+    windowed.data <- rolling.window(x$sensor, ...)
   } else {
     windowed.data <- manual.chunk.time(x$sensor, type = type)
   }
 
   return(sensor(windowed.data))
 }
 
+rolling.window <- function(data.in, n){
+  data.in$w <- rep(n, nrow(data.in))
+  return(data.in)
+}
+
 #' @importFrom stats window
 auto.chunk.time <- function(data.in){
 

README.Rmd

@@ -96,3 +96,14 @@ clean(data, 'x > 9999', 'persist(x) > 10', 'MAD(x) > 3', replace=NA)
 ```
 
 
+#### flagging data with a moving window
+The `MAD(x,w)` function can use a rolling window by leveraging the `RcppRoll` R package. 
+
+```{r}
+sensor <- read(file, format="wide_burst", date.format="%m/%d/%Y %H:%M")
+sensor = window(sensor, n=300, type='rolling')
+flag(sensor, 'x == 999999', 'persist(x) > 3', 'MAD(x,w) > 3', 'MAD(x) > 3')
+```
+
+
+

README.md

@@ -25,7 +25,7 @@ This package is still very much in development, so the API may change at any tim
 
 High-frequency aquatic sensor QAQC procedures. `sensorQC` imports data, and runs various statistical outlier detection techniques as specified by the user.
 
-### `sensorQC` Functions (as of v0.3.3)
+### `sensorQC` Functions (as of v0.4.0)
 
 | Function | Title                                                  |
 |----------|:-------------------------------------------------------|
@@ -226,3 +226,44 @@ clean(data, 'x > 9999', 'persist(x) > 10', 'MAD(x) > 3', replace=NA)
     ## 6  2
     ## 7  3
     ## 8  4
+
+#### flagging data with a moving window
+
+The `MAD(x,w)` function can use a rolling window by leveraging the `RcppRoll` R package.
+
+``` r
+sensor <- read(file, format="wide_burst", date.format="%m/%d/%Y %H:%M")
+```
+
+    ## number of observations:5100
+
+``` r
+sensor = window(sensor, n=300, type='rolling')
+flag(sensor, 'x == 999999', 'persist(x) > 3', 'MAD(x,w) > 3', 'MAD(x) > 3')
+```
+
+    ## Warning in MAD.roller(x, w): MAD.roller function has not been robustly
+    ## tested w/ NAs
+
+    ## object of class "sensor"
+    ##                  times     x
+    ## 1  2013-11-01 00:00:00 48.86
+    ## 2  2013-11-01 00:00:01 49.04
+    ## 3  2013-11-01 00:00:02 49.50
+    ## 4  2013-11-01 00:00:03 48.91
+    ## 5  2013-11-01 00:00:04 48.90
+    ## 6  2013-11-01 00:00:05 48.96
+    ## 7  2013-11-01 00:00:06 48.48
+    ## 8  2013-11-01 00:00:07 48.97
+    ## 9  2013-11-01 00:00:08 48.97
+    ## 10 2013-11-01 00:00:09 48.99
+    ## 11 2013-11-01 00:00:10 48.35
+    ## 12 2013-11-01 00:00:11 48.51
+    ## 13 2013-11-01 00:00:12 49.25
+    ## 14 2013-11-01 00:00:13 48.82
+    ## 15 2013-11-01 00:00:14 49.22
+    ##   ...
+    ## x == 999999 (15 flags)
+    ## persist(x) > 3 (4 flags)
+    ## MAD(x,w) > 3 (187 flags)
+    ## MAD(x) > 3 (91 flags)