feat(oak): rest-activity rhythm metrics (IS, IV, RA, cosinor) - #323
ceyhunolcan wants to merge 9 commits into
Conversation
|
(gotta stick this somewhere) |
|
Thanks, that's kind of you! For attribution: Ceyhun Olcan Affiliation (if the file takes one): Dartmouth College Happy with whatever format fits your contributors file, name alone is fine if ORCID/affiliation don't fit the schema. |
|
Pushed a follow-up: added tests for the |
|
(ok, I haven't gotten to a point where I can review this one yet) |
|
I did not realize this was an entirely new analysis, very cool. Since there are no conflicts we could be ready to go, but with my work on the ruff.... branch I wanted to go through this in full, my process is to read the code, do my formatting changes, then I add that formatting commit to the git blame ignorer. I'm currently doing this process on a new branch Review for inclusion - I probably cannot just review and OK this, I will need to task someon on our end is qualified to do that type of review. @hydawo This is the work I mentioned. Question: very roughly and non-rigorously, what is the performance like on this? |
|
Yo, I just pushed 8e2e124 on Tests pass. Why don't you do a sanity check? Then we can work out next step for some code review from me while we are working out the more academic code review. I'll propose you |
|
update: it is this commit a225e2c |
4875240 to
b7ead43
Compare
|
Done! Reset to match a225e2c and pushed. On the sanity check: the mypy failures were config/typing, not logic. After the mypy.ini → pyproject migration On performance: per participant, cost is dominated by reading and epoching the raw accelerometer CSVs, not the metrics. IS/IV/RA are single vectorized passes over the per-epoch vector and the cosinor is a 3-column least-squares, all linear in epochs, so a multi-week recording at 1-min epochs is well under a second of compute once loaded. I/O scales with study size; the math doesn't. |
|
mypy - mypy passes on the ruff-check-reflow-and-typing branch without modification, it did not involve numba. Numba should not be added to the mypy ignore module. Pandas should have been removed from that list and it is a bug that it was still present to begin with. I have now removed it in a commit on the ruff-check-reflow-and-typing branch. I think if you revert your commit there is no merge conflict. I'm also going to need a statement from you on use of ai tooling here, we are developing an ai tool policy but don't have one yet. |
|
Dropped the numba overrides so pyproject matches your ruff-check-reflow-and-typing config (no numba, no pandas; pandas type-checked strictly). The DatetimeIndex fix in _align_to_day_grid stays, since that's what makes rhythms.py pass under strict pandas typing. Should merge cleanly now. On AI tooling: the substance of this work is my own. The rest-activity rhythm methods (IS, IV, RA, cosinor) come from a library I wrote and validated against real NHANES data before this PR and that's my prior research, not anything generated here. The engineering calls were mine as well: what to fix, what to leave alone, how to structure the tests, and how to respond to review. I do use AI (Claude) in my workflow, mainly for drafting and thinking through options, but I direct it and verify everything it helps with for this PR that meant reproducing each issue myself, writing and running the tests, running the full ruff/mypy/pytest gauntlet locally, and confirming behaviour before pushing I don't submit anything I can't explain or haven't checked. So the methods, the judgment, and the verification are mine, and I stand behind all of the code. I think transparency here is the right default, so I'm glad you're putting a policy together and happy to follow whatever the lab lands on. |
|
(been busy with other issues, putting in some time over here. I'll be back at work Tuesday, so no need for an instant turnaround.) I had a colleague take a look at this pull request (academic perspective), I'm going to detail that and then respond to your comment and do further re-review. Requests, Recommendations
|
|
|
Code style request
|
|
Thanks and no rush understood. On references: the methods come from actrhythm, a small library I wrote implementing these non-parametric rest-activity metrics (IS, IV, RA with L5/M10) and a single-component cosinor: github.com/ceyhunolcan/actrhythm. Method references are Van Someren (1999) and Gonçalves (2014) for the nonparametric metrics and Cornelissen for the cosinor the ones cited in rhythms.md. I developed and applied these in a preprint of mine, a cross-sectional NHANES 2013-2014 study of olfactory dysfunction and 24-hour activity-rhythm fragmentation (n=2,327), on Research Square: https://doi.org/10.21203/rs.3.rs-9830931/v1 (analysis code at github.com/ceyhunolcan/od-activity-rhythm, Zenodo 10.5281/zenodo.20132927). In actrhythm the implementations are pinned to analytically-derived values (IV of an alternating series = 4.0, IS of identical days = 1.0, RA on a known profile = 9/11), the same ground-truth style as the tests in this PR. Glad to add any further refs into rhythms.md. On your colleague's review: thank you, these are careful and well-taken points. A couple (the IV missing-data normalization and coverage-dependence of the rhythm metrics) I'd rather handle properly than patch quickly, so I'll work through them, a minimum daily-coverage criterion before counting valid days, IV normalized over valid consecutive pairs, the average-day-profile fill that can bias L5/RA, and the acrophase-to-local-clock mapping, along with the |
|
Pushed the cosinor acrophase precondition note, the epoch_seconds/epochs_per_day document clarification, and the IV normalization correction (now dividing over valid consecutive pairs instead of the raw epoch count, with a gappy-data regression test). Since the two coverage-related concerns (#1 and #5) are essentially one question concerning the amount of coverage we need before reporting a metric, I would like to validate the policy with you before implementing them. #5 (average-day profile): As your colleague pointed out, _average_day_profile currently uses the global profile mean to fill completely empty clock bins, which biases L5/RA because a missing nighttime bin receives a daytime-inclusive average and appears overly busy. Under-covered bins will remain as NaN after I remove that fill. How should windows that straddle missing bins be handled by L5/M10? If we choose strict: L5/M10 return NaN when there is no complete window, and a 5h/10h window is only valid if all of its bins are observed. Metrics that are straightforward and conservative are only reported when they are fully supported and if we choose tolerant: if coverage exceeds a criterion (such as ≥2/3 or ≥80%), a window mean is calculated over its observed bins; otherwise, NaN. My lean is strict, it's the cleaner story for a statistical default, and it composes with the valid-day gate below, so we only report metrics we can stand behind. The tradeoff is more NaNs on gappy records, which I'd argue is the honest outcome. Either way I'll make L5/M10 NaN-aware and document the rule. #1 (valid-day counting): related, and I saw that the module already had the necessary machinery. The MIN_DAY_COVERAGE = 0.5 constant is already used by _daily_frame to gate daily rows every day. Although no single day is sufficiently covered, six half-observed days count as three valid days and pass min_valid_days=3 because the participant-level valid_days count in run() is (all non-NaN epochs) / epochs_per_day. If each one of them passes the coverage threshold, I'll adjust it to just count one day. The policy challenge there is whether to adopt a more stringent wear threshold or maintain the current 0.5, which unites both gates on a single constant. Although I might suggest using ≥16 h/day (~0.67) in previous actigraphy work, 0.5 is more in line with what is already included in the module. If you don't have a preference, I'd be happy to use my leans (strict windows + uniform 0.5). I just wanted to confirm the thresholds with you first because these are methodological decisions that should be agreed upon before I lock them in. |
|
Quick housekeeping note: GitHub is now flagging a pyproject.toml conflict. It's just the mypy overrides block, your shapely/type-check cleanup on the base (dropping holidays/ratelimit/shapely/timezonefinder and enabling those checks) overlaps the pandas-override line I'd removed earlier, so it's a textual overlap rather than a real disagreement; your version is the fuller cleanup and already covers it. Happy to rebase onto the current base and take your pyproject.toml, but since your review process rebases this branch anyway, I'll leave it for you to fold in unless you'd rather I push a rebase. (Just let me know and letting you know about it ) |
|
You are welcome to rebase or merge at your preference and pace, this is 99% new code so excepting glitches like the .toml details it should be low on conflicts. (there's always I'll review a final file list before merging.) I keep a There are other changes coming in, I'll mention or make post-merge cleanups if something comes up. |
a6d9d96 to
222b0e1
Compare
|
(the target branch has been merged into develop, retargeting this pull to develop) |
… very minor edits.
…view) - IV divides successive-diff and deviation sums by valid-pair and valid-point counts, not raw epoch count, so missing epochs don't bias it (reduces to the Van Someren n/(n-1) form when complete); add gappy-data regression test - clarify epoch_seconds vs epochs_per_day (86_400/epoch_seconds) in docs - document cosinor acrophase midnight-alignment precondition; use 86_400
222b0e1 to
cd6dce9
Compare
|
Re-rebased onto develop now that the base merged in the retarget was clean, no conflicts (my changes are additive to oak plus the docs and blame-ignore entries, and pyproject already matches develop). Full oak suite (26), mypy, and ruff all pass on top of develop's current tip. Ready for your mirror to re-sync, and for the final file-list review whenever you get to it. |
JinjooThank you for answering my questions and for making these revisions. I have two follow-ups from my colleague, @jinjoo
Similarly, pyActigraphy already implements non-parametric RAR metrics: Is there a specific reason to hard-code these methods within Forest instead of using or wrapping these existing packages?
[And the rest is me again] There have been some surrounding code changes
Typing issues
General
Documentation
Naming Changes
|
|
Thanks both, these are useful comments. Daily coverage threshold: you're correct, and I don't have a reference for 0.5. It was already in the module as the gate for generating daily rows and I repurposed it, which was the wrong choice for a participant-level inclusion Existing packages: I did look at both, and my reasoning is close to Eli's point. Each non-parametric metric is a short, well-defined formula, so a small named function per metric is easier to read, to test against analytically known values, and to reason about than wrapping a package with a much larger API. It also avoids adding a dependency for what amounts to a few dozen lines of arithmetic. CosinorPy is a good package and does considerably more than this module needs (multi-component fits, diagnostics, plotting), and pyActigraphy brings its own data model and readers for device formats that Forest does not use. If the lab later wants multi-component cosinor or the wider pyActigraphy feature set, wrapping one of them would be the sensible route, and I will note that in the documentation. Oak reorganization: my module splits along the same lines as the split in the tree. Preprocessing includes the ENMO calculation, the timestamp index, and the day-grid alignment; analysis includes IS, IV, RA with L5/M10, the average-day profile, and the cosinor; the runner includes the participant loop, frequency handling, and CSV writing. So the separation applies to me directly. My instinct is to keep it as a self-contained module that follows the same internal separation rhythm metrics are a parallel analysis to step counting rather than another layer of it, they share no preprocessing (ENMO per interval is a different operation from the 10 Hz resampling in preprocess_bout), and keeping them together makes the module readable on its own. You are deciding how tree-level components are contained, though, so tell me which you prefer. If you do want the fold-in, I would rather do it as a follow-up than in this branch, so the file list you review stays a single new module instead of three modified ones. Naming and constants: moving to forest.constants. SECONDS_IN_DAY, renaming rhythms.py to circadian_rhythms.py, expanding rar to the full words, removing the leading underscores, and renaming _parse_time_utc to state that it constructs the index, while keeping the coverage constant in the file for now. I'm substituting "interval" for "epoch" throughout, so interval_seconds and intervals_per_day. One thing to flag: that rename also changes the n_valid_epochs output column, which is user-facing, so say if you would rather the column names stay as they are. Also happy to add a file header in the same format as the other oak modules if you want authorship recorded that way. Typing: in line with the pattern in bonsai and jasmine, cosinor will return a dataclass with typed fields instead of a dict. Since the row dict is the one that has to remain, rar_metrics will take the list and append its row inside. I'll add element types to the list annotations. On the os.path discrepancy: I used pathjoin for the joins but left two isdir calls bare. I appreciate you pointing that out, I should have caught it. fixing both. I'll do the rewriting myself. Removing the installation section and the "executive summary" heading, spelling out ENMO and the rest on first use, writing in longer sentences, and adding worked output examples along with guidance on how to interpret the numbers, including what unusual values indicate about the data and what the literature links them to. And you're right, I don't speak English as my first language, so that's a sensible guess rather than an unreasonable one. The tone check would be very appreciated. When something is off, just let me know and I'll take care of it. |
|
Also related to the column-name question the output files are also rar_summary.csv and <backend_id>_rar_daily.csv, so the same question applies there since "rar" is one of the acronyms you asked me to expand. I've left both the columns and the filenames alone pending your preference. |
|
@ceyhunolcan please take a look at my recent comments on the #279 issue about the requirements for LLM contributions. |
I just read the guidelines after seeing your note on #279. I'm finishing up the remaining points. When I push I'll include the AI usage statement in addition to the output documentation. |
2 participants
Adds
oak/rhythms.py: rest-activity rhythm summaries from accelerometer data interdaily stability (IS), intradaily variability (IV), relative amplitude (RA) with L5/M10, and a 24 h cosinor (MESOR, amplitude, acrophase, R²). Same raw accelerometer stream as the gait/step output, a different (24 h patterning) question. Activity is ENMO per epoch;run()mirrorsoak.base.run.Placement: a module inside oak rather than a new tree, since RAR is an accelerometer analysis. It writes its own output (recording-level
rar_summary.csv, plus per-participant<backend_id>_rar_daily.csvat a daily frequency) rather than extending_gait_daily.csv, because IS/IV are recording-level while gait output is per-day.Notes / open questions
min_amp;gravityparam handles m/s².Tests:
tests/oak/test_rhythms.pyground-truth recovery (clean rhythm recovers acrophase/MESOR/amplitude, IS≈1; white noise → IS at 1/n_days, IV≈2; identical days → IS=1) + end to-end. ruff/flake8/mypy clean; full suite green.Docs:
docs/source/rhythms.md(+ toctree).