beast.irreg: Bayesian time series decomposition for changepoint, trend,...

View source: R/beast.irreg.R

beast.irregR 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 beastfunction 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[1]=sorder.minmax[2]), 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[1] is not equal to sorder.minmax[2].

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[1]=torder.minmax[2], 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[1] is not equal to torder.minmax[2].

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.

  1. '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.

  2. '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.

  3. '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.

  4. '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[2]+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[1] 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[1] 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[2]; 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[2]; the probabilities associated with the changepoints cp. Filled with NaNs for the remaining elements if ncp<tcp.max.

  • cpCI: [Changepoints Credible Interval] a matrix of dimension tcp.max x 2; the credibable intervals for the detected changepoints cp.

  • cpAbruptChange: [Abrupt change at Changepoints] a vector of length tcp.max; the jumps in the fitted trend curves at the detected changepoints cp.

  • Y: a vector of length N; the estimated trend component. It is the Bayesian model averaging of all the individual sampled trend.

  • SD: [Standard Deviation] a vector of length N; the estimated standard deviation of the estimated trend component.

  • CI: [Standard Deviation] a matrix of dimension N x 2; the estimated credible interval of the estimated trend. One vector of the matrix is for the upper envelope and another for the lower envelope.

  • slp: [Slope] a vector of length N; the time-varying slope of the fitted trend component .

  • slpSD: [Standar Deviation of Slope] a vector of length N; the SD of the slope for the trend component.

  • slpSgnPosPr: [PRobability of slope having a positive sign] a vector of length N; the probability of the slope being positive (i.e., increasing trend) for the trend component. For example, if slpSgnPosPr=0.80 at a given point in time, it means that 80% of the individual trend models sampled in the MCMC chain has a positive slope at that point.

  • slpSgnZeroPr: [PRobability of slope being zero] a vector of length N; the probability of the slope being zero (i.e., a flat constant line) for the trend component. For example, if slpSgnZeroPr=0.10 at a given point in time, it means that 10% of the individual trend models sampled in the MCMC chain has a zero slope at that point. The probability of slope being negative can be obtained from 1-slpSgnZeroPr-slpSgnPosPr.

  • pos_ncp:

  • neg_ncp:

  • pos_ncpPr:

  • neg_ncpPr:

  • pos_cpOccPr:

  • neg_cpOccPr:

  • pos_cp:

  • neg_cp:

  • pos_cpPr:

  • neg_cpPr:

  • pos_cpAbruptChange:

  • neg_cpAbruptChange:

  • pos_cpCI:

  • neg_cpCI: The above variables have the same outputs as those variables without the prefix 'pos' and 'neg', except that we differentiate the changepoints with a POStive jump in the trend from those changepoints with a NEGative jump. For example, pos_ncp refers to the average number of trend changepoints that jump up (i.e., positively) in the trend.

  • inc_ncp:

  • dec_ncp:

  • inc_ncpPr:

  • dec_ncpPr:

  • inc_cpOccPr:

  • dec_cpOccPr:

  • inc_cp:

  • dec_cp:

  • inc_cpPr:

  • dec_cpPr:

  • inc_cpAbruptChange:

  • dec_cpAbruptChange:

  • inc_cpCI:

  • dec_cpCI: The above variables have the same outputs as those variables without the prefix 'inc' and 'dec', except that we differentiate the changepoints at which the trend slope increases from those changepoints at which the trend slope decreases. For example, if the trend slopes before and after a chngpt is 0.4 and 2.5, then the changepoint is counted toward inc_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[2]+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[1] 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[1] 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[2]; 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<scp.max, the remaining of the vector is filled with NaNs.

  • cpPr: [Changepoints PRobability] a vector of length scp.max; the probabilities associated with the changepoints cp. Filled with NaNs for the remaining elements if ncp<scp.max.

  • cpCI: [Changepoints Credible Interval] a matrix of dimension scp.max x 2; the credibable intervals for the detected changepoints cp.

  • cpAbruptChange: [Abrupt change at Changepoints] a vector of length scp.max; the jumps in the fitted trend curves at the detected changepoints cp.

  • Y: a vector of length N; the estimated trend component. It is the Bayesian model averaging of all the individual sampled trend.

  • SD: [Standard Deviation] a vector of length N; the estimated standard deviation of the estimated trend component.

  • CI: [Standard Deviation] a matrix of dimension N x 2; the estimated credible interval of the estimated trend. One vector of the matrix is for the upper envelope and another for the lower envelope.

  • amp: [AMPlitude] a vector of length N; the time-varying amplitude of the estimated seasonality.

  • ampSD: [Standar Deviation of AMPlitude] a vector of length N; , the SD of the amplitude of the seasonality.

  • pos_ncp:

  • neg_ncp:

  • pos_ncpPr:

  • neg_ncpPr:

  • pos_cpOccPr:

  • neg_cpOccPr:

  • pos_cp:

  • neg_cp:

  • pos_cpPr:

  • neg_cpPr:

  • pos_cpAbruptChange:

  • neg_cpAbruptChange:

  • pos_cpCI:

  • neg_cpCI: The above variables have the same outputs as those variables without the prefix 'pos' and 'neg', except that we differentiate the changepoints with a POStive jump in the trend from those changepoints with a NEGative jump. For example, pos_ncp refers to the average number of trend changepoints that jump up (i.e., positively) in the trend.

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).

See Also

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 
# beast123() function instead.




Rbeast documentation built on May 18, 2022, 9:05 a.m.

Related to beast.irreg in Rbeast...