bfastpp {bfast}

Time Series Preprocessing for BFAST-Type Models

Description

Time series preprocessing for subsequent regression modeling. Based on a (seasonal) time series, a data frame with the response, seasonal terms, a trend term, (seasonal) autoregressive terms, and covariates is computed. This can subsequently be employed in regression models.

Usage

bfastpp(
  data,
  order = 3,
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1
)

Arguments

data

A time series of class ts, or another object that can be coerced to such. For seasonal components, a frequency greater than 1 is required.

order

numeric. Order of the harmonic term, defaulting to 3.

lag

numeric. Orders of the autoregressive term, by default omitted.

slag

numeric. Orders of the seasonal autoregressive term, by default omitted.

na.action

function for handling NAs in the data (after all other preprocessing).

stl

character. Prior to all other preprocessing, STL (season-trend decomposition via LOESS smoothing) can be employed for trend-adjustment and/or season-adjustment. The "trend" or "seasonal" component or both from stl are removed from each column in data. By default ("none"), no STL adjustment is used.

decomp

"stlplus" or "stl": use the NA-tolerant decomposition package or the reference package (which can make use of time series with 2-3 observations per year)

sbins

numeric. Controls the number of seasonal dummies. If integer > 1, sets the number of seasonal dummies to use per year. If <= 1, treated as a multiplier to the number of observations per year, i.e. ndummies = nobs/year * sbins.

Details

To facilitate (linear) regression models of time series data, bfastpp facilitates preprocessing and setting up regressor terms. It returns a data.frame containing the first column of the data as the response while further columns (if any) are used as covariates xreg. Additionally, a linear trend, seasonal dummies, harmonic seasonal terms, and (seasonal) autoregressive terms are provided.

Optionally, each column of data can be seasonally adjusted and/or trend-adjusted via STL (season-trend decomposition via LOESS smoothing) prior to preprocessing. The idea would be to capture season and/or trend nonparametrically prior to regression modelling.

Value

If no formula is provided, bfastpp returns a "data.frame" with the following variables (some of which may be matrices).

time

numeric vector of time stamps,

response

response vector (first column of data),

trend

linear time trend (running from 1 to number of observations),

season

factor indicating season period,

harmon

harmonic seasonal terms (of specified order),

lag

autoregressive terms (or orders lag, if any),

slag

seasonal autoregressive terms (or orders slag, if any),

xreg

covariate regressor (all columns of data except the first, if any).

If a formula is given, bfastpp returns a list with components X, y, and t, where X is the design matrix of the model, y is the response vector, and t represents the time of observations. X will only contain variables that occur in the formula. Columns of X have names as decribed above.

Author(s)

Achim Zeileis

References

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi: 10.1016/j.rse.2012.02.022, https://doi.org/10.1016/j.rse.2012.02.022.

See Also

bfastmonitor

Examples

## set up time series
ndvi <- as.ts(zoo::zoo(cbind(a = som$NDVI.a, b = som$NDVI.b), som$Time))
ndvi <- window(ndvi, start = c(2006, 1), end = c(2009, 23))

## parametric season-trend model
d1 <- bfastpp(ndvi, order = 2)
d1lm <- lm(response ~ trend + harmon, data = d1)
summary(d1lm)
# plot visually (except season, as it's a factor)
plot(zoo::read.zoo(d1)[,-3],
  # Avoid clipping plots for pretty output
  ylim = list(c(min(d1[,2]), max(d1[,2])),
            c(min(d1[,3]), max(d1[,3])),
            c(-1, 1), c(-1, 1), c(-1, 1), c(-1, 1),
            c(min(d1[,6]), max(d1[,6]))
       ))

## autoregressive model (after nonparametric season-trend adjustment)
d2 <- bfastpp(ndvi, stl = "both", lag = 1:2)
d2lm <- lm(response ~ lag, data = d2)
summary(d2lm)

## use the lower level lm.fit function
d3 <- bfastpp(ndvi, stl = "both", lag = 1:2)
d3mm <- model.matrix(response ~ lag, d3)
d3lm <- lm.fit(d3mm, d3$response)
d3lm$coefficients