Package 'season'

Description

Generates random linear surrogate data of a time series with the same second-order properties.

Usage

aaft(data, nsur)

Arguments

Details

The AAFT uses phase-scrambling to create a surrogate of the time series that has a similar spectrum (and hence similar second-order statistics). The AAFT is useful for testing for non-linearity in a time series, and is used by nonlintest.

Value

a matrix of the nsur surrogates.

Author(s)

Adrian Barnett [email protected]

References

Kugiumtzis D (2000) Surrogate data test for nonlinearity including monotonic transformations, Phys. Rev. E, vol 62 no. 1, 2000. doi:10.1103/PhysRevE.62.R25

Examples

surr <- aaft(CVD$cvd, nsur = 1)
plot(CVD$cvd, type = "l")
lines(surr[ ,1], col = "red")

Description

The data are: a) the monthly frequencies of birthdays and an expected number based on monthly birth statistics for 1975 to 1991. b) all 617 players' initials and birthdays (excluding non-Australian born players).

Usage

AFL

Format

A list with the following 5 variables:

• month: integer month (1 to 12)

• players: number of players born in each month (12 observations)

• expected: expected number of players born in each month (12 observations)

• initials: player initials (617 observations)

• dob: date of birth in date format (617 observations; year-month-day format)

Source

Dates of birth from Wikipedia.

Examples

barplot(AFL$players, names.arg=month.abb)

Plot the fitted sinusoid from a cosinor() model

Description

Returns a ggplot of the fitted sinusoid. The y-axis is on the probability scale for logit and cloglog link functions; the x-axis switches between month, year, and hour depending on the type of the model.

Usage

## S3 method for class 'Cosinor'
autoplot(object, ...)

Arguments

Value

a ggplot object.

Author(s)

Nicholas Tierney

See Also

cosinor(), summary.Cosinor(), seasrescheck()

Examples

res <- cosinor(
  cvd ~ 1,
  date = "month",
  data = CVD,
  type = "monthly",
  family = poisson(),
  offsetmonth = TRUE
)
autoplot(res)
autoplot(res) + ggplot2::theme_minimal()

Plot the monthly estimates from monthglm()

Description

Returns a ggplot of the per-month coefficient estimates from a monthglm() fit, with a confidence interval per month and (for Poisson/binomial families) a horizontal reference line at the rate / odds ratio of 1.

Usage

## S3 method for class 'monthglm'
autoplot(object, alpha = 0.05, ...)

Arguments

Value

a ggplot object.

Author(s)

Nicholas Tierney

See Also

monthglm()

Examples

mmodel <- monthglm(
  formula = cvd ~ 1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop / 100000),
  offsetmonth = TRUE,
  refmonth = 6
)
autoplot(mmodel)
autoplot(mmodel) + ggplot2::labs(x = "Month", y = "Rate ratio")

Plot the monthly mean estimates from monthmean()

Description

Returns a ggplot of the per-month means computed by monthmean().

Usage

## S3 method for class 'Monthmean'
autoplot(object, ...)

Arguments

Value

a ggplot object.

Author(s)

Nicholas Tierney

See Also

monthmean()

Examples

mmean <- monthmean(
  data = CVD,
  resp = "cvd",
  offsetpop = expression(pop / 100000),
  adjmonth = "average"
)
autoplot(mmean)
autoplot(mmean) + ggplot2::theme_minimal()

Description

Returns a ggplot contour plot of the third-order moment region that exceeds the bootstrap test limits computed by nonlintest(). Returns NULL invisibly with an informational message if no points exceed the limits.

Usage

## S3 method for class 'nonlintest'
autoplot(object, ...)

Arguments

Value

a ggplot contour plot, or NULL invisibly when no points exceed the test limits.

Author(s)

Nicholas Tierney

See Also

nonlintest()

Examples

## Not run: 
test_res <- nonlintest(data = CVD$cvd, n.lag = 4, n.boot = 1000)
autoplot(test_res)

## End(Not run)

Plot the trend and seasonal estimates from nscosinor()

Description

Returns a ggplot of the trend and seasonal components estimated by nscosinor(), faceted by component and with a ribbon for the confidence interval.

Usage

## S3 method for class 'nsCosinor'
autoplot(object, ...)

Arguments

Value

a ggplot object faceted by trend / season(s).

Author(s)

Nicholas Tierney

See Also

nscosinor()

Examples

## Not run: 
res <- nscosinor(
  data = CVD, response = "adj", cycles = 12,
  niters = 200, burnin = 100, tau = c(10, 50)
)
autoplot(res)
autoplot(res) + ggplot2::theme_minimal()

## End(Not run)

Plot the periodogram from peri()

Description

Produce a ggplot of the periodogram in both radians and cycles. The returned ggplot can be extended with + (e.g. + ggplot2::theme_minimal()).

Usage

## S3 method for class 'peri'
autoplot(object, ...)

Arguments

Value

a ggplot object.

Author(s)

Nicholas Tierney

See Also

peri()

Examples

p <- peri(CVD$cvd)
autoplot(p)
autoplot(p) + ggplot2::theme_minimal()

Plot the third-order moment from third()

Description

Produce a ggplot contour of the third-order moment over its non-redundant region.

Usage

## S3 method for class 'third'
autoplot(object, ...)

Arguments

Value

a ggplot contour plot.

Author(s)

Nicholas Tierney

See Also

third()

Examples

t <- third(CVD$cvd, n.lag = 4)
autoplot(t)

Description

Fits a time-stratified case–crossover to regularly spaced time series data. The function is not suitable for irregularly spaced individual data. The function only uses a time-stratified design, and other designs such as the symmetric bi-directional design, are not available.

Usage

casecross(
  formula,
  data,
  exclusion = 2,
  stratalength = 28,
  matchdow = FALSE,
  usefinalwindow = TRUE,
  matchconf = NULL,
  confrange = 0,
  stratamonth = FALSE
)

Arguments

Details

The case–crossover method compares "case" days when events occurred (e.g., deaths) with control days to look for differences in exposure that might explain differences in the number of cases. Control days are selected to be nearby to case days, which means that only recent changes in the independent variable(s) are compared. By only comparing recent values, any long-term or seasonal variation in the dependent and independent variable(s) can be eliminated. This elimination depends on the definition of nearby and on the seasonal and long-term patterns in the independent variable(s).

Control and case days are only compared if they are in the same stratum. The stratum is controlled by stratalength, the default value is 28 days, so that cases and controls are compared in four week sections. Smaller stratum lengths provide a closer control for season, but reduce the available number of controls. Control days that are close to the case day may have similar levels of the independent variable(s). To reduce this correlation it is possible to place an exclusion around the cases. The default is 2, which means that the smallest gap between a case and control will be 3 days.

To remove any confounding by day of the week it is possible to additionally match by day of the week (matchdow), although this usually reduces the number of available controls. This matching is in addition to the strata matching.

It is possible to additionally match case and control days by an important confounder (matchconf) in order to remove its effect. Control days are matched to case days if they are: i) in the same strata, ii) have the same day of the week if matchdow=TRUE, iii) have a value of matchconf that is within plus/minus confrange of the value of matchconf on the case day.

If the range is set too narrow then the number of available controls will become too small, which in turn means the number of case days with at least one control day is compromised.

The method uses conditional logistic regression (see survival::coxph()), so the parameter estimates are odds ratios.

The code assumes that the data frame contains a date variable (in Date() format) called "date".

Value

a list with the following elements:

• call: the original call to the casecross function.

• cox_model: conditional logistic regression model of class coxph.

• n_cases: total number of cases.

• n_case_days: number of case days with at least one control day.

• n_control_days: average number of control days per case day.

Author(s)

Adrian Barnett [email protected]

References

Janes, H., Sheppard, L., Lumley, T. (2005) Case-crossover analyses of air pollution exposure data: Referent selection strategies and their implications for bias. Epidemiology 16(6), 717–726. doi:10.1097/01.ede.0000181315.18836.9d.

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

See Also

summary.casecross, coxph

Examples

# cardiovascular disease data
# subset for example
CVDdaily <- subset(CVDdaily, date <= as.Date('1987-12-31'))
# Effect of ozone on CVD death
model1 <- casecross(
  cvd ~ o3mean + tmpd + Mon + Tue + Wed + Thu + Fri + Sat,
  data = CVDdaily
)
summary(model1)
# match on day of the week
model2 <- casecross(cvd ~ o3mean + tmpd, matchdow = TRUE, data = CVDdaily)
summary(model2)
# match on temperature to within a degree
model3 <- casecross(
  cvd ~ o3mean + Mon + Tue + Wed + Thu + Fri + Sat,
  data = CVDdaily,
  matchconf = "tmpd",
  confrange = 1
)
summary(model3)

Description

Calculates the mean and confidence interval for the phase based on a chain of MCMC samples.

Usage

ciPhase(theta, alpha = 0.05)

Arguments

Details

The estimates of the phase are rotated to have a centre of $\pi$, the point on the circumference of a unit radius circle that is furthest from zero. The mean and confidence interval are calculated on the rotated values, then the estimates are rotated back.

Value

a list with the following elements:

• mean: the estimated mean phase.

• lower: the estimated lower limit of the confidence interval.

• upper: the estimated upper limit of the confidence interval.

Author(s)

Adrian Barnett [email protected]

References

Fisher, N. (1993) Statistical Analysis of Circular Data. Cambridge University Press. Page 36.

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

Examples

# 2000 normal samples, centred on zero
theta <- rnorm(n = 2000, mean = 0, sd = pi / 50)
hist(theta, breaks = seq(-pi / 8, pi / 8, pi / 30))
ciPhase(theta)

Description

Fits a cosinor model as part of a generalized linear model.

Usage

cosinor(
  formula,
  date,
  data,
  family = stats::gaussian(),
  alpha = 0.05,
  cycles = 1,
  rescheck = FALSE,
  type = c("daily", "weekly", "monthly", "hourly"),
  offsetmonth = FALSE,
  offsetpop = NULL,
  text = TRUE
)

Arguments

Details

The cosinor model captures a seasonal pattern using a sinusoid. It is therefore suitable for relatively simple seasonal patterns that are symmetric and stationary. The default is to fit an annual seasonal pattern (cycle=1), but other higher frequencies are possible (e.g., twice per year: cycle=2). The model is fitted using a sine and cosine term that together describe the sinusoid. These parameters are added to a generalized linear model, so the model can be fitted to a range of dependent data (e.g., Normal, Poisson, Binomial). Unlike the nscosinor() model, the cosinor model can be applied to unequally spaced data.

Value

Returns an object of class "Cosinor" with the following parts:

• call: the original call to the cosinor function.

• glm: an object of class glm (see glm()).

• fitted: fitted values for intercept and cosinor only (ignoring other independent variables).

• fitted.plus: standard fitted values, including all other independent variables.

• residuals: residuals.

• date: name of the date variable (in Date format when type = daily).

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

See Also

summary.Cosinor(), plot.Cosinor()

Examples

## cardiovascular disease data (offset based on number of days in...
## ...the month scaled to an average month length)
res <- cosinor(
  cvd ~ 1,
  date = 'month',
  data = CVD,
  type = 'monthly',
  family = poisson(),
  offsetmonth = TRUE
)
summary(res)
seasrescheck(res$residuals) # check the residuals
## stillbirth data
res <- cosinor(
  stillborn ~ 1,
  date = 'dob',
  data = stillbirth,
  family = binomial(link = 'cloglog')
)
summary(res)
plot(res)
## hourly indoor temperature data
res <- cosinor(
  bedroom ~ 1,
  date = 'datetime',
  type = 'hourly',
  data = indoor
)

summary(res)
# to get the p-values for the sine and cosine estimates
summary(res$glm)

Description

Creates an adjacency matrix in a form suitable for use with nimble

Usage

createAdj(matrix, suffix = NULL)

Arguments

Details

Adjacency matrices are used by conditional autoregressive (CAR) models to smooth estimates according to some neighbourhood map. The basic idea is that neighbouring areas have more in common than non-neighbouring areas and so will be positively correlated.

As well as correlations in space it is possible to use CAR models to model similarities in time.

In this case the matrix represents those time points that we wish to assume to be correlated.

Value

A list of:

• num: the total number of neighbours

• adj: the index number of the adjacent neighbours

• weights: weights

Author(s)

Adrian Barnett [email protected]

Examples

# Nearest neighbour matrix for 5 time points
x <- c(NA,1,NA,NA,NA)
V <- toeplitz(x)
V
createAdj(V)

Description

Monthly number of deaths from cardiovascular disease in people aged 75 and over in Los Angeles for the years 1987 to 2000.

Usage

CVD

Format

A data frame with 168 observations on the following 8 variables:

• year: year of death

• month: month of death

• yrmon: a combination of year and month: $year+(month-1)/12$

• cvd: monthly number of CVD deaths

• tmpd: mean monthly temperature (degrees Fahrenheit)

• pop: Los Angeles population aged 75+ in the year 2000 (this value is constant as only one year was available, but in general the population will change over time)

• ndaysmonth: number of days in each month (used as an offset)

• adj: adjusted number of CVD deaths per month using a standardised month length. Monthly number of CVD deaths multiplied by (365.25/12)/ndaysmonth. So the standard month length is 30.4 days.

Source

From the NMMAPS study.

References

Samet JM, Dominici F, Zeger SL, Schwartz J, Dockery DW (2000). The National Morbidity, Mortality, and Air Pollution Study, Part I: Methods and Methodologic Issues. Research Report 94, Health Effects Institute, Cambridge MA.

Examples

plot(
  CVD$yrmon,
  CVD$cvd,
  type = 'o',
  xlab = 'Date',
  ylab = 'Number of CVD deaths per month'
  )

Description

Daily number of deaths from cardiovascular disease in people aged 75 and over in Los Angeles for the years 1987 to 2000.

Usage

CVDdaily

Format

A data frame with 5114 observations on the following 16 variables:

• date: date of death in date format (year-month-day)

• cvd: daily number of CVD deaths

• dow: day of the week (character)

• tmpd: daily mean temperature (degrees Fahrenheit)

• o3mean: daily mean ozone (parts per billion)

• o3tmean: daily trimmed mean ozone (parts per billion)

• Mon: indicator variable for Monday

• Tue: indicator variable for Tuesday

• Wed: indicator variable for Wednesday

• Thu: indicator variable for Thursday

• Fri: indicator variable for Friday

• Sat: indicator variable for Saturday

• month: month (integer from 1 to 12)

• winter: indicator variable for winter

• spring: indicator variable for spring

• summer: indicator variable for summer

• autumn: indicator variable for autumn

Source

From the NMMAPS study.

References

Samet JM, Dominici F, Zeger SL, Schwartz J, Dockery DW (2000). The National Morbidity, Mortality, and Air Pollution Study, Part I: Methods and Methodologic Issues. Research Report 94, Health Effects Institute, Cambridge MA.

Examples

plot(
  CVDdaily$date,
  CVDdaily$cvd,
  type = 'p',
  xlab = 'Date',
  ylab = 'Number of CVD deaths')

Description

Exercise data in longitudinal format from a physical activity intervention study in Logan, Queensland. Some subjects were lost to follow-up, so all three visits are not available for all subjects.

Usage

exercise

Format

A data frame with 1302 observations on the following 7 variables:

• id: subject number

• visit: visit number (1, 2 or 3)

• date: date of interview (year-month-day)

• year: year of interview

• month: month of interview

• bmi: body mass index at visit 1 (kg/m$^2$)

• walking: walking time per week (in minutes) at each visit

Source

From Prof Elizabeth Eakin and colleagues, The University of Queensland, Brisbane.

References

Eakin E, et al (2009) Telephone counselling for physical activity and diet in type 2 diabetes and hypertension, Am J of Prev Med, vol 36, pages 142–9

Examples

boxplot(exercise$walking ~ exercise$month)

Description

The data are indoor temperatures (in degrees C) for a bedroom and living room in a house in Brisbane, Australia for the dates 10 July 2013 to 3 October 2013. Temperatures were recorded using data loggers which recorded every hour to the nearest 0.5 degrees.

Format

A data.frame with the following 3 variables:

• datetime: date and time in POSIXlt format

• living: the living room temperature

• bedroom: the bedroom temperature

Source

Adrian G Barnett.

Examples

res <- cosinor(
 bedroom ~ 1,
 date = 'datetime',
 type = 'hourly',
 data = indoor
)
summary(res)

Description

Returns the day and month (for "daily") or fraction of the month (for "monthly") given a fraction of the year. Assumes a year length of 365.25 days for "daily". When using "monthly" the 1st of January is 1, the 1st of December is 12, and the 31st of December is 12.9. For "hourly" it returns the fraction of the 24-hour clock starting from zero (midnight).

Usage

invyrfraction(
  frac,
  type = c("daily", "monthly", "hourly", "weekly"),
  text = TRUE
)

Arguments

Value

the date (day and month for "daily"), fractional month (for "monthly"), or fraction of the 24-hour clock (for "hourly").

Author(s)

Adrian Barnett [email protected]

Examples

invyrfraction(c(0, 0.5, 0.99), type = "hourly")
invyrfraction(c(0, 0.5, 0.99), type = "daily")
invyrfraction(c(0, 0.5, 0.99), type = "weekly")
invyrfraction(c(0, 0.5, 0.99), type = "monthly")

Description

Fit a generalized linear model with a categorical variable of month.

Usage

monthglm(
  formula,
  data,
  family = stats::gaussian(),
  refmonth = 1,
  monthvar = "month",
  offsetmonth = FALSE,
  offsetpop = NULL
)

Arguments

Details

Month is fitted as a categorical variable as part of a generalized linear model. Other independent variables can be added to the right-hand side of formula.

This model is useful for examining non-sinusoidal seasonal patterns. For sinusoidal seasonal patterns see cosinor().

The data frame should contain the integer months and the year as a 4 digit number. These are used to calculate the number of days in each month accounting for leap years.

Value

a list with the following elements:

• call: the original call to the monthglm function.

• fit: GLM model.

• fitted: fitted values.

• residuals: residuals.

• out: details on the monthly estimates.

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

See Also

summary.monthglm, plot.monthglm

Examples

model <- monthglm(
  formula = cvd~1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop/100000),
  offsetmonth = TRUE
  )
summary(model)

Description

For time series recorded at monthly intervals it is often useful to examine (and plot) the average in each month. When using count data we should adjust the mean to account for the unequal number of days in the month (e.g., 31 in January and 28 or 29 in February).

Usage

monthmean(
  data,
  resp,
  offsetpop = NULL,
  adjmonth = c("none", "thirty", "average")
)

Arguments

Details

This function assumes that the data set (data) contains variables for the year and month called year and month, respectively.

Value

Returns an object of class "Monthmean" with the following parts:

• mean: a vector of length 12 with the monthly means.

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

See Also

plot.Monthmean()

Examples

# cardiovascular disease data
mmean <- monthmean(
  data=CVD,
  resp='cvd',
  offsetpop = expression(pop/100000),
  adjmonth = 'average'
  )
mmean
plot(mmean)

mmean <- monthmean(
  data=CVD,
  resp='cvd',
  offsetpop = expression(pop/100000),
  adjmonth = 'thirty'
  )

Description

A bootstrap test of non-linearity in a time series using the third-order moment.

Usage

nonlintest(data, n.lag, n.boot, alpha = 0.05)

Arguments

Details

The test uses aaft to create linear surrogates with the same second-order properties, but no (third-order) non-linearity. The third-order moments (third) of these linear surrogates and the actual series are then compared from lags 0 up to n.lag (excluding the skew at the co-ordinates (0,0)). The bootstrap test works on the overall area outside the limits, and gives an indication of the overall non-linearity. The plot using region shows those co-ordinates of the third order moment that exceed the null hypothesis limits, and can be a useful clue for guessing the type of non-linearity. For example, a large value at the co-ordinates (0,1) might be caused by a bi-linear series $X_t=\alpha X_{t-1}\varepsilon_{t-1} +\varepsilon_t$.

Value

Returns an object of class "nonlintest" with the following parts:

• region: the region of the third order moment where the test exceeds the limits (up to n.lag).

• n.lag: the maximum lag tested using the third-order moment.

• stats: a list of statistics for the area outside the test limits:

  • outside: the total area outside of limits (summed over the whole third-order moment).

  • stan: the total area outside the limits divided by its standard deviation to give a standardised estimate.

  • median: the median area outside the test limits.

  • upper: the (1-alpha)th percentile of the area outside the limits.

  • pvalue: bootstrap p-value of the area outside the limits to test if the series is linear.

  • test: reject the null hypothesis that the series is linear (TRUE/FALSE).

Author(s)

Adrian Barnett [email protected]

References

Barnett AG & Wolff RC (2005) A Time-Domain Test for Some Types of Nonlinearity, IEEE Transactions on Signal Processing, vol 53, pages 26–33 doi: 10.1109/TSP.2004.838942.

See Also

print.nonlintest() plot.nonlintest()

Examples

## Not run: 
test_res <- nonlintest(data = CVD$cvd, n.lag = 4, n.boot = 1000)

## End(Not run)

Description

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

Usage

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

Arguments

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. To give stable estimates a reasonably long sample should be used (niters), and the possibly poor initial estimates should be discarded (burnin).

Value

Returns an object of class "nsCosinor" with the following parts:

• call: the original call to the nscosinor function.

• time: the year and month for monthly data.

• trend: mean trend and 95\

• season: mean season(s) and 95\

• oseason: overall season(s) and 95\ be the same as season if there is only one seasonal cycle.

• fitted: fitted values and 95\ season(s).

• residuals: residuals based on mean trend and season(s).

• n: the length of the series.

• chains: MCMC 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).

• cycles: vector of cycles in units of time.

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

Barnett, A.G., Dobson, A.J. (2004) Estimating trends and seasonality in coronary heart disease Statistics in Medicine. 23(22) 3505–23.

See Also

plot.nsCosinor() summary.nsCosinor()

Examples

# model to fit an annual pattern to the monthly cardiovascular disease data
f <- c(12)
tau <- c(10,50)
## Not run: 
  res12 <- nscosinor(
    data = CVD,
    response = 'adj',
    cycles = f,
    niters = 200,
    burnin = 50,
    tau = tau
    )
summary(res12)
plot(res12)

## End(Not run)

Description

Estimated periodogram using the fast Fourier transform (fft).

Usage

peri(data, adjmean = TRUE, plot = lifecycle::deprecated())

Arguments

Value

an object of class "peri" (a list) with the following elements:

• peri: periodogram, I($\omega$).

• f: frequencies in radians, $\omega$.

• c: frequencies in cycles of time, $2\pi/\omega$.

• amp: amplitude periodogram.

• phase: phase periodogram.

Pass the result to autoplot() to draw the plot.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.peri()

Examples

p <- peri(CVD$cvd)
autoplot(p)

Description

Calculate the phase given the estimated sine and cosine values from a cosinor model. Returns the phase in radians, in the range $[0,2\pi)$. The phase is the peak in the sinusoid. Applies atan2() over a branching workflow for each coordinate. See https://en.wikipedia.org/wiki/Atan2 for more information.

Usage

phasecalc(cosine, sine)

Arguments

Value

the estimated phase in radians.

Author(s)

Adrian Barnett [email protected]

References

Fisher, N.I. (1993) Statistical Analysis of Circular Data. Cambridge University Press, Cambridge. Page 31.

Examples

# pi/2
phasecalc(cosine=0, sine=1)

Description

Circular plot of a monthly variable. This circular plot can be useful for estimates of an annual seasonal pattern. Darker shades of grey correspond to larger numbers. Pie chart where each segment is one month and the fill is shaded to indicate the value (darker = larger). The first value is assumed to be January. This circular plot can be useful for estimates of an annual seasonal pattern. Darker shades of grey correspond to larger numbers.

Usage

plot_circle(months, dp = 1)

Arguments

Value

a ggplot object.

Author(s)

Nicholas Tierney

See Also

plotCircle() (deprecated base R version)

Examples

plot_circle(months = seq(1, 12, 1), dp = 0)
plot_circle(months = c(1, 3, 5, 7, 9, 11, 12, 10, 8, 6, 4, 2)) +
  ggplot2::labs(title = "Example monthly pattern")

Description

Plots results by month. Assumes the data frame contains variables called year and month. Faceted or coloured line plot of a monthly response over time. Assumes the data frame contains year and month columns.

Usage

plot_month(data, resp, panels = 12)

Arguments

Value

a ggplot object.

See Also

plotMonth() (deprecated wrapper)

Examples

plot_month(data = CVD, resp = "cvd", panels = 12)
plot_month(data = CVD, resp = "cvd", panels = 1) +
  ggplot2::theme_minimal()

Description

Soft-deprecated in favour of autoplot.Cosinor(), which returns a ggplot object you can extend with +. The base R plot below still works.

Usage

## S3 method for class 'Cosinor'
plot(x, ...)

Arguments

Details

The code produces the fitted sinusoid based on the intercept and sinusoid. The y-axis is on the scale of probability if the link function is "logit" or "cloglog". If the analysis was based on monthly data then month is shown on the x-axis. If the analysis was based on daily data then time is shown on the x-axis.

Value

connected line plot of fitted sinusoid object produced by cosinor.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.Cosinor(), cosinor(), summary.Cosinor(), seasrescheck()

Examples

## cardiovascular disease data (offset based on number of days in...
## ...the month scaled to an average month length)
res <- cosinor(
  cvd ~ 1,
  date = 'month',
  data = CVD,
  type = 'monthly',
  family = poisson(),
  offsetmonth = TRUE
  )
# Recommended:
autoplot(res)
# Still works, but deprecated:
plot(res)

Description

Soft-deprecated in favour of autoplot.monthglm(), which returns a ggplot object. The base R plot() will still work, but we recommend using autoplot.monthglm().

Usage

## S3 method for class 'monthglm'
plot(x, alpha = 0.05, ylim = NULL, xlab = "", ylab = "", ...)

Arguments

Value

Plot of the estimated from a generalized linear model with a categorical variable of month.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.monthglm(), monthglm()

Examples

mmodel <- monthglm(
  formula = cvd ~ 1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop / 100000),
  offsetmonth = TRUE,
  refmonth = 6
)
# Recommended:
autoplot(mmodel)
# Still works, but deprecated:
plot(mmodel)

Description

Soft-deprecated in favour of autoplot.Monthmean(), which returns a ggplot object you can extend with +. The base R plot below still works.

Usage

## S3 method for class 'Monthmean'
plot(x, ...)

Arguments

Value

Connected dot plot of estimated monthly means.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.Monthmean(), monthmean()

Examples

mmean <- monthmean(
  data = CVD,
  resp = 'cvd',
  offsetpop = expression(pop / 100000),
  adjmonth = 'average'
)
mmean
# Recommended:
autoplot(mmean)
# Still works, but deprecated:
plot(mmean)

Description

Soft-deprecated in favour of autoplot.nonlintest(), which returns a ggplot object you can extend with +. The base R-style plot below still works.

Usage

## S3 method for class 'nonlintest'
plot(x, plot = lifecycle::deprecated(), ...)

Arguments

Value

ggplot Contour plot of the region of the third-order moment outside the test limits generated by nonlintest().

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.nonlintest(), nonlintest()

Examples

test_res <- nonlintest(data = CVD$cvd, n.lag = 4, n.boot = 1000)
# Recommended:
autoplot(test_res)
# Still works, but deprecated:
plot(test_res)

Description

Soft-deprecated in favour of autoplot.nsCosinor(), which is the same ggplot — just nudges the recommended idiom so users can + theme_bw() etc.

Usage

## S3 method for class 'nsCosinor'
plot(x, ...)

Arguments

Details

The code produces the season(s) and trend estimates.

Value

a plot of class ggplot.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.nsCosinor(), nscosinor()

Examples

# model to fit an annual pattern to the monthly cardiovascular disease data
f <- 12
tau <- c(10,50)
## Not run: 
  res12 <- nscosinor(
    data = CVD,
    response = 'adj',
    cycles = f,
    niters = 200,
    burnin = 100,
    tau = tau
    )
# Recommended:
autoplot(res12)
# Still works, but deprecated:
plot(res12)

## End(Not run)

Description

Circular plot of a monthly variable. This circular plot can be useful for estimates of an annual seasonal pattern. Darker shades of grey correspond to larger numbers.

Usage

plotCircle(months, dp = 1, ...)

Arguments

Details

Soft-deprecated in favour of plot_circle(), which returns a ggplot object you can extend with +. The base R plot below still works.

Value

a donut-type plot of a monthly variable

Author(s)

Adrian Barnett [email protected]

See Also

plot_circle()

Examples

# Recommended:
plot_circle(months = seq(1, 12, 1), dp = 0)
# Still works, but deprecated:
plotCircle(months = seq(1, 12, 1), dp = 0)

Description

A circular plot useful for visualising monthly or weekly data.

Usage

plotCircular(
  area1,
  area2 = NULL,
  spokes = NULL,
  scale = 0.8,
  labels,
  stats = TRUE,
  dp = 1,
  clockwise = TRUE,
  spoke.col = "black",
  lines = FALSE,
  centrecirc = 0.03,
  main = "",
  xlab = "",
  ylab = "",
  pieces.col = c("white", "gray"),
  length = FALSE,
  legend = TRUE,
  auto.legend = list(x = "bottomright", fill = NULL, labels = NULL, title = ""),
  ...
)

Arguments

Details

A circular plot can be useful for spotting the shape of the seasonal pattern. This function can be used to plot any circular patterns, e.g., weekly or monthly. The number of segments will be the length of the variable area1.

The plots are also called rose diagrams, with the segments then called "petals".

Value

a circular plot, also known as "rosebud", and "nightingale" plots.

Author(s)

Adrian Barnett [email protected]

References

Fisher, N.I. (1993) Statistical Analysis of Circular Data. Cambridge University Press, Cambridge.

Examples

weekfreq <- table(round(runif(100, min = 1, max = 7)))
# weeks (random data)
daysoftheweek <- c(
  'Monday',
  'Tuesday',
  'Wednesday',
  'Thursday',
  'Friday',
  'Saturday',
  'Sunday'
)
plotCircular(area1 = weekfreq, labels = daysoftheweek, dp = 0)
# Observed number of AFL players with expected values
plotCircular(
  area1 = AFL$players,
  area2 = AFL$expected,
  scale = 0.72,
  labels = month.abb,
  dp = 0,
  lines = TRUE,
  legend = FALSE
)
plotCircular(
  area1 = AFL$players,
  area2 = AFL$expected,
  scale = 0.72,
  labels = month.abb,
  dp = 0,
  lines = TRUE,
  pieces.col = c("green", "red"),
  auto.legend = list(labels = c("Obs", "Exp"), title = "# players"),
  main = "Observed and Expected AFL players"
)
# months (dummy data)
plotCircular(
  area1 = seq(1, 12, 1),
  scale = 0.7,
  labels = month.abb,
  dp = 0
)

Description

Plots results by month. Assumes the data frame contains variables called year and month.

Usage

plotMonth(data, resp, panels = 12, ...)

Arguments

Details

Soft-deprecated in favour of plot_month(), which returns a ggplot object you can extend with +. The existing function still works.

Value

Facetted lineplot of response over time, one facet per month.

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

See Also

plot_month()

Examples

# Recommended:
  plot_month(data = CVD, resp = "cvd", panels = 12)
  # Still works, but deprecated:
  plotMonth(data = CVD, resp = "cvd", panels = 12)

Description

The default print method for a casecross() object produced by casecross().

Usage

## S3 method for class 'casecross'
print(x, ...)

Arguments

Details

Uses print.coxph().

Value

printed output of casecross().

Author(s)

Adrian Barnett [email protected]

See Also

casecross(), summary.casecross(), coxph()

Examples

# subset for example
CVDdaily <- subset(CVDdaily, date <= as.Date('1987-12-31'))
# Effect of ozone on CVD death
model1 <- casecross(
  cvd ~ o3mean + tmpd + Mon + Tue + Wed + Thu + Fri + Sat,
  data = CVDdaily
)
model1

Prints basic results from monthglm()

Description

Prints basic results from monthglm()

Usage

## S3 method for class 'monthglm'
print(x, ...)

Arguments

Value

basic summary using the glm method for print()..

Examples

mmodel <- monthglm(
  formula = cvd~1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop/100000),
  offsetmonth = TRUE
  )
mmodel

Description

Print the monthly mean estimates from a Monthmean object produced by monthmean().

Usage

## S3 method for class 'Monthmean'
print(x, digits = 1, ...)

Arguments

Value

a table of values of Month and means

Author(s)

Adrian Barnett [email protected]

See Also

monthmean

Examples

mmean <- monthmean(
  data = CVD,
  resp = 'cvd',
  offsetpop = expression(pop / 100000),
  adjmonth = 'average'
)
mmean

Description

The default print method for a nonlintest object produced by nonlintest().

Usage

## S3 method for class 'nonlintest'
print(x, ...)

Arguments

Value

summary of Results of the Non-linear Test.

Author(s)

Adrian Barnett [email protected]

See Also

nonlintest() plot.nonlintest()

Examples

## Not run: 
test_res <- nonlintest(data = CVD$cvd, n.lag = 4, n.boot = 1000)
test_res

## End(Not run)

Description

The default print method for a nsCosinor object produced by nscosinor().

Usage

## S3 method for class 'nsCosinor'
print(x, ...)

Arguments

Details

Prints out the model call, number of MCMC samples, sample size and residual summary statistics.

Value

Prints out the model call, number of MCMC samples, sample size and residual summary statistics.

Author(s)

Adrian Barnett [email protected]

See Also

nscosinor() summary.nsCosinor()

Examples

# model to fit an annual pattern to the monthly cardiovascular disease data
f <- 12
tau <- c(10,50)
## Not run: 
  res12 <- nscosinor(
    data = CVD,
    response = 'adj',
    cycles = f,
    niters = 250,
    burnin = 100,
    tau = tau
    )
res12
summary(res12)

## End(Not run)

Description

Printing a summary of a month.glm

Usage

## S3 method for class 'summary.monthglm'
print(x, ...)

Arguments

Value

summary of monthglm results. This includes, number of observations, and a table of the mean, lower, upper, z value, and p value for each of the months.

Examples

mmodel <- monthglm(
  formula = cvd~1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop/100000),
  offsetmonth = TRUE
  )
summary(mmodel)

Print a summary of an nscosinor() object

Description

Print a summary of an nscosinor() object

Usage

## S3 method for class 'summary.nsCosinor'
print(x, ...)

Arguments

Value

summary - Statistics for non-stationary cosinor based on MCMC chains.

Examples

# model to fit an annual pattern to the monthly cardiovascular disease data
f <- c(12)
tau <- c(10,50)
## Not run: 
  res12 <- nscosinor(
    data = CVD,
    response = 'adj',
    cycles = f,
    niters = 200,
    burnin = 100,
    tau = tau
    )
summary(res12)

## End(Not run)

Description

Monthly number of babies born with schizophrenia in Australia from 1930 to 1971. The national number of births and number of cases are missing for January 1960 are missing.

Usage

schz

Format

A data frame with 504 observations on the following 6 variables:

• year: year of birth

• month: month of birth

• yrmon: a combination of year and month: $year+(month-1)/12$

• NBirths: monthly number of births in Australia, used as an offset

• SczBroad: monthly number of schizophrenia births using the broad diagnostic criteria

• SOI: southern oscillation index

Source

From Prof John McGrath and colleagues, The University of Queensland, Brisbane.

Examples

plot(
  schz$yrmon,
  schz$SczBroad,
  type = 'o',
  xlab = 'Date',
  ylab = 'Number of schizophrenia births'
)

Description

Tests the residuals for any remaining seasonality.

Usage

seasrescheck(res)

Arguments

Details

Plots: i) histogram of the residuals, ii) a scatter plot against residual order, iii) the autocovariance, iv) the cumulative periodogram (see cpgram())

Value

four plots: 1. The histogram of the residuals; 2. A scatter plot against residual order; 3. The autocovariance; 4. The cumulative periodogram (see cpgram()).

Author(s)

Adrian Barnett [email protected]

Examples

# cardiovascular disease data
# (use an offset of the scaled number of days in a month)
model <- cosinor(
  cvd ~ 1,
  date = 'month',
  data = CVD,
  type = 'monthly',
  family = poisson(),
  offsetmonth = TRUE
  )
seasrescheck(resid(model))

Description

Plots a sinusoid over 0 to 2$\pi$.

Usage

sinusoid(amplitude, frequency, phase, ...)

Arguments

Details

Sinusoidal curves are useful for modelling seasonal data. A sinusoid is plotted using the equation: $A\cos(ft-P), t=0,\ldots,2 \pi$, where $A$ is the amplitude, $f$ is the frequency, $t$ is time and $P$ is the phase.

Value

plot of sinusoidal wave of a given amplitude, frequency, and phase.

Author(s)

Adrian Barnett [email protected]

References

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

Examples

sinusoid(
  amplitude = 1,
  frequency = 1,
  phase = 1
  )

Description

Monthly number of stillbirths in Australia from 1998 to 2000. It is a rare event; there are 352 stillbirths out of 60,110 births. To preserve confidentiality the day of birth has been randomly re-ordered.

Usage

stillbirth

Format

A data frame with 60,110 observations on the following 7 variables:

• dob: date of birth (year-month-day)

• year: year of birth

• month: month of birth

• yrmon: a combination of year and month: $year+(month-1)/12$

• seifa: SEIFA score, an area level measure of socioeconomic status in quintiles

• gestation: gestation in weeks

• stillborn: stillborn (yes/no); 1=Yes, 0=No

Source

From Queensland Health.

Examples

table(stillbirth$month, stillbirth$stillborn)

Description

The default summary method for an object produced by casecross(). Shows the number of control days, the average number of control days per case days, and the parameter estimates.

Usage

## S3 method for class 'casecross'
summary(object, ...)

Arguments

Value

The number of control days, the average number of control days per case days, and the parameter estimates.

Author(s)

Adrian Barnett [email protected]

See Also

casecross()

Examples

# cardiovascular disease data
# subset for example
CVDdaily <- subset(CVDdaily, date <= as.Date('1987-12-31'))
# Effect of ozone on CVD death
model1 <- casecross(
  cvd ~ o3mean + tmpd + Mon + Tue + Wed + Thu + Fri + Sat,
  data = CVDdaily
)
summary(model1)
# match on day of the week
model2 <- casecross(cvd ~ o3mean + tmpd, matchdow = TRUE, data = CVDdaily)
summary(model2)
# match on temperature to within a degree
model3 <- casecross(
  cvd ~ o3mean + Mon + Tue + Wed + Thu + Fri + Sat,
  data = CVDdaily,
  matchconf = 'tmpd',
  confrange = 1
)
summary(model3)

Description

Print and summary methods for a Cosinor object produced by cosinor(). summary() summarises the sinusoidal seasonal pattern and tests whether there is a statistically significant seasonal or circadian pattern (assuming a smooth sinusoidal pattern). The amplitude describes the average height of the sinusoid, and the phase describes the location of the peak. The scale of the amplitude depends on the link function: for logistic regression the amplitude is given on a probability scale; for Poisson regression the amplitude is given on an absolute scale. print() uses the glm method for print() on the underlying model.

Usage

## S3 method for class 'Cosinor'
summary(object, digits = 2, ...)

## S3 method for class 'Cosinor'
print(x, ...)

## S3 method for class 'summary.Cosinor'
print(x, ...)

Arguments

Value

summary.Cosinor() returns a list with the following named elements:

• n: sample size.

• amp: estimated amplitude.

• amp.scale: the scale of the estimated amplitude (empty for standard regression; "probability scale" for logistic regression; "absolute scale" for Poisson regression).

• phase: estimated peak phase on a time scale.

• lphase: estimated low phase on a time scale (half a year after/before phase).

• significant: statistically significant sinusoid (TRUE/FALSE).

• alpha: statistical significance level.

• digits: minimal number of significant digits.

• text: add explanatory text to the returned phase value (TRUE) or return a number (FALSE).

• type: type of data (yearly/monthly/weekly/hourly).

• ctable: table of regression coefficients.

print.Cosinor() and print.summary.Cosinor() are called for their side effect of printing to the console and invisibly return x.

summary of output from cosinor().

Methods (by generic)

• print(summary.Cosinor): Print a summary.Cosinor object: amplitude, phase, statistical significance, and the regression coefficient table.

Functions

• print(Cosinor): Print basic results from cosinor() using the glm print method.

Author(s)

Adrian Barnett [email protected]

See Also

cosinor() plot.Cosinor() invyrfraction() glm()

Examples

## cardiovascular disease data (offset based on number of days in...
## ...the month scaled to an average month length)
res <- cosinor(
  cvd ~ 1,
  date = 'month',
  data = CVD,
  type = 'monthly',
  family = poisson(),
  offsetmonth = TRUE
  )
res
summary(res)


## cardiovascular disease data (offset based on number of days in...
## ...the month scaled to an average month length)
res <- cosinor(
  cvd ~ 1,
  date = 'month',
  data = CVD,
  type = 'monthly',
  family = poisson(),
  offsetmonth = TRUE
  )
res
summary(res)

## hourly indoor temperature data
res <- cosinor(
  bedroom ~ 1,
  date = 'datetime',
  type = 'hourly',
  data = indoor
)
summary(res)
# to get the p-values for the sine and cosine estimates
summary(res$glm)

Description

The default summary method for an object produced by monthglm().

Usage

## S3 method for class 'monthglm'
summary(object, ...)

Arguments

Details

The estimates are the mean, 95\ p-value (comparing each month to the reference month). If Poisson regression was used then the estimates are shown as rate ratios. If logistic regression was used then the estimates are shown as odds ratios.

Value

a list with the following elements:

• n: sample size.

• month.ests: parameter estimates for the intercept and months.

• month.effect: scale of the monthly effects. RR for rate ratios, OR for odds ratios, or empty otherwise.

Author(s)

Adrian Barnett [email protected]

See Also

monthglm, plot.monthglm

Examples

mmodel <- monthglm(
  formula = cvd~1,
  data = CVD,
  family = poisson(),
  offsetpop = expression(pop/100000),
  offsetmonth = TRUE
  )
summary(mmodel)

Description

The default summary method for a nsCosinor object produced by nscosinor().

Usage

## S3 method for class 'nsCosinor'
summary(object, ...)

Arguments

Details

The amplitude describes the average height of each seasonal cycle, and the phase describes the location of the peak. The results for the phase are given in radians (0 to 2$\pi$), they can be transformed to the time scale using the invyrfraction() making sure to first divide by 2$\pi$.

The larger the standard deviation for the seasonal cycles, the greater the non-stationarity. This is because a larger standard deviation means more change over time.

Value

a list with the following elements:

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

• burnin: number of MCMC samples discarded as a burn-in.

• tau: vector of smoothing parameters, tau[1] for trend, tau[2] for 1st seasonal parameter, tau[3] for 2nd seasonal parameter, etc.

• stats: summary statistics (mean and confidence interval) for the residual standard deviation, the standard deviation for each seasonal cycle, and the amplitude and phase for each cycle.

Author(s)

Adrian Barnett [email protected]

See Also

nscosinor() plot.nsCosinor()

Examples

# model to fit an annual pattern to the monthly cardiovascular disease data
f <- c(12)
tau <- c(10,50)
## Not run: 
  res12 <- nscosinor(
    data = CVD,
    response = 'adj',
    cycles = f,
    niters = 5000,
    burnin = 1000,
    tau = tau
    )
summary(res12)
plot(res12)

## End(Not run)

Description

Estimated third order moment for a time series.

Usage

third(
  data,
  n.lag,
  centre = TRUE,
  outmax = TRUE,
  plot = lifecycle::deprecated()
)

Arguments

Details

The third-order moment is the extension of the second-order moment (essentially the autocovariance). The equation for the third order moment at lags (j,k) is: $n^{-1}\sum X_t X_{t+j} X_{t+k}$. The third-order moment is useful for testing for non-linearity in a time series, and is used by nonlintest().

Value

an object of class "third" (a list) with the following elements:

• waxis: the axis from -n.lag to n.lag.

• third: the estimated third order moment in the range -n.lag to n.lag, including the symmetries.

• n.lag: the maximum lag.

Pass the result to autoplot() to draw the contour plot.

Author(s)

Adrian Barnett [email protected]

See Also

autoplot.third()

Examples

t <- third(CVD$cvd, n.lag = 12)
autoplot(t)

Description

Tests for a seasonal pattern in Binomial data. A test of whether monthly data has a sinusoidal seasonal pattern. The test has low power compared with the cosinor() test.

Usage

wtest(cases, offset, data, alpha = 0.05)

Arguments

Value

a list with the following elements:

• test: test statistic.

• pvalue: p-value.

Author(s)

Adrian Barnett [email protected]

References

Walter, S.D., Elwood, J.M. (1975) A test for seasonality of events with a variable population at risk. British Journal of Preventive and Social Medicine 29, 18–21.

Barnett, A.G., Dobson, A.J. (2010) Analysing Seasonal Health Data. Springer. doi:10.1007/978-3-642-10748-1

Examples

# tabulate the total number of births and the number of stillbirths
freqs <- table(stillbirth$month, stillbirth$stillborn)
data <- list()
data$trials <- as.numeric(freqs[, 1] + freqs[, 2])
data$success <- as.numeric(freqs[, 2])
# test for a seasonal pattern in stillbirth
test <- wtest(
  cases = 'success',
  offset = 'trials',
  data = data
)

Description

Calculate the fraction of the year for a date variable (after accounting for leap years) or for month.

Usage

yrfraction(date, type = c("daily", "weekly", "monthly"))

Arguments

Details

Returns the fraction of the year in the range [0,1).

Value

the fraction of the year.

Author(s)

Adrian Barnett [email protected]

Examples

# create fractions for the start, middle and end of the year
date <- as.Date(c(0, 181, 364), origin = '1991-01-01')
# create fractions based on these dates
yrfraction(date)
yrfraction(1:12, type = 'monthly')