--- title: "Getting started with actiRhythm" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Getting started with actiRhythm} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: REFERENCES.bib csl: apa.csl link-citations: true --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, dpi = 110 ) ``` You have an ActiGraph recording for one person and a practical question. Does this person have a stable daily rhythm, how strong is it, and is it real or could it be noise? This article answers that question end to end on one bundled recording. By the end you will be able to read the file, describe the rhythm without assuming a shape, test whether it is real, locate and bound its period, probe its finer structure, and write out a one-row summary you can stack across people. The method assumes a multi-day activity-count series with timestamps, which is what an ActiGraph `.agd` file holds and what counts exported from any source can supply. Every number and figure below is reproducible from the installed package. ```{r libraries, message = FALSE, warning = FALSE} library(actiRhythm) library(ggplot2) ``` ## The recording `read.agd()` opens the `.agd` SQLite file and auto-detects the epoch length and axes; `agd.counts()` returns a tidy data frame with `$axis1` (the vertical-axis counts) and `$timestamp`. One recording ships with the package, reached through `example_agd()`. ```{r read} agd <- agd.counts(read.agd(example_agd(1), verbose = FALSE)) head(agd[, c("timestamp", "axis1")]) ``` Read the actogram top to bottom before computing anything. It is the fastest way to see whether there is a rhythm to measure at all. ```{r actogram, fig.height = 5, fig.cap = "Double-plotted actogram: each row is one day shown twice across 48 hours, so a stable rhythm forms a vertical band of activity at the same clock time each day; scattered fill marks a fragmented or shifting rhythm."} plot_actogram(agd$axis1, agd$timestamp) ``` In this recording the active rows do not stack into one clean vertical band; they slide a little later or earlier from day to day, which is the low interdaily stability the next section puts a number on. ## The nonparametric rhythm `circadian.rhythm()` describes the rest-activity pattern without assuming a waveform. It returns interdaily stability (IS) and intradaily variability (IV), introduced by @witting1990, together with the relative amplitude (RA) and the least- and most-active windows L5 and M10 of @vansomeren1999, each with its onset time. ```{r rhythm} rhythm <- circadian.rhythm(agd$axis1, agd$timestamp) rhythm ``` Read the three headline numbers together. IS ranges 0 to 1 and rises as the daily pattern repeats more exactly from day to day; IV runs from about 0 (a smooth sine) to about 2 (noisy or split into ultradian bouts); RA ranges 0 to 1 and grows with the contrast between an active day and a restful night. Here a high RA (about 0.98) sits with a low IS (about 0.23): a strong day-night contrast carried on irregular timing. The header also notes that six of the eight calendar days qualified, because `circadian.rhythm()` sets aside days with too little valid recording before computing the metrics. ## Cosinor and the rhythmicity test `cosinor.analysis()` fits a 24-hour cosine and returns the MESOR (rhythm-adjusted mean), amplitude, and acrophase, the clock time of the peak [@cornelissen2014]. ```{r cosinor} cos <- cosinor.analysis(agd$axis1, agd$timestamp, period = 24) c(mesor = unname(cos$mesor), amplitude = unname(cos$amplitude), acrophase = unname(cos$acrophase)) ``` The acrophase is **circular**: it is a clock time, so 23.5 and 0.5 are an hour apart, not 23 hours apart, and it is only interpretable once the amplitude is distinguishable from zero. `rhythmicity.test()` is that check, the zero-amplitude F-test of @nelson1979, reported with the percent of variance the rhythm explains. ```{r rhythmicity} rhythmicity.test(agd$axis1, agd$timestamp, cosinor_result = cos) ``` A small p-value rejects the flat-line null. Here the test is clearly significant (p well under 0.05) and the single cosine explains a little over forty percent of the variance, so the acrophase above is worth reading. The percent-rhythm also reminds you that the rest of the variance, more than half, is a mix of noise and structure a single cosine cannot capture. Overlay the fit to see what it is and is not catching. ```{r cosinor-plot, fig.cap = "Hourly activity with the fitted 24-hour cosine overlaid; the peak of the curve is the acrophase and its height above the MESOR line is the amplitude. Systematic departures from the curve are where a single cosine is too smooth for real rest-activity data."} plot_extended_cosinor(agd$axis1, agd$timestamp) ``` Look where the fitted curve and the hourly points part company: the single cosine runs too smooth through the early hours, and that gap is the structure the percent-rhythm says one cosine leaves on the table. When you need the joint uncertainty in amplitude and acrophase together, `cosinor.confidence.ellipse(cos)` returns the @bingham1982 ellipse; an ellipse that encloses the origin is the geometric form of "no detectable rhythm". ## The period A cosinor assumes exactly 24 hours. `circadian.period()` does not: it runs a Lomb-Scargle periodogram over the 18 to 30 hour band and returns the dominant period `tau`, with a false-alarm probability in `$p_value` [@lomb1976; @scargle1982; @baluev2008]. ```{r period} per <- circadian.period(agd$axis1, agd$timestamp) c(tau = per$tau, p_value = per$p_value) ``` The strongest cycle here sits near 24.9 hours, a little longer than a solar day, and the false-alarm probability is essentially zero, so the cycle itself is not in doubt. Its exact length is less certain. `period.ci()` puts a confidence interval around that same peak with a circular block bootstrap [@kunsch1989; @politis1992], which respects the autocorrelation in activity data. The interval it returns is wide enough to span 24 hours. On a single recording the period is clearly present but not pinned to the hour. Passing `seed = 1` makes the bootstrap reproducible; the 50 replicates here keep the vignette fast, and a real interval would use more. ```{r period-ci} period.ci(agd$axis1, agd$timestamp, n_boot = 50, seed = 1) ``` ```{r periodogram, fig.cap = "Lomb-Scargle spectral power across candidate periods; the tallest peak is the dominant rhythm, and a peak rising above the significance line rejects the no-rhythm null."} plot_periodogram(agd$axis1, agd$timestamp) ``` The chi-square periodogram of @sokolove1978 asks the same question a different way, scoring each trial period by how much of the variance folds onto it. A second estimator that clears its own threshold at the same period is reassurance that the peak is real rather than an artifact of one method. ```{r chisq, fig.cap = "Chi-square (Sokolove-Bushell) periodogram: the Qp statistic against trial period, with the significance threshold. The peak should fall at the period the Lomb-Scargle periodogram found."} plot_chisq(agd$axis1, agd$timestamp) ``` ## Watching the period drift A single periodogram averages over the whole recording. `circadian.spectrogram()` slides a window across it and recomputes the periodogram in each, so you can see the dominant period move, the chronobiology equivalent of a music spectrogram. ```{r spectrogram, fig.height = 4, fig.cap = "Period (vertical) against time (horizontal), coloured by spectral power. A horizontal band of colour is a stable period; a band that bends upward or downward is a rhythm lengthening or shortening across the recording."} circadian.spectrogram(agd$axis1, agd$timestamp, step_hours = 24)$plot ``` Follow the band of colour left to right: here the dominant period wanders rather than holding steady at 24 hours. Read the final window with particular caution, since it rides the upper edge of the search band, where a weakening rhythm and a genuinely longer period can look alike. ## Nonlinear structure Two recordings can share the same IS and amplitude yet differ in their moment-to-moment dynamics. `fractal.dfa()` estimates the detrended-fluctuation exponent alpha, the long-range temporal correlation in the series [@peng1994; @hu2001]. ```{r dfa} fractal.dfa(agd$axis1)$alpha ``` An alpha near 1 is the "1/f" scaling typical of healthy physiological output; near 0.5 is uncorrelated noise; well above 1 is a smoother, more random-walk-like signal. Here alpha comes out near 0.93, close to the 1/f end, so this series carries genuine long-range correlation rather than looking like noise. The log-log fit shows whether one exponent holds across all time scales. ```{r dfa-plot, fig.cap = "Detrended fluctuation on log-log axes; the slope is the DFA exponent alpha. One straight line is monofractal scaling, a visible kink is a crossover between short- and long-range dynamics."} plot_dfa(agd$axis1) ``` Two further questions use the same series and a single call each. `multiscale.entropy()` asks how the complexity of the signal changes as you coarse-grain it across time scales [@costa2002], and `mfdfa()` asks whether one scaling exponent is enough or the series is genuinely multifractal [@kantelhardt2002]. Their help pages carry the math, and the [Beyond the basics](https://rdazadda.github.io/actiRhythm/articles/beyond-the-basics.html) article runs them on this recording; here it is enough to know they refine, not replace, the single alpha above. ## Rest-activity transitions `state.transitions()` summarises how readily the subject switches between rest and activity, the kRA (rest-to-active) and kAR (active-to-rest) rates that capture fragmentation a single amplitude cannot [@lim2011]. ```{r transitions} st <- state.transitions(agd$axis1) c(kRA = st$kRA, kAR = st$kAR) ``` For this subject kAR runs well above kRA: once active they settle back to rest quickly, but once at rest they stay there. That asymmetry is the fragmentation a single amplitude cannot see. ## Pinpointing sleep and wake `state.transitions()` gives one pair of rates for the whole recording. `sleep.changepoints()` instead locates the sleep-onset and wake-onset of each night directly from the counts, with no scored sleep required: a 24-hour cosinor bounds each rest and active span roughly, and a change point inside each bound places the precise transition [@chensun2024]. ```{r changepoints} cp <- sleep.changepoints(agd$axis1, agd$timestamp) cp ``` ```{r cp-plot, fig.height = 3.5, fig.cap = "Each detected sleep onset (blue) and wake onset (orange) by date and clock time. The detector reads the per-night sleep and wake timing directly from the counts."} cps <- cp$changepoints cps$clock <- as.numeric(format(cps$time, "%H")) + as.numeric(format(cps$time, "%M")) / 60 ggplot(cps, aes(time, clock, colour = type)) + geom_point(size = 3) + scale_colour_manual(values = c("sleep onset" = "#236192", "wake onset" = "#E69F00")) + labs(x = "Date", y = "Clock hour", colour = NULL) + theme_actiRhythm() ``` For this recording the detector finds seven nights, sleep onset late in the evening and wake onset spread from early morning to around midday, with most rest spans in the eleven-to-twelve-hour range. That is the per-night sleep and wake timing the single transition rate above cannot localise. ## Every rest bout, naps included `sleep.changepoints()` finds the one main rest bout of each cycle. `rest.periods()` takes the complementary view, consolidating every spell of low activity into a rest bout, however many a day holds, daytime naps and fragmented rest included [@roenneberg2015]. Each epoch is compared to a fraction of its own 24-hour activity level, and runs that stay below it grow into consolidated bouts. ```{r restperiods} rp <- rest.periods(agd$axis1, agd$timestamp) rp ``` ```{r restperiods-plot, fig.height = 3, fig.cap = "Each consolidated rest bout as a bar from onset to offset. Detecting more than one bout a day, including naps, is what separates this from the one-per-cycle change-point detector."} ggplot(rp$rest_periods, aes(onset, bout, xend = offset, yend = bout, colour = type)) + geom_segment(linewidth = 3) + scale_colour_manual(values = c(main = "#236192", nap = "#E69F00")) + labs(x = "Time", y = "Rest bout", colour = NULL) + theme_actiRhythm() ``` For this recording it consolidates nine rest bouts across six days, more than one a day, so a few days carry a second rest beyond the main night. `rest.crespo()` reaches a comparable result by a different route, rank-order filtering and binary morphology rather than seed-and-grow consolidation [@crespo2012]; running both gives an independent cross-check of where a recording's rest bouts fall. Neither applies a wear-time filter, so gate the counts on valid wear first if a recording has device-off stretches. ## Scoring sleep and its regularity The metrics so far describe the rhythm; the next need a per-epoch sleep or wake label. `sleep.cole.kripke()` produces that label directly from the counts, the standard count-based classifier for adults [@cole1992], with `sleep.sadeh()`, validated on adults and adolescents, as the alternative. ```{r scoring} state <- sleep.cole.kripke(agd$axis1) table(state) ``` Count-based scoring labels every low-activity epoch as sleep, so a mostly sedentary recording scores a high sleep fraction; what matters downstream is how regular that label is from day to day. `sleep.regularity.index()` measures exactly that, the probability that two epochs 24 hours apart share the same state [@phillips2017]. ```{r sri} sleep.regularity.index(state, agd$timestamp) ``` When you also have explicit sleep periods, `social.jet.lag()` returns the work-day to free-day mid-sleep difference [@wittmann2006] and `lids()` the locomotor inactivity during sleep that tracks ultradian sleep structure [@winnebeck2018]. ## Pooling the evidence Each test above answers "is there a rhythm?" in its own language. `consensus.rhythmicity()` runs three of them, the cosinor F-test, the Lomb-Scargle false-alarm probability, and the chi-square periodogram, and reports both a majority vote and a Cauchy-combined p-value [@liu2020]. ```{r consensus} consensus.rhythmicity(agd$axis1, agd$timestamp) ``` For this recording a majority of the tests find a rhythm (two of the three) and the combined p-value is effectively zero, so the evidence leans firmly toward a real rhythm even where one method on its own hesitates. The three tests run on the same series, so they are not independent; the Cauchy combination stays valid under that dependence, where Fisher's method would not. ## A one-row summary A one-row summary is then easy to assemble and `rbind()` across subjects. ```{r summary} data.frame( IS = rhythm$IS, IV = rhythm$IV, RA = rhythm$RA, tau = per$tau, amplitude = cos$amplitude, acrophase = cos$acrophase, dfa_alpha = fractal.dfa(agd$axis1)$alpha ) ``` Read across, this row is one person: a strongly amplitude-modulated rhythm whose timing is irregular from day to day, running on a period a little over 24 hours. ## A whole study at once You have now run the whole arc on one recording, ending with a single summary row. Most studies, though, are a folder of recordings rather than one, and that is what `circadian.batch()` is for. It runs this entire arc over every `.agd` file in a directory and returns one row per recording, so a cohort becomes a single data frame you can model. A file that fails to read is reported in an `error` column instead of stopping the run. The package ships two recordings, so this runs as a real, if small, batch. ```{r batch} batch <- circadian.batch( system.file("extdata", package = "actiRhythm"), verbose = FALSE ) batch[, c("file", "IS", "IV", "RA", "period_tau", "rhythm_p_value")] ``` Each row is one subject, described by the same metrics as the walkthrough above, and already the two differ: the second is steadier from day to day (higher IS) but less amplitude-modulated (lower RA), and its rhythm falls just short of significance (p = 0.05). Point the call at your own folder and nothing else changes. When you want the full analysis for a single recording rather than a summary row, `circadian.workbook()` writes every table to a multi-sheet Excel file: a summary row plus the hourly profile, both periodograms, the fluctuation and multifractal spectra, and a data dictionary. ```{r workbook, eval = FALSE} circadian.workbook(agd$axis1, agd$timestamp, file = "subject01.xlsx") ``` ## Where next - Decomposition, the multi-harmonic profile, and the finer nonparametric and phase metrics are demonstrated in [Beyond the basics](https://rdazadda.github.io/actiRhythm/articles/beyond-the-basics.html). - Reading raw `.gt3x`/`.cwa`/`.bin` files, calibration, the ENMO/MAD/z-angle metrics, and diary-free posture-based sleep are in [From raw acceleration](https://rdazadda.github.io/actiRhythm/articles/raw-pipeline.html). - When a rhythm drifts, shifts phase, or fragments across the recording, see [Nonstationary and complex rhythms](https://rdazadda.github.io/actiRhythm/articles/nonstationary-rhythms.html). - Every method family has its own walkthrough as well: cosinor models, period and spectral analysis, sleep and rest detection, phase and regularity, fractal and nonlinear dynamics, choosing a method, and validation against reference packages. The full set is in the [articles index](https://rdazadda.github.io/actiRhythm/articles/index.html). - `?circadian.rhythm`, `?cosinor.analysis`, `?circadian.period`, `?fractal.dfa`, `?consensus.rhythmicity` give the per-function help with full references. - `?circadian.batch` and `?circadian.workbook` are the study-scale entry points. - For a multi-subject design, `population.cosinor()` and `cosinor.compare()` give the mean rhythm and a two-group test (Bingham/Hotelling); for covariates, repeated measures, or nesting, feed the per-subject metrics this package emits (MESOR, amplitude, acrophase, IS, IV, RA, period) into a mixed model with `lme4`/`nlme` alongside your design. The package gives you the metrics; you supply the multilevel model. - Cite the package with `citation("actiRhythm")`, and report `packageVersion("actiRhythm")` with your results so the analysis stays reproducible across releases. - Questions and bugs: . ## References