nscosinor: Non-stationary cosinor

Description

Decompose a time series using a non-stationary cosinor for the seasonal pattern.

Usage

nscosinor(data, response, cycles, niters=1000, burnin=500, tau, inits,
                 lambda=1/12, div=50, monthly=TRUE, alpha=0.05)

Arguments

data

a data frame.

response

response variable.

cycles

vector of cycles in units of time, e.g., for a six and twelve month pattern cycles=c(6,12).

niters

total number of MCMC samples (default=1000).

burnin

number of MCMC samples discarded as a burn-in (default=500).

tau

vector of smoothing parameters, tau[1] for trend, tau[2] for 1st seasonal parameter, tau[3] for 2nd seasonal parameter, etc. Larger values of tau allow more change between observations and hence a greater potential flexibility in the trend and season.

inits

vector of initial value(s) for the MCMC chain(s) for std.season. Initial values should be given for each seasonal cycle.

lambda

distance between observations (lambda=1/12 for monthly data, default).

div

divisor at which MCMC sample progress is reported (default=50).

monthly

TRUE for monthly data.

alpha

Statistical significance level used by the confidence intervals.

...

further arguments passed to or from other methods.

Value

Returns an object of class nsCosinor with the following parts:
callthe original call to the nscosinor function.
timethe year and month for monthly data.
trendmean trend and 95% confidence interval.
seasonmean season(s) and 95% confidence interval(s).
oseasonoverall season(s) and 95% confidence interval(s). This will be the same as season if there is only one seasonal cycle.
fittedfitted values, based on trend + season(s).
residualsresiduals based on mean trend and season(s).
nthe length of the series.
chainsMCMC chains (of class mcmc) of variance estimates: standard error for overall noise (std.error), standard error for season(s) (std.season), phase(s) and amplitude(s)
cyclesvector of cycles in units of time.

Details

This model is designed to decompose an equally spaced time series into a trend, season(s) and noise. A seasonal estimate is estimated as $s_t=A_t\cos(\omega_t-P_t)$, where t is time, $A_t$ is the non-stationary amplitude, $P_t$ is the non-stationary phase and $\omega_t$ is the frequency. A non-stationary seasonal pattern is one that changes over time, hence this model gives potentially very flexible seasonal estimates. The frequency of the seasonal estimate(s) are controlled by cycle. The cycles should be specified in units of time. If the data is monthly, then setting lambda=1/12 and cycles=12 will fit an annual seasonal pattern. If the data is daily, then setting lambda= 1/365.25 and cycles=365.25 will fit an annual seasonal pattern. Specifying cycles= c(182.6,365.25) will fit two seasonal patterns, one with a twice-annual cycle, and one with an annual cycle. The estimates are made using a forward and backward sweep of the Kalman filter. Repeated estimates are made using Markov chain Monte Carlo (MCMC). For this reason the model can take a long time to run (we aim to improve this in the next version). To give stable estimates a reasonably long sample should be used (niters), and the possibly poor initial estimates should be discarded (burnin).

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. Barnett, A.G., Dobson, A.J. (2004) Estimating trends and seasonality in coronary heart disease Statistics in Medicine. 23(22) 3505--23.

Examples

Run this code

data(CVD)
# model to fit an annual pattern to the monthly cardiovascular disease data
f = c(12)
inits = c(1)
tau = c(130,10)
res12 = nscosinor(data=CVD, response=adj, cycles=f, niters=5000,
         burnin=1000, tau=tau, inits=inits)
summary(res12)
plot(res12)
plot(res12$chains$amp)
res12 = nscosinor(data=CVD, response=adj, cycles=f, niters=50, burnin=10, tau=tau, inits=inits)

Run the code above in your browser using DataLab