#' @title BFM With Magnitude Threshold
#' @description bfm with automatic magn thresholding and re-monitoring
#'
#' @param x Numeric. Time series vector
#' @param dates Time (date) index of \code{x}
#' @param start Numeric. Start of monitoring period. See \code{\link{bfastmonitor}}
#' @param monend Numeric. Optional: end of the monitoring period. See \code{\link{bfmSpatial}}
#' @param magnThresh Numeric. Threshold for breakpoints. All breakpoints with magnitude < \code{magnThresh} are ignored, all data \code{adjust} years preceding this breakpoint are deleted, and monitoring continues with a revised history period.
#' @param magnPeriod Numeric. Period of time (in years) from which to compute magnitude (as median of residuals)
#' @param adjust Numeric. Delete all observations within this many years preceding a breakpoint, if that breakpoint has a magnitude less than \code{magnThresh}
#' @param verbose Logical. Display error messages (if any)?
#' @param formula See \code{\link{bfastmonitor}}
#' @param order See \code{\link{bfastmonitor}}
#' @param returnBFM Logical. Return a regular \code{bfastmonitor} object?
#' @param ... Arguments to pass to \code{\link{bfastmonitor}}
#'
#' @return Numeric vector with output: break date, magnitude, length of history period (in years), R-squared of fitted model, Adjusted R-squared of fitted model, Intercept of fitted model, cosine coefficient(s) (if applicable), sine coefficient(s) (if applicable), error flag (1 means an error was encountered in the \code{try()} statement). If \code{returnBFM} is \code{FALSE}, a list with the above values as first item and a regular BFM object as second item is returned.
#'
#' @details
#' The magnitude is computed as the median of the residuals within the period \code{magnPeriod} years after the break date. If this magnitude is greater than \code{magnThresh}, all observations preceding the break date by \code{adjust} years are deleted (including the observation at the break date), the start of the monitoring period is reset to the break date, the model is re-fit and monitoring continues until the end of the time series is reached, or a breakpoint with magnitude > \code{magnThresh} is found.
#'
#' The default value for \code{magnPeriod} is 1, assuming that no follow-up land use or other dynamics take place within a 1-year period.
#'
#' The amplitude of the fitted 1st order model can be computed by taking \code{sqrt{harmoncos^2 + harmonsin^2}}
#'
#' @author Ben DeVries
#' @import bfastSpatial
#'
#' @export
#'
#' @examples
#' library(bfmThresh)
#' data(tura)
#' s <- getSceneinfo(names(tura))
#' xy <- c(820692.4, 829564.4)
#' z <- zoo(tura[cellFromXY(tura, xy)][1, ], s$date)
#' z <- window(z, start = "1999-01-01")
#' bfm <- bfmThresh(z, dates = time(z), start = c(2005, 1), plot = TRUE)
#' print(bfm)
#' bfm2 <- bfmThresh(z, dates = time(z), start = c(2005, 1), magnThresh = -500, plot = TRUE)
#' print(bfm2)
bfmThresh <- function(x, dates, start, monend = NULL, magnThresh = 0, magnPeriod = 1, adjust = 1, verbose = FALSE, formula = response ~ harmon, order = 1, returnBFM = FALSE, ...) {
bts <- bfastts(x, dates, type = 'irregular')
# trim time series by monend
if(!is.null(monend))
bts <- window(bts, end = monend)
# magnitude monitor
iterate <- TRUE
err <- NA
# determine length of coefficient vector
# = intercept [+ trend] [+ harmoncos*order] [+ harmonsin*order]
coef_len <- 1 # intercept
modterms <- attr(terms(formula), "term.labels")
if("trend" %in% modterms)
coef_len <- coef_len + 1
if("harmon" %in% modterms)
coef_len <- coef_len + (order * 2) # sin and cos terms
while(iterate) {
bfm <- try(bfastmonitor(bts, start = start, formula = formula, order = order, ...), silent = !verbose)
if(class(bfm) == 'try-error') {
bkp <- magn <- hist <- rsq <- adjrsq <- NA
coefs <- rep(NA, coef_len)
err <- 1
break
} else {
bkp <- bfm$breakpoint
hist <- bfm$history[2] - bfm$history[1]
rsq <- summary(bfm$model)$r.squared
adjrsq <- summary(bfm$model)$adj.r.squared
coefs <- coef(bfm$model)
}
# breakpoint magnitude
if(!is.na(bkp)) {
postpp <- subset(bfm$tspp, time >= bkp & time <= (bkp + magnPeriod))
magn <- median(postpp$response - postpp$prediction)
} else {
magn <- NA
}
if(!is.na(magn) & magn >= magnThresh) {
bts[(time(bts) <= bkp) & (time(bts) > bkp - adjust)] <- NA
start <- c(floor(bkp), 1)
} else {
iterate <- FALSE
}
}
res <- c(bkp, magn, hist, rsq, adjrsq, coefs, err)
names(res) <- c("breakpoint", "magnitude", "history", "r.squared", "adj.r.squared", names(coefs), "error")
if(returnBFM)
res <- list(res, bfm)
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.