| Type: | Package |
| Title: | Seasonal Analysis of Health Data |
| Version: | 0.3.16 |
| Description: | Routines for the seasonal analysis of health data, including regression models, time-stratified case-crossover, plotting functions and residual checks, see Barnett and Dobson (2010) ISBN 978-3-642-10748-1. Thanks to Yuming Guo for checking the case-crossover code. |
| License: | GPL (≥ 3) |
| URL: | https://github.com/agbarnett/season |
| BugReports: | https://github.com/agbarnett/season/issues |
| Depends: | ggplot2 (≥ 0.9.3), MASS, R (≥ 3.0.1), survival |
| Imports: | coda, graphics, grDevices, methods, stats, stringr |
| Suggests: | dlnm, knitr, rmarkdown, testthat (≥ 3.0.0) |
| VignetteBuilder: | knitr |
| Config/testthat/edition: | 3 |
| Encoding: | UTF-8 |
| LazyData: | true |
| RoxygenNote: | 7.3.2 |
| NeedsCompilation: | no |
| Packaged: | 2025-08-21 05:28:47 UTC; barnetta |
| Author: | Adrian Barnett [aut, cre], Peter Baker [aut], Oliver Hughes [aut] |
| Maintainer: | Adrian Barnett <a.barnett@qut.edu.au> |
| Repository: | CRAN |
| Date/Publication: | 2025-08-21 07:50:02 UTC |
season: Tools for Uncovering and Estimating Seasonal Patterns.
Description
The package contains graphical methods for displaying seasonal data andregression models for detecting and estimating seasonal patterns.
Details
The regression models can be applied to normal, Poisson or binomialdependent data distributions. Tools are available for both time series data(equally spaced in time) and survey data (unequally spaced in time).
Sinusoidal (parametric) seasonal patterns are available(cosinor,nscosinor), as well as models formonthly data (monthglm), and the case-crossover method tocontrol for seasonality (casecross).
season aims to fill an important gap in the software by providing arange of tools for analysing seasonal data. The examples are based on healthdata, but the functions are equally applicable to any data with a seasonalpattern.
Author(s)
Adrian Barnett <a.barnett@qut.edu.au>
Peter Baker
Oliver Hughes
Maintainer: Adrian Barnett <a.barnett@qut.edu.au>
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
Australian Football League (AFL) Players' Birthdays for the 2009 Season
Description
The data are: a) the monthly frequencies of birthdays and an expected numberbased on monthly birth statistics for 1975 to 1991. b) all 617 players'initials and birthdays (excluding non-Australian born players).
Usage
AFLFormat
A list with the following 5 variables.
- month
integer month (1 to 12)
- players
numberof players born in each month (12 observations)
- expected
expected number of players born in each month (12observations)
- 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
data(AFL)barplot(AFL$players, names.arg=month.abb)Cardiovascular Deaths in Los Angeles, 1987–2000
Description
Monthly number of deaths from cardiovascular disease in people aged 75 andover in Los Angeles for the years 1987 to 2000.
Usage
CVDFormat
A data frame with 168 observations on the following 8 variables.
- year
year of death
- month
month ofdeath
- 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 (thisvalue is constant as only one year was available, but in general thepopulation will (of course) change over time)
- ndaysmonth
number of days in each month (used as an offset)
- adj
adjusted number of CVD deaths per month using astandardised 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 EffectsInstitute, Cambridge MA.
Examples
data(CVD)plot(CVD$yrmon, CVD$cvd, type='o', xlab='Date', ylab='Number of CVD deaths per month')Daily Cardiovascular Deaths in Los Angeles, 1987–2000
Description
Daily number of deaths from cardiovascular disease in people aged 75 andover in Los Angeles for the years 1987 to 2000.
Usage
CVDdailyFormat
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
dailymean temperature (degrees Fahrenheit)
- o3mean
daily meanozone (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 EffectsInstitute, Cambridge MA.
Examples
data(CVDdaily)plot(CVDdaily$date, CVDdaily$cvd, type='p', xlab='Date', ylab='Number of CVD deaths')Amplitude Adjusted Fourier Transform (AAFT)
Description
Generates random linear surrogate data of a time series with the samesecond-order properties.
Usage
aaft(data, nsur)Arguments
data | a vector of equally spaced numeric observations (time series). |
nsur | the number of surrogates to generate (1 or more). |
Details
The AAFT uses phase-scrambling to create a surrogate of the time series thathas a similar spectrum (and hence similar second-order statistics). The AAFTis useful for testing for non-linearity in a time series, and is used bynonlintest.
Value
surrogates | a matrix of the |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Kugiumtzis D (2000) Surrogate data test for nonlinearityincluding monotonic transformations,Phys. Rev. E, vol 62
Examples
data(CVD)surr = aaft(CVD$cvd, nsur=1)plot(CVD$cvd, type='l')lines(surr[,1], col='red')Case–crossover Analysis to Control for Seasonality
Description
Fits a time-stratified case–crossover to regularly spaced time series data.The function is not suitable for irregularly spaced individual data. Thefunction only uses a time-stratified design, and other designs such as thesymmetric bi-directional design, are not available.
Usage
casecross( formula, data, exclusion = 2, stratalength = 28, matchdow = FALSE, usefinalwindow = FALSE, matchconf = "", confrange = 0, stratamonth = FALSE)Arguments
formula | formula. The dependent variable should be an integer count(e.g., daily number of deaths). |
data | data set as a data frame. |
exclusion | exclusion period (in days) around cases, set to 2(default). Must be greater than or equal to zero and smaller than |
stratalength | length of stratum in days, set to 28 (default). |
matchdow | match case and control days using day of the week(TRUE/default=FALSE). This matching is in addition to the strata matching. |
usefinalwindow | use the last stratum in the time series, which islikely to contain less days than all the other strata (TRUE/default=FALSE). |
matchconf | match case and control days using an important confounder(optional; must be in quotes). |
confrange | range of the confounder within which case and control dayswill be treated as a match (optional). Range = |
stratamonth | use strata based on months, default=FALSE. Instead of afixed strata size when using |
Details
The case–crossover method compares “case” days when events occurred(e.g., deaths) with control days to look for differences in exposure thatmight explain differences in the number of cases. Control days are selectedto be nearby to case days, which means that only recent changes in theindependent variable(s) are compared. By only comparing recent values, anylong-term or seasonal variation in the dependent and independent variable(s)can be eliminated. This elimination depends on the definition of nearby andon 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. Thestratum is controlled bystratalength, the default value is 28 days,so that cases and controls are compared in four week sections. Smallerstratum lengths provide a closer control for season, but reduce theavailable number of controls. Control days that are close to the case daymay have similar levels of the independent variable(s). To reduce thiscorrelation it is possible to place anexclusion around the cases.The default is 2, which means that the smallest gap between a case andcontrol will be 3 days.
To remove any confounding by day of the week it is possible to additionallymatch by day of the week (matchdow), although this usually reducesthe number of available controls. This matching is in addition to the stratamatching.
It is possible to additionally match case and control days by an importantconfounder (matchconf) in order to remove its effect. Control daysare matched to case days if they are: i) in the same strata, ii) have thesame day of the week ifmatchdow=TRUE, iii) have a value ofmatchconf that is within plus/minusconfrange of the value ofmatchconf on the case day. If the range is set too narrow then thenumber of available controls will become too small, which in turn means thenumber of case days with at least one control day is compromised.
The method uses conditional logistic regression (seecoxph andso the parameter estimates are odds ratios.)
The code assumes that the data frame contains a date variable (inDate format) called ‘date’.
Value
call | the original call to the casecross function. |
c.model | conditional logistic regression model of class |
ncases | total number of cases. |
ncasedays | number of case dayswith at least one control day. |
ncontroldayss | average number ofcontrol days per case day. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Janes, H., Sheppard, L., Lumley, T. (2005) Case-crossoveranalyses of air pollution exposure data: Referent selection strategies andtheir implications for bias.Epidemiology 16(6), 717–726.
Barnett, A.G., Dobson, A.J. (2010)Analysing Seasonal Health Data.Springer.
See Also
summary.casecross,coxph
Examples
# cardiovascular disease datadata(CVDdaily)CVDdaily = subset(CVDdaily, date<=as.Date('1987-12-31')) # subset for example# Effect of ozone on CVD deathmodel1 = casecross(cvd ~ o3mean+tmpd+Mon+Tue+Wed+Thu+Fri+Sat, data=CVDdaily)summary(model1)# match on day of the weekmodel2 = casecross(cvd ~ o3mean+tmpd, matchdow=TRUE, data=CVDdaily)summary(model2)# match on temperature to within a degreemodel3 = casecross(cvd ~ o3mean+Mon+Tue+Wed+Thu+Fri+Sat, data=CVDdaily, matchconf='tmpd', confrange=1)summary(model3)Mean and Confidence Interval for Circular Phase
Description
Calculates the mean and confidence interval for the phase based on a chainof MCMC samples.
Usage
ciPhase(theta, alpha = 0.05)Arguments
theta | chain of Markov chain Monte Carlo (MCMC) samples of the phase. |
alpha | the confidence level (default = 0.05 for a 95% confidenceinterval). |
Details
The estimates of the phase are rotated to have a centre of\pi, thepoint on the circumference of a unit radius circle that is furthest fromzero. The mean and confidence interval are calculated on the rotated values,then the estimates are rotated back.
Value
mean | the estimated mean phase. |
lower | the estimatedlower limit of the confidence interval. |
upper | the estimated upperlimit of the confidence interval. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
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.
Examples
theta = rnorm(n=2000, mean=0, sd=pi/50) # 2000 normal samples, centred on zerohist(theta, breaks=seq(-pi/8, pi/8, pi/30))ciPhase(theta)Cosinor Regression Model for Detecting Seasonality in Yearly Data orCircadian Patterns in Hourly Data
Description
Fits a cosinor model as part of a generalized linear model.
Usage
cosinor( formula, date, data, family = gaussian(), alpha = 0.05, cycles = 1, rescheck = FALSE, type = "daily", offsetmonth = FALSE, offsetpop = NULL, text = TRUE)Arguments
formula | regression formula. |
date | a date variable if type=“daily”, or an integer between 1and 53 if type=“weekly”, or an integer between 1 and 12 iftype=“monthly”, or aPOSIXct date if type=“hourly”. |
data | data set as a data frame. |
family | a description of the error distribution and link function tobe used in the model. Available link functions: identity, log, logit,cloglog. Note, it must have the parentheses. |
alpha | significance level, set to 0.05 (default). |
cycles | number of seasonal cycles per year if type=“daily”,“weekly” or “monthly”; number of cycles per 24 hours iftype=“hourly” |
rescheck | plot the residual checks (TRUE/FALSE), see |
type | “daily” for daily data (default), or “weekly” forweekly data, or “monthly” for monthly data, or “hourly” forhourly data. |
offsetmonth | include an offset to account for the uneven number ofdays in the month (TRUE/FALSE). Should be used for monthly counts(type=“monthly”) (with |
offsetpop | include an offset for the population (optional), thisshould be a variable in the data frame. Do not log-transform this offset, asthe transform is applied by the code. |
text | add explanatory text to the returned phase value (TRUE) orreturn a number (FALSE). Passed to the |
Details
The cosinor model captures a seasonal pattern using a sinusoid. It istherefore suitable for relatively simple seasonal patterns that aresymmetric and stationary. The default is to fit an annual seasonal pattern(cycle=1), but other higher frequencies are possible (e.g., twice peryear:cycle=2). The model is fitted using a sine and cosine term thattogether describe the sinusoid. These parameters are added to a generalizedlinear model, so the model can be fitted to a range of dependent data (e.g.,Normal, Poisson, Binomial). Unlike thenscosinor model, the cosinormodel can be applied to unequally spaced data.
Value
Returns an object of class “Cosinor” with the followingparts:
call | the original call to the cosinor function. |
glm | anobject of class |
fitted | fittedvalues for intercept and cosinor only (ignoring other independentvariables). |
fitted.plus | standard fitted values, including all otherindependent variables. |
residuals | residuals. |
date | name of thedate variable (in Date format when type=‘daily’). |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
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)data(CVD)res = cosinor(cvd~1, date='month', data=CVD, type='monthly', family=poisson(), offsetmonth=TRUE)summary(res)seasrescheck(res$residuals) # check the residuals## stillbirth datadata(stillbirth)res = cosinor(stillborn~1, date='dob', data=stillbirth, family=binomial(link='cloglog'))summary(res)plot(res)## hourly indoor temperature datares = cosinor(bedroom~1, date='datetime', type='hourly', data=indoor)summary(res)# to get the p-values for the sine and cosine estimatessummary(res$glm)Creates an Adjacency Matrix
Description
Creates an adjacency matrix in a form suitable for using inBRugs orWinBUGS.
Usage
createAdj(matrix, filename = "Adj.txt", suffix = NULL)Arguments
matrix | square matrix with 1's for neighbours and NA's fornon-neighbours. |
filename | filename that the adjacency matrix file will be written to(default=‘Adj.txt’). |
suffix | string to be appended to ‘num’, ‘adj’ and‘weights’ object names |
Details
Adjacency matrices are used by conditional autoregressive (CAR) models tosmooth estimates according to some neighbourhood map. The basic idea isthat neighbouring areas have more in common than non-neighbouring areas andso will be positively correlated.
As well as correlations in space it is possible to use CAR models to modelsimilarities in time.
In this case the matrix represents those time points that we wish to assumeto be correlated.
Value
Creates a text file namedfilename that contains the totalnumber of neighbours (num), the index number of the adjacent neighbours(adj) and the weights (weights).
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
# Nearest neighbour matrix for 5 time pointsx = c(NA,1,NA,NA,NA)(V = toeplitz(x))createAdj(V)Exercise Data from Queensland, 2005–2007
Description
Exercise data in longitudinal format from a physical activity interventionstudy in Logan, Queensland. Some subjects were lost to follow-up, so allthree visits are not available for all subjects.
Usage
exerciseFormat
A data frame with 1302 observations on the following 7 variables.
- id
subject number
- visit
visitnumber (1, 2 or 3)
- date
date of interview (year-month-day)
- year
year of interview
- month
month ofinterview
- 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 ofQueensland, Brisbane.
References
Eakin E, et al (2009) Telephone counselling for physicalactivity and diet in type 2 diabetes and hypertension,Am J of PrevMed, vol 36, pages 142–9
Examples
data(exercise)boxplot(exercise$walking ~ exercise$month)Count the Number of Days in the Month
Description
Counts the number of days per month given a range of dates. Used to adjustmonthly count data for the at-risk time period. For internal use only.
Usage
flagleap(data, report = TRUE, matchin = FALSE)Arguments
data | data. |
report | produce a brief report on the range of time used(default=TRUE). |
matchin | expand the result to match the start and end dates, otherwiseonly dates in the data will be returned (default=FALSE). |
Details
The data should contain the numeric variable called ‘year’ as a 4digit year (e.g., 1973).
Value
year | year (4 digits). |
month | month (2 digits). |
ndaysmonth | number of days in the month (either 28, 29, 30 or 31). |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Indoor Temperature Data
Description
The data are indoor temperatures (in degrees C) for a bedroom and livingroom in a house in Brisbane, Australia for the dates 10 July 2013 to 3October 2013. Temperatures were recorded using data loggers which recordedevery hour to the nearest 0.5 degrees.
Format
Adata.frame with the following 3 variables.
- datetime
date and time in
POSIXltformat- living
the living room temperature
- bedroom
the bedroom temperature
Source
Adrian G Barnett.
Examples
data(indoor)res = cosinor(bedroom~1, date='datetime', type='hourly', data=indoor)summary(res)Inverse Fraction of the Year or Hour
Description
Inverts a fraction of the year or hour to a useful time scale.
Usage
invyrfraction(frac, type = "daily", text = TRUE)Arguments
frac | a vector of fractions of the year, all between 0 and 1. |
type | “ |
text | add an explanatory text to the returned value (TRUE) or return anumber (FALSE). |
Details
Returns the day and month (fordaily) or fraction of the month (formonthly) given a fraction of the year. Assumes a year length of365.25 days fordaily. When usingmonthly the 1st of Januaryis 1, the 1st of December is 12, and the 31st of December is 12.9. Forhourly it returns the fraction of the 24-hour clock starting fromzero (midnight).
Value
daym | date (day and month for |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
invyrfraction(c(0, 0.5, 0.99), type='daily')invyrfraction(c(0, 0.5, 0.99), type='monthly')invyrfraction(c(0, 0.5, 0.99), type='hourly')Forward and Backward Sweep of the Kalman Filter
Description
Internal function to do a forward and backward sweep of the Kalman filter,used bynscosinor. For internal use only.
Usage
kalfil(data, f, vartheta, w, tau, lambda, cmean)Arguments
data | a data frame. |
f | vector of cycles in units of time. |
vartheta | variance for noise. |
w | variance of seasonal component. |
tau | controls flexibility of trend and season. |
lambda | distance between observations. |
cmean | used to give a vague prior for the starting values. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Fit a GLM with Month
Description
Fit a generalized linear model with a categorical variable of month.
Usage
monthglm( formula, data, family = gaussian(), refmonth = 1, monthvar = "month", offsetmonth = FALSE, offsetpop = NULL)Arguments
formula | regression model formula, e.g., |
data | a data frame. |
family | a description of the error distribution and link function tobe used in the model (default= |
refmonth | reference month, must be between 1 and 12 (default=1 forJanuary). |
monthvar | name of the month variable which is either an integer (1 to12) or a character or factor (‘Jan’ to ‘Dec’ or ‘January’ to ‘December’)(default='month'). |
offsetmonth | include an offset to account for the uneven number ofdays in the month (TRUE/FALSE). Should be used for monthly counts (with |
offsetpop | include an offset for the population (optional), thisshould be a variable in the data frame. Do not log-transform the offset asthe log-transform is applied by the function. |
Details
Month is fitted as a categorical variable as part of a generalized linearmodel. Other independent variables can be added to the right-hand side offormula.
This model is useful for examining non-sinusoidal seasonal patterns. Forsinusoidal seasonal patterns seecosinor.
The data frame should contain the integer months and the year as a 4 digitnumber. These are used to calculate the number of days in each monthaccounting for leap years.
Value
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 Barnetta.barnett@qut.edu.au
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
See Also
summary.monthglm,plot.monthglm
Examples
data(CVD)mmodel = monthglm(formula=cvd~1 ,data=CVD, family=poisson(), offsetpop=expression(pop/100000), offsetmonth=TRUE)summary(mmodel)Monthly Means
Description
Calculate the monthly mean or adjusted monthly mean for count data.
Usage
monthmean(data, resp, offsetpop = NULL, adjmonth = FALSE)Arguments
data | data set as a data frame. |
resp | response variable in the data set for which the means will becalculated. |
offsetpop | optional population, used as an offset (default=NULL). |
adjmonth | adjust monthly counts and scale to a 30 day month(‘ |
Details
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 adjustthe mean to account for the unequal number of days in the month (e.g., 31 inJanuary and 28 or 29 in February).
This function assumes that the data set (data) contains variables forthe year and month called year and month, respectively.
Value
Returns an object of class “Monthmean” with the followingparts:
mean | a vector of length 12 with the monthly means. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
See Also
plot.Monthmean
Examples
# cardiovascular disease datadata(CVD)mmean = monthmean(data=CVD, resp='cvd', offsetpop=expression(pop/100000), adjmonth='average')mmeanplot(mmean)Test of Non-linearity of a Time Series
Description
A bootstrap test of non-linearity in a time series using the third-ordermoment.
Usage
nonlintest(data, n.lag, n.boot, alpha = 0.05)Arguments
data | a vector of equally spaced numeric observations (time series). |
n.lag | the number of lags tested using the third-order moment, maximum= length of time series. |
n.boot | the number of bootstrap replications (suggested minimum of100; 1000 or more would be better). |
alpha | statistical significance level of test (default=0.05). |
Details
The test usesaaft to create linear surrogates with the samesecond-order properties, but no (third-order) non-linearity. The third-ordermoments (third) of these linear surrogates and the actual series arethen compared from lags 0 up ton.lag (excluding the skew at theco-ordinates (0,0)). The bootstrap test works on the overall area outsidethe limits, and gives an indication of the overall non-linearity. The plotusingregion shows those co-ordinates of the third order moment thatexceed the null hypothesis limits, and can be a useful clue for guessing thetype of non-linearity. For example, a large value at the co-ordinates (0,1)might be caused by a bi-linear seriesX_t=\alphaX_{t-1}\varepsilon_{t-1} +\varepsilon_t.
Value
Returns an object of class “nonlintest” with the followingparts:
region | the region of the third order moment where the testexceeds the limits (up to |
n.lag | the maximum lagtested using the third-order moment. |
stats | a list of followingstatistics for the area outside the test limits: |
outside | the totalarea outside of limits (summed over the whole third-order moment). |
stan | the total area outside the limits divided by its standarddeviation to give a standardised estimate. |
median | the median areaoutside the test limits. |
upper | the (1- |
pvalue | Bootstrap p-value of the areaoutside the limits to test if the series is linear. |
test | Reject thenull hypothesis that the series is linear (TRUE/FALSE). |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Barnett AG & Wolff RC (2005) A Time-Domain Test for Some Typesof Nonlinearity,IEEE Transactions on Signal Processing, vol 53,pages 26–33
See Also
print.nonlintest,plot.nonlintest
Examples
data(CVD)## Not run: test.res = nonlintest(data=CVD$cvd, n.lag=4, n.boot=1000)Non-stationary Cosinor
Description
Decompose a time series using a non-stationary cosinor for the seasonalpattern.
Usage
nscosinor( data, response, cycles, niters = 1000, burnin = 500, tau, 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 twelvemonth pattern |
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 1stseasonal parameter, tau[3] for 2nd seasonal parameter, etc. Larger values oftau allow more change between observations and hence a greater potentialflexibility in the trend and season. |
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 confidenceintervals. |
Details
This model is designed to decompose an equally spaced time series into atrend, season(s) and noise. A seasonal estimate is estimated ass_t=A_t\cos(\omega_t-P_t), wheret is time,A_t is thenon-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 thismodel gives potentially very flexible seasonal estimates.
The frequency of the seasonal estimate(s) are controlled bycycle.The cycles should be specified in units of time. If the data is monthly,then settinglambda=1/12 andcycles=12 will fit an annualseasonal pattern. If the data is daily, then settinglambda=1/365.25 andcycles=365.25 will fit an annual seasonalpattern. Specifyingcycles=c(182.6,365.25) will fit twoseasonal patterns, one with a twice-annual cycle, and one with an annualcycle.
The estimates are made using a forward and backward sweep of the Kalmanfilter. Repeated estimates are made using Markov chain Monte Carlo (MCMC).For this reason the model can take a long time to run. To give stableestimates a reasonably long sample should be used (niters), and thepossibly poor initial estimates should be discarded (burnin).
Value
Returns an object of class “nsCosinor” with the followingparts:
call | the original call to the nscosinor function. |
time | the year and month for monthly data. |
trend | mean trendand 95% confidence interval. |
season | mean season(s) and 95%confidence interval(s). |
oseason | overall season(s) and 95%confidence interval(s). This will be the same as |
fitted | fitted values and 95% confidenceinterval, based on trend + season(s). |
residuals | residuals based onmean trend and season(s). |
n | the length of the series. |
chains | MCMC chains (of class mcmc) of variance estimates: standarderror for overall noise (std.error), standard error for season(s)(std.season), phase(s) and amplitude(s) |
cycles | vector of cycles inunits of time. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
Barnett, A.G., Dobson, A.J. (2004) Estimating trends and seasonality incoronary heart diseaseStatistics in Medicine. 23(22) 3505–23.
See Also
plot.nsCosinor,summary.nsCosinor
Examples
data(CVD)# model to fit an annual pattern to the monthly cardiovascular disease dataf = 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)Initial Values for Non-stationary Cosinor
Description
Creates initial values for the non-stationary cosinor decompositionnscosinor. For internal use only.
Arguments
data | a data frame. |
response | response variable. |
tau | vector of smoothing parameters, tau[1] for trend, tau[2] for 1stseasonal parameter, tau[3] for 2nd seasonal parameter, etc. Larger values oftau allow more change between observations and hence a greater potentialflexibility in the trend and season. |
lambda | distance between observations (lambda=1/12 for monthly data,default). |
n.season | number of seasons. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Periodogram
Description
Estimated periodogram using the fast Fourier transform (fft).
Usage
peri(data, adjmean = TRUE, plot = TRUE)Arguments
data | a data frame. |
adjmean | subtract the mean from the series before calculating theperiodogram (default=TRUE). |
plot | plot the estimated periodogram (default=TRUE). |
Value
peri | periodogram, I( |
f | frequencies inradians, |
c | frequencies in cycles of time, |
amp | amplitude periodogram. |
phase | phaseperiodogram. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
data(CVD)p = peri(CVD$cvd)Phase from Cosinor Estimates
Description
Calculate the phase given the estimated sine and cosine values from acosinor model.
Usage
phasecalc(cosine, sine)Arguments
cosine | estimated cosine value from a cosinor model. |
sine | estimated sine value from a cosinor model. |
Details
Returns the phase in radians, in the range[0,2\pi). The phase is thepeak in the sinusoid.
Value
phaser | Estimated phase in radians. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Fisher, N.I. (1993)Statistical Analysis of CircularData. Cambridge University Press, Cambridge. Page 31.
Examples
phasecalc(cosine=0, sine=1) # pi/2Plot the Results of a Cosinor
Description
Plots the fitted sinusoid from aCosinor object produced bycosinor.
Usage
## S3 method for class 'Cosinor'plot(x, ...)Arguments
x | a |
... | additional arguments passed to the sinusoid plot. |
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 monthlydata then month is shown on the x-axis. If the analysis was based on dailydata then time is shown on the x-axis.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
cosinor,summary.Cosinor,seasrescheck
Plot of Monthly Mean Estimates
Description
Plots estimated monthly means.
Usage
## S3 method for class 'Monthmean'plot(x, ...)Arguments
x | a |
... | additional arguments passed to the plot. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
monthmean
Plot of Monthly Estimates
Description
Plots the estimated from a generalized linear model with a categoricalvariable of month.
Usage
## S3 method for class 'monthglm'plot(x, alpha = 0.05, ylim = NULL, ...)Arguments
x | a |
alpha | statistical significance level of confidence intervals. |
ylim | y coordinates ranges (the default is NULL, and the limits areautomatically calculated). |
... | additional arguments passed to the plot. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
monthglm
Examples
data(CVD)mmodel = monthglm(formula=cvd~1, data=CVD, family=poisson(), offsetpop=expression(pop/100000), offsetmonth=TRUE, refmonth=6)plot(mmodel)Plot the Results of the Non-linear Test
Description
Creates a contour plot of the region of the third-order moment outside thetest limits generated bynonlintest.
Usage
## S3 method for class 'nonlintest'plot(x, plot = TRUE, ...)Arguments
x | a |
plot | display plot (TRUE) or return plot (FALSE) |
... | additional arguments to |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
nonlintest
Plot the Results of a Non-stationary Cosinor
Description
Plots the trend and season(s) from ansCosinor object produced bynscosinor.
Usage
## S3 method for class 'nsCosinor'plot(x, ...)Arguments
x | a |
... | further arguments passed to or from other methods. |
Details
The code produces the season(s) and trend estimates.
Value
gplot | A plot of class |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
nscosinor
Print the Results of a Cosinor
Description
The default print method for aCosinor object produced bycosinor.
Usage
## S3 method for class 'Cosinor'print(x, ...)Arguments
x | a |
... | optional arguments to |
Details
Usesprint.glm.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
cosinor,summary.Cosinor,glm
Print the Results from Monthmean
Description
Print the monthly means from aMonthmean object produced bymonthmean.
Usage
## S3 method for class 'Monthmean'print(x, digits = 1, ...)Arguments
x | a |
digits | minimal number of significant digits, see |
... | additional arguments passed to the print. |
Details
The code prints the monthly mean estimates.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
monthmean
Print the Results of a Case-Crossover Model
Description
The default print method for acasecross object produced bycasecross.
Usage
## S3 method for class 'casecross'print(x, ...)Arguments
x | a |
... | optional arguments to |
Details
Usesprint.coxph.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
casecross,summary.casecross,coxph
Printmonthglm
Description
Printmonthglm
Usage
## S3 method for class 'monthglm'print(x, ...)Arguments
x | Object of class |
... | further arguments passed to or from other methods. |
Print the Results of the Non-linear Test
Description
The default print method for anonlintest object produced bynonlintest.
Usage
## S3 method for class 'nonlintest'print(x, ...)Arguments
x | a |
... | additional arguments to |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
nonlintest,plot.nonlintest
Print the Results of a Non-stationary Cosinor
Description
The default print method for ansCosinor object produced bynscosinor.
Usage
## S3 method for class 'nsCosinor'print(x, ...)Arguments
x | a |
... | further arguments passed to or from other methods. |
Details
Prints out the model call, number of MCMC samples, sample size and residualsummary statistics.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
nscosinor,summary.nsCosinor
printing a summary of a Cosinor
Description
printing a summary of a Cosinor
Usage
## S3 method for class 'summary.Cosinor'print(x, ...)Arguments
x | a |
... | optional arguments to |
printing a summary of a month.glm
Description
printing a summary of a month.glm
Usage
## S3 method for class 'summary.monthglm'print(x, ...)Arguments
x | a |
... | further arguments passed to or from other methods. |
printing a summary of an nscosinor
Description
printing a summary of an nscosinor
Usage
## S3 method for class 'summary.nscosinor'print(x, ...)Arguments
x | a |
... | further arguments passed to or from other methods. |
Random Inverse Gamma Distribution
Description
Internal function to simulate a value from an inverse Gamma distribution,used bynscosinor. See the MCMCpack library. For internal use only.
Arguments
n | number of observations. |
shape | Gamma shape parameter. |
scale | Gamma scale parameter (default=1). |
Schizophrenia Births in Australia, 1930–1971
Description
Monthly number of babies born with schizophrenia in Australia from 1930 to1971. The national number of births and number of cases are missing forJanuary 1960 are missing.
Usage
schzFormat
A data frame with 504 observations on the following 6 variables.
- year
year of birth
- month
month ofbirth
- yrmon
a combination of year and month:
year+(month-1)/12- NBirths
monthly number of births inAustralia, used as an offset
- SczBroad
monthly number ofschizophrenia births using the broad diagnostic criteria
- SOI
southern oscillation index
Source
From Prof John McGrath and colleagues, The University of Queensland,Brisbane.
Examples
data(schz)plot(schz$yrmon, schz$SczBroad, type='o', xlab='Date', ylab='Number of schizophrenia births')Seasonal Residual Checks
Description
Tests the residuals for any remaining seasonality.
Usage
seasrescheck(res)Arguments
res | residuals from some time series regression model. |
Details
Plots: i) histogram of the residuals, ii) a scatter plot against residualorder, iii) the autocovariance, iv) the cumulative periodogram (seecpgram)
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
# cardiovascular disease data# (use an offset of the scaled number of days in a month)data(CVD)model = cosinor(cvd~1, date='month', data=CVD, type='monthly', family=poisson(), offsetmonth=TRUE)seasrescheck(resid(model))Plot a Sinusoid
Description
Plots a sinusoid over 0 to 2\pi.
Usage
sinusoid(amplitude, frequency, phase, ...)Arguments
amplitude | the amplitude of the sinusoid (its maximum value). |
frequency | the frequency of the sinusoid in 0 to 2 |
phase | the phase of the sinusoid (location of the peak). |
... | additional arguments passed to the plot. |
Details
Sinusoidal curves are useful for modelling seasonal data. A sinusoid isplotted using the equation:A\cos(ft-P), t=0,\ldots,2 \pi, whereA is the amplitude,f is the frequency,t is time andP is the phase.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Barnett, A.G., Dobson, A.J. (2010)Analysing SeasonalHealth Data. Springer.
Examples
sinusoid(amplitude=1, frequency=1, phase=1)Stillbirths in Queensland, 1998–2000
Description
Monthly number of stillbirths in Australia from 1998 to 2000. It is a rareevent; there are 352 stillbirths out of 60,110 births. To preserveconfidentiality the day of birth has been randomly re-ordered.
Usage
stillbirthFormat
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 levelmeasure of socioeconomic status in quintiles
- gestation
gestation in weeks
- stillborn
stillborn (yes/no); 1=Yes, 0=No
Source
From Queensland Health.
Examples
data(stillbirth)table(stillbirth$month, stillbirth$stillborn)Summary of the Results of a Case-crossover Model
Description
The default summary method for acasecross object produced bycasecross.
Usage
## S3 method for class 'casecross'summary(object, ...)Arguments
object | a |
... | further arguments passed to or from other methods. |
Details
Shows the number of control days, the average number of control days percase days, and the parameter estimates.
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
casecross
Summary for a Monthglm
Description
The default summary method for amonthglm object produced bymonthglm.
Usage
## S3 method for class 'monthglm'summary(object, ...)Arguments
object | a |
... | further arguments passed to or from other methods. |
Details
The estimates are the mean, 95% confidence interval, Z-value and associatedp-value (comparing each month to the reference month). If Poisson regressionwas used then the estimates are shown as rate ratios. If logistic regressionwas used then the estimates are shown as odds ratios.
Value
n | sample size. |
month.ests | parameter estimates for theintercept and months. |
month.effect | scale of the monthly effects.‘RR’ for ‘rate ratios’, ‘OR’ for ‘odds ratios’,or empty otherwise. |
x | object of class |
object | object of class |
list() | furtherarguments passed to or from other methods. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
monthglm,plot.monthglm
Summary for a Non-stationary Cosinor
Description
The default summary method for ansCosinor object produced bynscosinor.
Usage
## S3 method for class 'nsCosinor'summary(object, ...)Arguments
object | a |
... | further arguments passed to or from other methods. |
Details
The amplitude describes the average height of each seasonal cycle, and thephase describes the location of the peak. The results for the phase aregiven in radians (0 to 2\pi), they can be transformed to the timescale using theinvyrfraction making sure to first divide by2\pi.
The larger the standard deviation for the seasonal cycles, the greater thenon-stationarity. This is because a larger standard deviation means morechange over time.
Value
cycles | vector of cycles in units of time, e.g., for a six andtwelve month pattern |
niters | total number ofMCMC samples. |
burnin | number of MCMC samples discarded as a burn-in. |
tau | vector of smoothing parameters, tau[1] for trend, tau[2] for 1stseasonal parameter, tau[3] for 2nd seasonal parameter, etc. |
stats | summary statistics (mean and confidence interval) for theresidual standard deviation, the standard deviation for each seasonal cycle,and the amplitude and phase for each cycle. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
See Also
nscosinor,plot.nsCosinor
Third-order Moment
Description
Estimated third order moment for a time series.
Usage
third(data, n.lag, centre = TRUE, outmax = TRUE, plot = TRUE)Arguments
data | a vector of equally spaced numeric observations (time series). |
n.lag | the number of lags, maximum = length of time series. |
centre | centre series by subtracting mean (default=TRUE). |
outmax | display the (x,y) lag co-ordinates for the maximum and minimumvalues (default=TRUE). |
plot | contour plot of the third order moment (default=TRUE). |
Details
The third-order moment is the extension of the second-order moment(essentially the autocovariance). The equation for the third order moment atlags (j,k) is:n^{-1}\sum X_t X_{t+j} X_{t+k}. The third-order momentis useful for testing for non-linearity in a time series, and is used bynonlintest.
Value
waxis | The axis – |
third | The estimated third order moment in the range –n.lag to n.lag,including the symmetries. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
data(CVD)third(CVD$cvd, n.lag=12)Walter and Elwood's Test of Seasonality
Description
Tests for a seasonal pattern in Binomial data.
Usage
wtest(cases, offset, data, alpha = 0.05)Arguments
cases | variable name for cases (“successes”). |
offset | variable name for at-risk population (“trials”). |
data | data frame (optional). |
alpha | significance level (default=0.05). |
Details
A test of whether monthly data has a sinusoidal seasonal pattern. The testhas low power compared with thecosinor test.
Value
test | test statistic. |
pvalue | p-value. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
References
Walter, S.D., Elwood, J.M. (1975) A test for seasonality ofevents with a variable population at risk.British Journal ofPreventive and Social Medicine 29, 18–21.
Barnett, A.G., Dobson, A.J. (2010)Analysing Seasonal Health Data.Springer.
Examples
data(stillbirth)# tabulate the total number of births and the number of stillbirthsfreqs = 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 stillbirthtest = wtest(cases='success', offset='trials', data=data)Fraction of the Year
Description
Calculate the fraction of the year for a date variable (after accounting forleap years) or for month.
Usage
yrfraction(date, type = "daily")Arguments
date | a date variable if type=‘ |
type | ‘daily’ for dates, or ‘ |
Details
Returns the fraction of the year in the range [0,1).
Value
yrfrac | Fraction of the year. |
Author(s)
Adrian Barnetta.barnett@qut.edu.au
Examples
# create fractions for the start, middle and end of the yeardate = as.Date(c(0, 181, 364), origin='1991-01-01')# create fractions based on these datesyrfraction(date)yrfraction(1:12, type='monthly')