bfastlite: Detect multiple breaks in a time series

In bfast: Breaks for Additive Season and Trend

Detect multiple breaks in a time series

Description

A combination of bfastpp and breakpoints to do light-weight detection of multiple breaks in a time series while also being able to deal with NA values by excluding them via bfastpp.

Usage

bfastlite(
  data,
  formula = response ~ trend + harmon,
  order = 3,
  breaks = "LWZ",
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1,
  h = 0.15,
  level = 0,
  type = "OLS-MOSUM",
  ...
)

bfast0n(
  data,
  formula = response ~ trend + harmon,
  order = 3,
  breaks = "LWZ",
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1,
  h = 0.15,
  level = 0,
  type = "OLS-MOSUM",
  ...
)

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.
formula	a symbolic description for the model in which breakpoints will be estimated.
order	numeric. Order of the harmonic term, defaulting to 3.
breaks	either a positive integer specifying the maximal number of breaks to be calculated, or a string specifying the information criterion to use to automatically determine the optimal number of breaks (see also logLik). "all" means the maximal number allowed by h is used. NULL is treated as the default of the breakpoints function (i.e. BIC).
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.
h	minimal segment size either given as fraction relative to the sample size or as an integer giving the minimal number of observations in each segment.
level	numeric; threshold value for the sctest.efp test for at least one break in a time series, to save processing time for no-change areas. The test is skipped if set to <= 0.
type	character, indicating the type argument to efp.
...	Additional arguments to breakpoints.

Value

An object of class bfastlite, with three elements:

breakpoints	output from breakpoints, containing information about the estimated breakpoints.
data_pp	preprocessed data as output by bfastpp.
sctest	output from sctest.efp, contatining information about the likelihood of the time series having at least one break.

Author(s)

Dainius Masiliunas, Jan Verbesselt

References

\insertRef

dainiusbfastlitebfast

Examples

plot(simts) # stl object containing simulated NDVI time series
datats <- ts(rowSums(simts$time.series))
# sum of all the components (season,abrupt,remainder)
tsp(datats) <- tsp(simts$time.series) # assign correct time series attributes
plot(datats)

# Detect breaks
bp = bfastlite(datats)

# Default method of estimating breakpoints
bp[["breakpoints"]][["breakpoints"]]

# Plot
plot(bp)

# Custom method of estimating number of breaks (request 2 breaks)
strucchangeRcpp::breakpoints(bp[["breakpoints"]], breaks = 2)

# Plot including magnitude based on RMSD for the cos1 component of harmonics
plot(bp, magstat = "RMSD", magcomp = "harmoncos1", breaks = 2)

# Try with a structural change test
bp <- bfastlite(datats, level=0.05)
print(bp)
plot(bp)

# Details of the structural change test with the type RE
bfastlite(datats, level=0.05, type="RE")$sctest


## Run bfastlite() on a raster
f <- system.file("extdata/modisraster.tif", package="bfast")
modisbrick <- terra::rast(f)

# Run on pixel 10
data <- unlist(modisbrick[10])
ndvi <- bfastts(data, dates, type = c("16-day"))
bfl <- bfastlite(ndvi, breaks = "BIC")
# Get max magnitude by RMSD
max(magnitude(bfl[["breakpoints"]])$Mag[,"RMSD"])

# Wrapper function
bflSpatial <- function(pixels)
{
    ts <- bfastts(pixels, dates, type = c("16-day"))
    bfl <- bfastlite(ts, breaks="BIC")
    bp <- bfl[["breakpoints"]]
    # Return number of breakpoints and max breakpoint magnitude
    if (length(bp[["breakpoints"]]) == 1 && is.na(bp[["breakpoints"]]))
        return(c(0, 0))
    
    return(c(length(bp[["breakpoints"]]), max(magnitude(bp)$Mag[,"RMSD"])))
}

# Run function on each raster pixel
rastbfl <- terra::app(modisbrick, bflSpatial)
terra::plot(rastbfl)

bfast documentation built on Oct. 22, 2024, 5:07 p.m.