# beast.irreg: Bayesian time series decomposition for changepoint, trend,... In Rbeast: Bayesian Change-Point Detection and Time Series Decomposition

 beast.irreg R Documentation

## Bayesian time series decomposition for changepoint, trend, and periodicity or seasonality

### Description

A Bayesian model averaging algorithm called BEAST to decompose time series or 1D sequential data into individual components, such as abrupt changes, trends, and periodic/seasonal variations. BEAST is useful for changepoint detection (e.g., breakpoints or structural breaks), nonlinear trend analysis, time series decomposition, and time series segmentation.

### Usage

```  beast.irreg(
y,
time, deltat,
freq         = NA,
season       = c('harmonic','dummy','none'),
scp.minmax   = c(0,10),  sorder.minmax = c(0,5),  sseg.min = NULL,
tcp.minmax   = c(0,10),  torder.minmax = c(0,1),  tseg.min = NULL,
detrend      = FALSE,
deseasonalize= FALSE,
mcmc.seed    = 0,
mcmc.burnin   = 200,
mcmc.chains  = 3,
mcmc.thin    = 5,
mcmc.samples = 8000,
ci             = FALSE,
precValue      = 1.5,
precPriorType  = c('componentwise','uniform','constant','orderwise'),
print.options  = TRUE,
print.progress = TRUE,
gui            = FALSE,
...
)
```

### Arguments

 `y` a vector for an irregular or unordered time series. Missing values such as NA and NaN are allowed. If `y` is regular (i.e., evenly-spaced in time), use the `beast`function instead. If `y` is a matrix or 3D array consisting of multiple regular or irregular time series (e.g., stacked images), use `beast123` instead. If a list of multiple time series is provided for `y`, the multivariate version of the BEAST algorithm will be invoked to decompose the multiple time series and detect common changepoints altogether. This feature is `experimental` only and under further development. Check `ohio` for a working example. `time` a vector of the same length as `y`; the time vector of datapoints of `y`. It can be a numeric vector or a vector of Dates. If your time is formatted as strings or vectors of year, months, and days, use `beast123()` instead where more formats are supported. `deltat` numeric; a user-specified time interval to aggregate `y` based on `time` into a regular time series. The BEAST model is currently formulated for regular data only, so internally, the `beast.irreg` function will aggregate/re-bin irregular data into regular ones; for the aggregation, the `deltat` MUST be provided to specify the desired bin size or time interval. The unit of delta needs to be consistent with `time`. If `time` takes a numeric vector, the unit is arbitrary and irrelevant to beast. If `time` takes a vector of Dates, the unit is assumed to YEAR; for example, if the desired time interval is 1 month (or 1 day), deltat should be 1/12 (or 1/365). `freq` numeric. Needed only for data with a periodic/cyclic component (i.e., season=`'harmonic'` or `'dummy'`) and ignored for trend-only data (i.e., `season='none'`). It specifies the number of samples/datapoints per cycle (e.g, a time series of monthly observations with an annual period has a frequency of 12); it may be a decimal real number (e.g., a time series of bi-hourly observations with a period of 37.5 hrs has a freq of 37.5/2=18.75). The period of the cyclic component in the unit of time is `period=deltat*freq`. `freq` is not a model parameter of BEAST and it has to be specified by the user. If `freq` is missing, BEAST first attempts to guess its value via auto-correlation before fitting the model, but in this case, `freq` is assumed to be an integer. `season` characters (default to 'harmonic'); specify if `y` has a periodic component or not. Three strings are possible. `'none'`: `y` is trend-only; no periodic components are present in the time series. The args for the seasonal component (i.e.,sorder.minmax, scp.minmax and sseg.max) will be ignored. `'harmonic'`: `y` has a periodic/seasonal component. The term 'season' is a misnomer, being used here to broadly refer to any periodic variations present in `y`. The periodicity is NOT a model parameter estimated by beast but a known constant given by the user through `freq`. By default, the periodic component is modeled as a harmonic curve–a combination of sins and cosines. If `'dummy'`, the same as `'harmonic'` except that the periodic/seasonal component is modeled as a non-parametric curve. The arg `sorder.minmax` is irrelevant and is ignored. `sorder.minmax` a vector of 2 integers (>=1); the min and max harmonic orders considered to fit the seasonal component. `sorder.minmax` is used only used if the time series has a seasonal component (i.e., `season='harmonic'`) and ignored for trend-only data or when `season='dummy'`. If the min and max orders are equal (`sorder.minmax=sorder.minmax`), BEAST assumes a constant harmonic order used and won't infer the posterior probability of harmonic orders. `scp.minmax` a vector of 2 integers (>=0); the min and max number of seasonal changepoints (scp) allowed in segmenting the seasonal component. `scp.minmax` is used only if `y` has a seasonal component (i.e., `season='harmonic' or 'dummy' `) and ignored for trend-only data. If the min and max changepoint numbers are equal, BEAST assumes a constant number of scp and won't infer the posterior probability of the number of changepoints, but it still estimates the occurrence probability of the changepoints over time (i.e., the most likely times at which these changepoints occur). If both the min and max numbers are set to 0, no changepoints are allowed; then a global harmonic model is used to fit the seasonal component, but still, the most likely harmonic order will be inferred if sorder.minmax is not equal to sorder.minmax. `sseg.min` an integer; the min segment length allowed between two neighboring season changepoints. That is, when fitting a piecewise harmonic seasonal model, no two changepoints are allowed to occur within a time window of length `sseg.min`. `sseg.min` must be an unitless integer–the number of time intervals/data points so that the time window in the original unit is `sseg.min*deltat`. `sseg.min` defaults to NULL and its value will be given a default value in reference to `freq`.
 `torder.minmax` a vector of 2 integers (>=0); the min and max orders of the polynomials considered to fit the trend component. The 0-th order corresponds to a constant term/a flat line and the 1st order is a line. If `torder.minmax=torder.minmax`, BEAST assumes a constant polynomial order used and won't infer the posterior probability of polynomial orders. `tcp.minmax` a vector of 2 integers; the min and max number of trend changepoints (tcp) allowed in segmenting the trend component. If the min and max changepoint numbers are equal, BEAST assumes a constant number of changepoints and won't infer the posterior probability of the number of changepoints for the trend, but it still estimates the occurrence probability of the changepoints over time (i.e., the most likely times at which these changepoints occur in the trend). If both the min and max numbers are set to 0, no changepoints are allowed; then a global polynomial trend is used to fit the trend component, but still, the most likely polynomial order will be inferred if torder.minmax is not equal to torder.minmax. `tseg.min` an integer; the min segment length allowed between two neighboring trend changepoints. That is, when fitting a piecewise polynomial trend model, no two changepoints are allowed to occur within a time window of length tseg.min. `tseg.min` must be an unitless integer–the number of time intervals/data points so that the time window in the original unit is `tseg.min*deltat`. `tseg.min` defaults to NULL and its value will be given a default value in reference to `freq` if the time series has a cyclic component.
 `detrend` logical; If `TRUE`, a global trend is first fitted and removed from the time series before running BEAST; after BEAST finishes, the global trend is added back to the BEAST result. `deseasonalize` logical; If `TRUE`, a global seasonal model is first fitted and removed from the time series before running BEAST; after BEAST finishes, the global seasonal curve is added back to the BEAST result. `deseasonalize` is ignored if `season='none'` (i.e., trend-only data). `mcmc.seed` integer (>=0); the seed for the random number generator used for Monte Carlo Markov Chain (mcmc). If `mcmc.seed=0`, an arbitrary seed is picked and the fitting results vary across runs. If fixed to the same non-zero integer, the result can be re-produced for different runs. But the results from the same seed may still vary if run on different computers because the random generator library depends on CPU's instruction sets. `mcmc.chains` integer (>0); the number of MCMC chains. `mcmc.thin` integer (>0); a factor to thin chains (e.g., if thinningFactor=5, samples will be taken every 3 iterations) `mcmc.burnin` integer (>0); the number of burn-in samples discarded at the start of each chain `mcmc.samples` integer (>=0); the number of samples collected per MCMC chain. The total number of iterations is `(burnin+samples*thin)*chains`. `ci` boolean; If `TRUE`, credible intervals (i.e., out\$season\$CI or out\$trend\$CI) will be computed for the estimated seasonal and trend components. Computing CI is time-consuming, due to sorting, so set `ci` to FALSE if a symmetric credible interval (i.e., out\$trend\$SD and out\$season\$SD) suffices. `precValue` numeric (>0); the hyperparameter of the precision prior; the default value is 1.5. `precValue` is useful only when `precPriorType`='constant', as further explained below `precPriorType` characters. It takes one of 'constant', 'uniform', 'componentwise' (the default), and 'orderwise'. Below are the differences between them. `'constant'`: the precision parameter used to parameterize the model coefficients is fixed to a constant specified by `precValue`. In other words, `precValue` is a user-defined hyperparameter and the fitting result may be sensitive to the chosen values of `precValue`. `'uniform'`: the precision parameter used to parameterize the model coefficients is a random variable; its initial value is specified by `precValue`. In other words, `precValue` will be inferred by the MCMC, so the fitting result will be insensitive to the choice in `precValue`. `'componentwise'`: multiple precision parameters are used to parameterize the model coefficients for individual components (e.g., one for season and another for trend); their initial values is specified by `precValue`. In other words, `precValue` will be inferred by the MCMC, so the fitting result will be insensitive to the choice in `precValue`. `'orderwise'`: multiple precision parameters are used to parameterize the model coefficients not just for individual components but also for individual orders of each component; their initial values is specified by `precValue`. In other words, `precValue` will be inferred by the MCMC, so the fitting result will be insensitive to the choice in `precValue`. `print.options` boolean. If `TRUE`,the full list of input parameters to BEAST will be printed out prior to the MCMC inference; the naming for this list (e.g., metadata, prior, and mcmc) differs slightly from the input to `beast`, but there is a one-to-one correspondence (e.g., prior\$trendMinSepDist=tseg.min). Internally, beast converts the input parameters to the forms of metadata, prior,and mcmc. Type 'View(beast)' to see the details or check the `beast123` function. `print.progress` boolean;If `TRUE`, a progressbar will be displayed. `gui` boolean. If `TRUE`, BEAST will be run in a GUI demonstration mode, with a GUI window to show an animation of the MCMC sampling in the model space step by step. Note that "`gui=TRUE`" works only for Windows x64 systems not Windows 32 or Linux/Mac systems. `...` additional parameters. There are many more settings for the implementation but not made available in the beast() interface; please use the function `beast123()` instead

### Value

The output is an object of class "beast". It is a list, consisting of the following variables. In the explanations below, we assume the input `y` is a single time series of length `N`:

 `time` a vector of size `1xN`: the times at the `N` sampled locations. By default, it is simply set to `1:N` `data` a vector, matrix, or 3D array; this is a copy of the input `data` if extra\$dumpInputData = TRUE. If extra\$dumpInputData=FALSE, it is set to NULL. If the original input `data` is irregular, the copy here is the regular version aggregated from the original at the time interval specified by metadata\$deltaTime. `marg_lik` numeric; the average of the model marginal likelihood; the larger marg_lik, the better the fitting for a given time series. `R2` numeric; the R-square of the model fitting. `RMSE` numeric; the RMSE of the model fitting. `sig2` numeric; the estimated variance of the model error. `trend` a list object consisting of various outputs related to the estimated trend component: `ncp`: [Number of ChangePoints]. a numeric scalar; the mean number of trend changepoints. Individual models sampled by BEAST has a varying dimension (e.g., number of changepoints or knots), so several alternative statistics (e.g., `ncp_mode`, `ncp_median`, and `ncp_pct90`) are also given to summarize the number of changepoints. For example, if `mcmc\$samples=10`, the numbers of changepoints for the 10 sampled models are assumed to be c(0, 2, 4, 1, 1, 2, 7, 6, 6, 1). The mean ncp is 3.1 (rounded to 3), the median is 2.5 (2), the mode is 1, and the 90th percentile (ncp_pct90) is 6.5. `ncp_mode`: [Number of ChangePoints]. a numeric scalar; the mode for number of changepoints. See the above for explanations. `ncp_median`: [Number of ChangePoints]. a numeric scalar; the median for number of changepoints. See the above for explanations. `ncp_pct90`: [Number of ChangePoints]. a numeric scalar; the 90th percentile for number of changepoints. See the above for explanations. `ncpPr`: [Probability of the Number of ChangePoints]. A vector of length `(tcp.minmax+1)=tcp.max+1`. It gives a probability distribution of having a certain number of trend changepoints over the range of [0,tcp.max]; for example, `ncpPr` is the probability of having no trend changepoint; `ncpPr[i]` is the probability of having (i-1) changepoints: Note that it is `ncpPr[i]` not `ncpPr[i-1]` because ncpPr is used for having zero changepoint. `cpOccPr`: [ChangePoint OCCurence PRobability]. a vector of length N; it gives a probability distribution of having a changepoint in the trend at each point of time. Plotting `cpOccPr` will depict a continious curve of probability-of-being-changepoint. Of particular note, in the curve, a higher peak indicates a higher chance of being a changepoint only at that particular SINGLE point in time and does not necessarily mean a higher chance of observing a changepoint AROUND that time. For example, a window of cpOccPr values `c(0,0,0.5,0,0)` (i.e., the peak prob is 0.5 and the summed prob is 0.5) is less likely to be a changepoint compared to another window `c(0.1,0.2,0.21,0.2,0.1)` (i.e., the peak prob is 0.21 but the summed prob is 0.71). `order`: a vector of length N; the average polynomial order needed to approximate the fitted trend. As an average over many sampled individual piece-wise polynomial trends, `order` is not necessarily an integer. `cp`: [Changepoints] a vector of length `tcp.max=tcp.minmax`; the most possible changepoint locations in the trend component. The locations are obtained by first applying a sum-filtering to the `cpOccPr` curve with a filter window size of tseg.min and then picking up to a total `prior\$MaxKnotNum`/`tcp.max` of the highest peaks in the filtered curve. NaNs are possible if no enough changepoints are identified. `cp` records all the possible changepoints identified and many of them are bound to be false positives. Do not blindly treat all of them as actual changepoints. `cpPr`: [Changepoints PRobability] a vector of length `tcp.max=tcp.minmax`; the probabilities associated with the changepoints `cp`. Filled with NaNs for the remaining elements if `ncp
 `season` a list object consisting of various outputs related to the estimated seasonal/periodic component: `ncp`: [Number of ChangePoints]. a numeric scalar; the mean number of seasonal changepoints. `ncpPr`: [Probability of the Number of ChangePoints]. A vector of length `(scp.minmax+1)=scp.max+1`. It gives a probability distribution of having a certain number of seasonal changepoints over the range of [0,scp.max]; for example, `ncpPr` is the probability of having no seasonal changepoint; `ncpPr[i]` is the probability of having (i-1) changepoints: Note that the index is i rather than (i-1) because ncpPr is used for having zero changepoint. `cpOccPr`: [ChangePoint OCCurence PRobability]. a vector of length N; it gives a probability distribution of having a changepoint in the seasonal component at each point of time. Plotting `cpOccPr` will depict a continious curve of probability-of-being-changepoint over the time. Of particular note, in the curve, a higher value at a peak indicates a higher chance of being a changepoint only at that particular SINGLE point in time, and does not necessarily mean a higher chance of observing a changepoint AROUND that time. For example, a window of cpOccPr values `c(0,0,0.5,0,0)` (i.e., the peak prob is 0.5 and the summed prob is 0.5) is less likely to be a changepoint compared to another window values `c(0.1,0.2,0.3,0.2,0.1)` (i.e., the peak prob is 0.3 but the summed prob is 0.8). `order`: a vector of length N; the average harmonic order needed to approximate the seasonal component. As an average over many sampled individual piece-wise harmonic curves, `order` is not necessarily an integer. `cp`: [Changepoints] a vector of length `scp.max=scp.minmax`; the most possible changepoint locations in the seasonal component. The locations are obtained by first applying a sum-filtering to the `cpOccPr` curve with a filter window size of sseg.min and then picking up to a total `ncp` of the highest peaks in the filtered curve. If `ncp

### References

1. Zhao, K., Wulder, M.A., Hu, T., Bright, R., Wu, Q., Qin, H., Li, Y., Toman, E., Mallick, B., Zhang, X. and Brown, M., 2019. Detecting change-point, trend, and seasonality in satellite time series data to track abrupt changes and nonlinear dynamics: A Bayesian ensemble algorithm. Remote Sensing of Environment, 232, p.111181 (the beast algorithm paper).

2. Zhao, K., Valle, D., Popescu, S., Zhang, X. and Mallick, B., 2013. Hyperspectral remote sensing of plant biochemistry using Bayesian model averaging with variable and band selection. Remote Sensing of Environment, 132, pp.102-119 (the Bayesian MCMC scheme used in beast).

3. Hu, T., Toman, E.M., Chen, G., Shao, G., Zhou, Y., Li, Y., Zhao, K. and Feng, Y., 2021. Mapping fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).

`beast`, `beast123`, `minesweeper`, `tetris`, `geeLandsat`

### Examples

```
library(Rbeast)

######################################################################################
# Note that the BEAST algorithm is currently implemented to handle only regular time
# series. 'beast.irreg' accepts irregular time series but internally it aggregates them
# into regular ones prior to applying the BEAST model. For the aggregation, both the
# "time" and "deltat" args are needed to specify indvidial times of data points and the
# regular time interval desired. If there is  a cyclic componet, freq in the unit of
# time intervals/datapoints should also be given; if not, a possible value is guessed
# via auto-correlation

######################################################################################
# 'ohio' is a data.frame on an irregular Landsat time series of reflectances & ndvi
# (e.g., surface greenness) at an Ohio site. It has multiple columns of alternative date
# formats, such as year, month, day, doy (date of year), rdate (R's date class), and
# time (fractional year)

data(ohio)
str(ohio)

######################################################################################
# Below, 'time' is given as numeric values, which can be of any arbitray unit. Although
# here 1/12 can be interepreted as 1/12 year or 1 month, BEAST itself doesn't care about
# the time unit. So, the unit of 1/12 is irrelevant for BEAST. 'freq' is missing and a
# guess of it is used.

o=beast.irreg(ohio\$ndvi, time=ohio\$time,deltat=1/12)
plot(o)
print(o)

######################################################################################
# Aggregrate the time series at a monthly interval (deltat=1/12) and xplictly provide
# the 'freq' arg

o=beast.irreg(ohio\$ndvi, time=ohio\$time,deltat=1/12, freq=12)

######################################################################################
# Aggregate the time series at a half-monthly time interval, and the freq becomes 24,
# that is, PERIOD (1 year) = deltat (1/24 year ) X  freq (24)

o=beast.irreg(ohio\$ndvi, time=ohio\$time,deltat=1/24, freq=24)

######################################################################################
# 'time' is given as R's dates. The unit is YEAR. 1/12 refers to 1/12 year or 1 month

o=beast.irreg(ohio\$ndvi, time=ohio\$rdate,deltat=1/12)

######################################################################################
# If the input time has formats (e.g., date strings) other than the above two, use the