library(Rbeast)
# \donttest{
#------------------------------------Example 1----------------------------------------#
# 'googletrend_beach' is the monthly Google search popularity of 'beach' in the US
# from 2004 to 2022, obtained from Google Trends. Sudden changes in the time series
# coincide with known extreme weather events (e.g., 2006 North American Blizzard, 2011
# US hottest summer on record, Record warm January in 2016) or the covid19 outbreak
o <- beast(googletrend_beach)
plot(o)
plot(o, vars=c('t','slpsgn') ) # plot the trend and probability of slope sign only.
# In the slpsgn panel, the upper red portion refers to
# probability of trend slope being positive, the middle
# green to the prob of slope being zero, and the lower
# blue to the probability of slope being negative.
#------------------------------------Example 2----------------------------------------#
# Yellowstone is a half-monthly satellite time series of 774 NDVI(vegetation greeness)
# observations starting from July 1-15,1981(i.e., start=c(1981,7,7)) at a Yellowstone
# site. It has 24 data points per year (i.e., freq=24). Note that the beast function
# hanldes only evenly-spaced regular time series. Irregular data need to be first
# aggegrated at a regular time interval of your choice--the aggregation functionality
# is implemented in beast.irreg() and beast123().
data(Yellowstone)
plot(Yellowstone)
# Yellowstone is not a object of class 'ts' but a pure vector without time attributes.
# Below, no extra argument is supplied, so default values (i.e.,start=1, deltat=1) are
# used and the time is 1:774. 'freq' is missing and so is guessed via auto-correlation.
# Use of auto-correlation to compute the period of a cyclic time series is not always
# reliable, so it is suggested to always supply 'freq' directly, as in Example 2 and
# Example 3.
out = beast(Yellowstone) # the times assumed to be 1:length(Yellowstone) by default
plot(out)
#------------------------------------Example 3----------------------------------------#
# The time info such as start,delta,and freq is explicitly provided. 'start' can be
# given as (1) a vector comprising year, month,& day, (2) a fractional number, or (3)
# a R's Date. The unit of start and deltat does not necessarily refer to time and can
# be arbitrary (e.g., a sequence of data observed at evenly-spaced distaces along a
# transect or a elevation gradient)
o=beast(Yellowstone, start = c(1981,7,7), deltat=1/24, freq=24)
#o=beast(Yellowstone, start = 1981.5137, deltat=1/24, freq=24)
#o=beast(Yellowstone, start = as.Date('1981-7-7'), deltat=1/24, freq=24)
print(o) # o is a R LIST object with many fields
plot(o) # plot many variables
plot(o, vars=c('y','s','t') ) # plot the Y, seasonal, and trend components only
plot(o, vars=c('s','scp','samp','t','tcp','tslp'))# plot some selected variables in
# 'o'. Type "?plot.beast" to see
# more about vars
plot(o, vars=c('s','t'),col=c('red','blue') ) # specify colors of selected subplots
plot(o$time, o$season$Y,type='l') # directly plot output: the fitted season
plot(o$time, o$season$cpOccPr) # directly plot output: season chgpt prob
plot(o$time, o$trend$Y,type='l') # directly plot output: the fitted trend
plot(o$time, o$trend$cpOccPr) # directly plot output: trend chgpt occurrence prob
plot(o$time, o$season$order) # directly plot output: avg harmonic order
plot(o, interactive=TRUE) # manually choose which variables to plot
# }
# \donttest{
#------------------------------------Example 4----------------------------------------#
# Specify other arguments explicitly. Default values are used for missing parameters.
# Note that beast(), beast.irreg(), and beast123() call the same internal C/C++ library,
# so in beast(), the input parameters will be converted to metadata, prior, mcmc, and
# extra parameters as explained for the beast123() function. Or type 'View(beast)' to
# check the parameter assignment in the code.
out = beast(Yellowstone, freq=24,season='dummy', mcmc.samples=5000, tseg.min=20)
plot(out)
out=beast(
Yellowstone, # Yellowstone: a pure numeric vector wo time info
start=1981.51, deltat=1/24, freq = 24,
season = 'harmonic', # periodic compnt exisits,fitted as a harmonic curve
scp.minmax = c(0,3), # min and max numbers of seasonal changpts allowed
sorder.minmax = c(1,5), # min and max harmonic orders allowed
sseg.min = 24, # the min length of segments btw neighboring chnpts
tcp.minmax = c(0,10),# min and max numbers of changpts allowed in the trend
torder.minmax = c(0,1), # min and maxx polynomial orders to fit trend
tseg.min = 24, # the min length of segments btw neighboring trend chnpts
deseasonalize = TRUE, # remove the global seasonality before fitting the beast model
detrend = TRUE, # remove the global trend before fitting the beast model
mcmc.seed = 0, # a seed for mcmc's random nummber generator; use a
# non-zero integer to reproduce results across runs
mcmc.burnin = 500, # number of initial iterations discarded
mcmc.chains = 2, # number of chains
mcmc.thin = 3, # include samples every 3 iterations
mcmc.samples = 6000 # number of samples taken per chain
# total iteration: (500+3*6000)*2
)
plot(out)
plot(out, interactive=TRUE)
# }
#------------------------------------Example 5----------------------------------------#
# Run an interactive GUI to visualize how BEAST is samplinig from the possible model
# spaces in terms of the numbers and timings of seasonal and trend changepoints.
# The GUI inferface allows changing the option parameters interactively. This GUI is
# only available on Win x64 machines, not Mac or Linux.
if (FALSE) {
beast(Yellowstone, freq=24, gui=TRUE)
}
#------------------------------------Example 6----------------------------------------#
# Apply beast to trend-only data. 'Nile' is the ANNUAL river flow of the river
# Nile at Aswan since 1871. It is a 'ts' object; its time attributes (start=1871,
# end=1970,frequency=1) are used to replace the user-supplied start,deltat, and freq,
# if any.
# \donttest{
data(Nile)
plot(Nile)
attributes(Nile) # a ts object with time attributes (i.e., tsp=(start,end,freq)
o = beast(Nile) # start=1871, delta=1, and freq=1 taken from Nile itself
plot(o)
o = beast(Nile, # the same as above. The user-supplied values (i.e., 2011,
start=2011, # 24434) are ignored bcz Nile has its own time attributes.
freq =24434, # Its frequency tag is 1 (i.e., trend-only), so season='none'
season='harmonic' # is used instead of the supplied 'harmonic'
)
# }
#------------------------------------Example 7----------------------------------------#
# NileVec is a pure data vector. The first run below is WRONG bcz NileVec was assumed
# to have a perodic component by default and beast gets a best estimate of freq=6 while
# the true value is freq=1. To fit a trend-only model, season='none' has to be explicitly
# specified, as in the 2nd & 3rd funs.
# \donttest{
NileVec = as.vector(Nile) # NileVec is not a ts obj but a pure numeric data vector
o = beast(NileVec) # WRONG WAY to call: No time attributes available to interpret
# NileVec. By default, beast assumes season='harmonic', start=1,
# & deltat=1. 'freq' is missing and guessed to be 6 (WRONG).
plot(o) # WRONG Results: The result has a suprious seasonal component
o=beast(NileVec,season='none') # Correct WAY to call: Use season='none' for trend-only
# analysis; the default time is the integer indices
# "1:length(NileVec)'.
o=beast(NileVec, # Recommended WAY to call: The true time attributes are
start = 1871, # given explicitly through start and deltat (or freq if
deltat = 1, # there is a cyclic/seasonal cmponent).
season = 'none')
plot(o)
# }
#------------------------------------Example 8----------------------------------------#
# beast can handle missing data. co2 is a monthly time series (i.e.,freq=12) starting
# from Jan 1959. We generate some missing values at random indices
if (FALSE) {
data(co2)
attributes(co2) # A ts object with time attributes (i.e., tsp)
badIdx = sample( 1:length(co2), 50) # Get a set of random indices
co2[badIdx] = NA # Insert some data gaps
out=beast(co2) # co2 is a ts object and its 'tsp' time attributes are used to get the
# true time info. No need to specify 'start','deltat', & freq explicity.
out=beast(co2, # The supplied time/freq values will be ignored bcz
start = c(1959,1,15),# co2 is a ts object; the correct freq = 12 will be
deltat = 1/12, # used.
freq = 365)
print(out)
plot(out)
}
#------------------------------------Example 9----------------------------------------#
# Apply beast to time seris-like sequence data: the unit of sequences is not
# necessarily time.
# \donttest{
data(CNAchrom11) # DNA copy number alterations in Chromesome 11 for cell line GM05296
# The data is orderd by genomic position (not time), and the values
# are the log2-based intensity ratio of copy numbers between the sample
# the reference. A value of zero means no gain or loss in copy number.
o = beast(CNAchrom11,season='none') # season is a misnomer here bcz the data has nothing
# to do with time. Regardless, we fit only a trend.
plot(o)
# }
#------------------------------------Example 10---------------------------------------#
# Apply beast to time seris-like data: the unit of sequences is not necessarily time.
# \donttest{
# Age of Death of Successive Kings of England
# If the data link is deprecated, install the time series data library instead,
# which is available at https://pkg.yangzhuoranyang.com/tsdl/
# install.packages("devtools")
# devtools::install_github("FinYang/tsdl")
# kings = tsdl::tsdl[[293]]
kings = scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)
out = beast(kings,season='none')
plot(out)
# }
#------------------------------------Example 11---------------------------------------#
# Another example from the tsdl data library
# \donttest{
# Number of monthly births in New York from Jan 1946 to Dec 1959
# If the data link becomes invalid, install the time series data package instead
# install.packages("devtools")
# devtools::install_github("FinYang/tsdl")
# kings = tsdl::tsdl[[534]]
births = scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
out = beast(births,start=c(1946,1,15), deltat=1/12 )
plot(out) # the result is wrong bcz the guessed freq via auto-correlation by beast
# is 2 rather than 12, so we recommend always specifying 'freq' explicitly
# for those time series with a periodic component, as shown below.
out = beast(births,start=c(1946,1,15), deltat=1/12, freq=12 )
plot(out)
# }
#------------------------------------Example 12---------------------------------------#
# Daily confirmed COVID-19 new cases and deaths across the globe
if (FALSE) {
data(covid19)
newcases = covid19$newcases
# This ts varies periodically every 7 days. 7 days can't be precisely represented
# in the unit of year bcz some years has 365 days and others has 366. So, here we
# use the date number as the time unit--the num of days lapsed since 1970-01-01.
datenum = as.numeric(covid19$date)
o = beast(newcases, start=min(datenum), deltat=1, freq=7)
o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date.
plot(o)
# Apply BEAST to the square root-transformed time series
o = beast(sqrt(newcases), start=min(datenum), deltat=1, freq=7)
o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date.
plot(o)
}
#------------------------------------Example 13---------------------------------------#
# The old API interface of beast is still made available but NOT recommended. It is
# kept mainly to ensure the working of the sample code on Page 475 in the text
# Ecological Metods by Drs. Southwood and Henderson.
if (FALSE) {
# The interface as shown here will be deprecated and NOT recommended.
beast(Yellowstone, 24) #24 is the freq: number of datapoints per period
# Specify the model or MCMC parameters through opt as in Rbeast v0.2
opt=list() #Create an empty list to append individual model parameters
opt$period=24 #Period of the cyclic component (i.e.,freq in the new version)
opt$minSeasonOrder=2 #Min harmonic order allowed in fitting season component
opt$maxSeasonOrder=8 #Max harmonic order allowed in fititing season component
opt$minTrendOrder=0 #Min polynomial order allowed to fit trend (0 for constant)
opt$maxTrendOrder=1 #Max polynomial order allowed to fit trend (1 for linear term)
opt$minSepDist_Season=20#Min separation time btw neighboring season changepoints
opt$minSepDist_Trend=20 #Min separation time btw neighboring trend changepoints
opt$maxKnotNum_Season=4 #Max number of season changepoints allowed
opt$maxKnotNum_Trend=10 #Max number of trend changepoints allowed
opt$omittedValue=NA #A customized value to indicate bad/missing values in the time
#series, in additon to those NA or NaN values.
# The following parameters used to configure the reverisible-jump MCMC (RJMCC) sampler
opt$chainNumber=2 #Number of parallel MCMC chains
opt$sample=1000 #Number of samples to be collected per chain
opt$thinningFactor=3 #A factor to thin chains
opt$burnin=500 #Number of burn-in samples discarded at the start
opt$maxMoveStepSize=30 #For the move proposal, the max window allowed in jumping from
#the current changepoint
opt$resamplingSeasonOrderProb=0.2 #The probability of selecting a re-sampling proposal
#(e.g., resample seasonal harmonic order)
opt$resamplingTrendOrderProb=0.2 #The probability of selecting a re-sampling proposal
#(e.g., resample trend polynomial order)
opt$seed=65654 #A seed for the random generator: If seed=0,random numbers differ
#for different BEAST runs. Setting seed to a chosen non-zero integer
#will allow reproducing the same result for different BEAST runs.
beast(Yellowstone, opt)
}
Run the code above in your browser using DataLab