Circadian rhythm analysis with the cosinor method: mesor, amplitude, peak time and comparing groups

Quick answer. For a measurement taken at several times of day, in different animals at each time point, fit a cosine with a 24-hour period by least squares (the cosinor model). It gives the mesor (rhythm-adjusted mean), the amplitude (half the peak-to-trough range of the fitted cosine) and the peak time (acrophase), with confidence intervals for amplitude and peak time, and a zero-amplitude F-test of whether there is a rhythm at all. Read amplitude and peak time only when that test shows a rhythm. To compare genotypes or treatments, fit one model with group × cosine and group × sine terms and test those terms to see whether the rhythm differs, then report the differences in amplitude and peak time as estimates with their intervals. A t-test at every time point is the wrong tool: with 12 time points and no real difference, at least one came out “significant” in about 46% of our simulated studies, and correcting the 12 tests for multiplicity left them far less powerful than the one model.

Terms used below. Circadian rhythm: a roughly 24-hour cycle. Least squares: fitting a model by making the squared differences between data and fitted values as small as possible. Zeitgeber time (ZT): time since lights on (ZT0) in a light–dark cycle; circadian time (CT) is the equivalent in constant darkness, where CT12 marks activity onset in nocturnal animals. Free-running period: the cycle length an animal keeps without light cues. Cosinor model: y = mesor + amplitude × cos(2π(t − peak time)/24), fitted by least squares as a linear regression on cos(2πt/24) and sin(2πt/24). Single-component cosinor: one 24-hour cosine; adding a second harmonic means adding a 12-hour cosine and sine as well. Mesor: the rhythm-adjusted mean. Amplitude: half the difference between the peak and trough of the fitted cosine. Acrophase (peak time): the time of the fitted peak. Zero-amplitude test: an F-test of whether the cosine terms improve on a flat line. F-test: compares the fit of two nested models; F(2, 66) means 2 terms tested and 66 residual degrees of freedom. Mixed model: a regression that adds a random effect for each animal, so repeated measurements of one animal are not treated as independent. Bonferroni correction: testing each of m comparisons at 0.05/m. Post hoc test: a comparison at each time point after an overall test. False discovery rate: the expected share of false positives among the genes called rhythmic. Delta method: a way to get approximate standard errors for quantities, such as amplitude, calculated from regression coefficients. 95% confidence interval (CI): the range of values compatible with the data. Power: the chance of detecting an effect that is really there. False positive: a “significant” result when there is no real effect.

Which analysis for your question

Question	Analysis	Report
Is there a 24-hour rhythm?	Cosinor zero-amplitude F-test	Mesor, amplitude and peak time with 95% CIs; F and p
Does the rhythm differ between groups?	One cosinor model with group × cosine and group × sine terms, assuming the same 24-h period; F-test for those terms	F and p; then differences in amplitude and peak time with CIs
Does the mean level differ?	The group term in the same model (mesor difference)	Mesor difference with CI
Which of thousands of genes cycle?	Rhythm detection methods for omics, such as JTK_CYCLE or RAIN, with false discovery rate control	Number of rhythmic genes at a stated false discovery rate
The same animals sampled repeatedly	A mixed-effects cosinor model with animal as a random effect, checking how correlation changes with time apart (repeated measures)	As above, with the random-effect SD

Cornelissen defines the parameters: “M is the MESOR (Midline Statistic Of Rhythm, a rhythm-adjusted mean), A is the amplitude (a measure of half the extent of predictable variation within a cycle)”, with the acrophase giving the timing of the peak, and describes the zero-amplitude F-test and the confidence intervals for amplitude and acrophase (Cornelissen 2014). For comparing groups, Parsons and colleagues’ CircaCompare parametrises the model so that it allows “the explicit estimation of rhythmic parameters” and their differences (Parsons et al. 2020). For sampling, Hughes and colleagues “recommend collecting at least 2 complete cycles of data when detecting rhythmicity” and “sampling at least 12 time points per cycle across 2 full cycles to optimize statistical power” (Hughes et al. 2017); Cornelissen cites designs “recommending at least 6 timepoints per cycle”. For repeated measurements of the same individuals, Hou and colleagues extend the cosinor “to a mixed-effect model framework, allowing for random effects” (Hou et al. 2021). JTK_CYCLE (Hughes et al. 2010) and RAIN (Thaben & Westermark 2014) are designed for detecting rhythms across many genes.

Worked example: a clock gene in wild-type and knockout mice

Constructed data: expression of a clock gene (log₂ scale) every 4 h from ZT0 to ZT44 (two cycles), in 3 different mice per time point and genotype (one sample per mouse). By construction the wild-type (WT) rhythm has amplitude 1.0 and peaks at ZT10; the knockout (KO) has amplitude 0.5 and peaks at ZT12; both have mesor 5.0. The random seed was the first one tried and was not chosen for its result. Stored code and data: rhythm.csv, make_data.py, analysis.py, analysis.R and the analysis output; analysis.py and analysis.R print the same output.

Statistic	WT	KO
Mesor	5.02	4.93
Amplitude (95% CI)	0.853 (0.674 to 1.033)	0.530 (0.390 to 0.671)
Peak time (95% CI)	ZT 10.1 (9.3 to 10.9)	ZT 13.0 (12.0 to 14.0)
Zero-amplitude test	F(2, 33) = 46.7, p = 2.4 × 10⁻¹⁰	F(2, 33) = 29.5, p = 4.4 × 10⁻⁸

Comparison	Result
Rhythm differs between genotypes (genotype × cosine and genotype × sine terms)	F(2, 66) = 13.58, p = 1.1 × 10⁻⁵
Amplitude, KO − WT	−0.323 (95% CI −0.547 to −0.099)
Peak time, KO − WT	2.8 h later (95% CI 1.5 to 4.2 h)
Separate t-tests at each of 12 time points	p < 0.05 at ZT 8, 12, 16, 28 and 36
The same, Bonferroni-corrected (p < 0.05/12)	ZT 16 only

Both genotypes are rhythmic; the knockout’s rhythm is smaller and peaks later. The estimates (amplitude difference −0.323, peak 2.8 h later) differ from the constructed truth (−0.5 and 2 h), but both true values lie inside their intervals. The joint F-test is the test of whether the rhythm differs; the amplitude and peak-time differences describe how, as estimates with intervals rather than as two further tests. The uncorrected time-point tests flag five scattered times, and the Bonferroni-corrected tests flag one; neither set describes either change. In 10,000 simulated studies of this design (simulation.py, the simulation output):

Scenario	Cosinor F-test for a difference, p < 0.05	At least one of 12 time-point t-tests, p < 0.05	The same, Bonferroni-corrected
Same rhythm in both genotypes	4.8% (false positives)	46.1% (false positives)	4.8% (false positives)
KO amplitude halved (1.0 vs 0.5)	91.9% (power)	82.9%	16.8% (power)
KO peak 3 h later (not the 2 h of the example)	99.9% (power)	97.2%	34.1% (power)

Uncorrected, the time-point tests flagged something in almost half the studies with no difference at all. A Bonferroni correction brings false positives back to 5% but then detects a halved amplitude in only 16.8% of studies, against 91.9% for the one model. Neither version says whether a change is in amplitude or timing.

The code expects one row per sample with columns genotype, zt, mouse and expr; rename your columns and group labels to match.

# Condensed, using statsmodels; analysis.py (NumPy and SciPy only) also gives the confidence intervals (delta method)
import numpy as np, pandas as pd
import statsmodels.formula.api as smf

d = pd.read_csv("rhythm.csv")                        # one row per mouse: genotype, zt (hours), mouse, expr
d["c"], d["s"] = np.cos(2 * np.pi * d.zt / 24), np.sin(2 * np.pi * d.zt / 24)
fit = smf.ols("expr ~ c + s", d[d.genotype == "WT"]).fit()
b1, b2 = fit.params["c"], fit.params["s"]
print(fit.params["Intercept"], np.hypot(b1, b2), (np.arctan2(b2, b1) * 24 / (2 * np.pi)) % 24)   # mesor, amplitude, peak (ZT)
print(fit.compare_f_test(smf.ols("expr ~ 1", d[d.genotype == "WT"]).fit()))   # zero-amplitude test: F, p, df
full = smf.ols("expr ~ (c + s) * genotype", d).fit(); red = smf.ols("expr ~ c + s + genotype", d).fit()
print(full.compare_f_test(red))                      # F, p, df: does the rhythm differ between genotypes?

d <- read.csv("rhythm.csv"); d$c <- cos(2 * pi * d$zt / 24); d$s <- sin(2 * pi * d$zt / 24)
fit <- lm(expr ~ c + s, data = subset(d, genotype == "WT")); b <- coef(fit)
c(mesor = b[[1]], amplitude = sqrt(b[["c"]]^2 + b[["s"]]^2), peak_ZT = (atan2(b[["s"]], b[["c"]]) * 24 / (2 * pi)) %% 24)
anova(lm(expr ~ 1, data = subset(d, genotype == "WT")), fit)                            # zero-amplitude test
anova(lm(expr ~ c + s + genotype, data = d), lm(expr ~ (c + s) * genotype, data = d))   # does the rhythm differ?
# analysis.R adds the confidence intervals and the amplitude and peak-time differences

Before you analyse

1. Sample across two full cycles, every 2 to 4 hours. Hughes et al. (2017) recommend 12 time points per cycle (every 2 h) for genome-scale studies; at least 6 per cycle (every 4 h), as in the example, is the older minimum and gives less precise peak times. Use the same times for every group, and record time as ZT (or circadian time in constant darkness).
2. Know whether samples are independent. Different mice at each time point (tissue sampling) fit the ordinary cosinor; repeated measurements of the same animals need a mixed model with animal as a random effect (Hou et al. 2021; repeated measures).
3. Check that a cosine fits. Plot the data with the fitted curve, and the residuals against fitted values; least-squares tests assume roughly normal residuals with similar spread, so analyse expression on a log scale, and use a generalised cosinor model for counts or proportions. Sharply peaked or two-peaked patterns need a second harmonic (extra 12-hour cosine and sine terms) or a different method.
4. Interpret amplitude and peak time only when the rhythm is clear. With a non-significant zero-amplitude test, the peak time is poorly defined, and delta-method intervals are unreliable when the amplitude is near zero. If a peak lies near ZT0, read its interval around the clock (ZT 23 to 1 is a 2-hour interval).
5. Use the period you assume consistently. The cosinor fixes the period (here 24 h); estimating the free-running period in constant conditions is a different analysis.

Check what an AI assistant ran

1. The assistant ran a t-test or post hoc test at every time point. Ask for one cosinor model with group terms.
2. The assistant reported the time of the highest raw mean as the peak. Ask for the fitted peak time with its interval.
3. The assistant treated repeated measurements of the same animals as independent. Ask for a mixed model.
4. The assistant reported an amplitude without a test of whether a rhythm exists. Ask for the zero-amplitude test.

Report it

Rhythmicity was assessed with a single-component cosinor model (24-h period) fitted by least squares, with a zero-amplitude F-test; amplitude and peak-time confidence intervals were obtained by the delta method. Genotypes were compared in one model including genotype × cosine and genotype × sine terms (F-test for the two terms; R [4.3.3]). Samples were collected every [4] h for [48] h, with [3] mice per time point and genotype.

Both genotypes showed a 24-h rhythm (WT: F(2, [33]) = [46.7]; KO: F(2, [33]) = [29.5]; p < 0.001 for both). The rhythm differed between genotypes (F(2, [66]) = [13.6], p < 0.001): amplitude, knockout minus wild type, [−0.32] (95% CI [−0.55] to [−0.10]); peak time [2.8] h later in the knockout (95% CI [1.5] to [4.2] h).

Replace the bracketed values with your own.

References

• Cornelissen G (2014). Cosinor-based rhythmometry. Theoretical Biology and Medical Modelling 11:16. doi:10.1186/1742-4682-11-16
• Parsons R, Parsons R, Garner N, Oster H, Rawashdeh O (2020). CircaCompare: a method to estimate and statistically support differences in mesor, amplitude and phase, between circadian rhythms. Bioinformatics 36:1208–1212. doi:10.1093/bioinformatics/btz730
• Hughes ME, Abruzzi KC, Allada R, et al. (2017). Guidelines for genome-scale analysis of biological rhythms. Journal of Biological Rhythms 32:380–393. doi:10.1177/0748730417728663
• Hou R, Tomalin LE, Suárez-Fariñas M (2021). cosinoRmixedeffects: an R package for mixed-effects cosinor models. BMC Bioinformatics 22:553. doi:10.1186/s12859-021-04463-3
• Hughes ME, Hogenesch JB, Kornacker K (2010). JTK_CYCLE: an efficient nonparametric algorithm for detecting rhythmic components in genome-scale data sets. Journal of Biological Rhythms 25:372–380. doi:10.1177/0748730410379711
• Thaben PF, Westermark PO (2014). Detecting rhythms in time series with RAIN. Journal of Biological Rhythms 29:391–400. doi:10.1177/0748730414553029