Modeling Physical Activity Change as Smooth Transformations: Temporal and Amplitude Patterns Associated with Physical Function in Older Women

The OPACH study was a Women’s Health Initiative (WHI, initiated in 1992 [40]) ancillary study that enrolled 7,048 ambulatory, community-dwelling women aged 63–99 years between 2012 and 2014 [41, 4, 7]. Among surviving OPACH participants, 2,356 consented to participate in the accelerometer sub-study of the WHISH trial (initially funded from 2015-2020) [42, 43, 44], with accelerometer deployments scheduled at 6, 18, and 36 months post-randomization. The focus of WHISH was to test whether increasing PA reduced the risk of primary CVD events in older women. Together with the OPACH baseline assessment, this design afforded up to four waves of accelerometer-measured PA data per participant. All participants provided written informed consent, and the studies were approved by the Institutional Review Board at the Fred Hutchinson Cancer Center [41, 43].

The present analyses focus on the first three accelerometer assessment waves: OPACH baseline (original $n_{0}$ = 6,489), WHISH wave 1 at 6 months post-randomization (W1; $n_{1}$ = 1,934), and WHISH wave 2 at 18 months post-randomization (W2; $n_{2}$ = 1,677). The average interval between baseline to W1 was 3 years and the average between W1 to W2 was 1 year. We refer to this analytic sample as the OPACH–WHISH sub-cohort. Data from the WHISH wave 3 at 36-month (W3; $n_{3}$ = 1,398) were excluded due to substantial loss to follow-up, to minimize potential bias in FPCA estimation and subsequent association analyses. At the time of the present analysis, WHISH randomization assignments (usual activity versus PA intervention) for OPACH–WHISH sub-cohort members had not yet been made available to the study team; implications are addressed in the Concluding Remarks (5).

Participants were provided with an ActiGraph GT3X+ (Pensacola, Florida) tri-axial accelerometer and instructed to wear it over the right hip, 24 hours per day for 7 days at each assessment wave, only to remove it for bathing and other water-based activities [41, 42]. The device captures acceleration along three orthogonal axes: vertical, horizontal, and perpendicular (VA, HA, PPA) [45, 46]. 15-second epoch count files were generated via ActiLife software (v.6) [41, 42]. Counts from each axis were summed over four consecutive 15-second intervals to obtain minute-level counts, then combined across axes to compute the vector magnitude (VM) as $\sqrt{\mbox{VA}^{2}+\mbox{HA}^{2}+\mbox{PPA}^{2}}$ [47]. Non-wear windows were identified using the Choi algorithm [48, 49].

Inclusion criteria and valid-day definitions. For each wave, minute-level counts from 6:00 a.m. to midnight (18 hours) were retained, as activity outside this window predominantly reflected sleep or in-bed time [32]. Days with $\geq$4 hours of non-wear between 6:00 a.m. and midnight were excluded (corresponding to the $90^{\mathrm{th}}$ percentile of daily non-wear), and a valid day required $\geq$14 hours of wear time within this window. Participants were required to have at least four valid days per visit to be included [46, 50, 5]. These criteria yielded $N_{1}=1{,}646$ participants with valid data at both baseline and W1, and $N_{2}=1{,}362$ with valid data at both W1 and W2 (84 of whom lacked qualified baseline data). The final accelerometer analytic sample included $N^{+}=1{,}730$ $(1,646+84)$ participants, each with at least two consecutive non-missing visits (mean valid days per visit: 6.3; median: 7; range: 4–10).

The retained VM for each participant were averaged across days within each visit, yielding a single diurnal activity profile per participant per visit. VM counts were aligned to chronological time ($t=1,\dots,1080$, representing minutes from 6 a.m. to midnight), and a cubic smoothing spline (25 degrees of freedom) was applied at each participant-visit to reduce measurement noise. Figure 1 illustrates one participant’s PA profile from daily minute-level measurements (1(a)) to the smoothed curve (1(b)) at baseline, with smoothed curves at W1 (1(c)) and W2 (1(d)).

The grand mean and standard deviation (SD) of the smoothed PA counts were computed pooling across all participants and visits. Each participant–visit curve was then centered by this grand mean and scaled by $4\times\text{SD}$ so that the majority of scaled curves lie approximately within $[-1,1]$. The time domain $[1,1080]$ was likewise rescaled to $[-1,1]$, balancing the penalty weights on deformation energies along the amplitude and time axes [32]. For each participant at each visit, the signed Area Under the scaled PA Curve (AUC) was defined as $\text{net-AUC}\;=\;\int_{\{t:\,y(t)>0\}}y(t)\,dt\;-\;\int_{\{t:\,y(t)<0\}}\lvert y(t)\rvert\,dt$, where $y(t)$ denotes the scaled PA curve. Although this reduces to the ordinary integral $\int y(t)\,dt$ that can be estimated by trapezoid rule, the decomposition emphasizes the directional interpretation. The between-visit difference, $\Delta\text{net-AUC}^{(\eta)}$, then served as a signed summary of net PA change, where $\eta\in\{0,1\}$ indexes the transition period ($\eta=0$: baseline to W1; $\eta=1$: W1 to W2). Positive values indicate an increase in overall PA; negative values indicate a decrease.

Physical function (PF) reflects the capacity to perform daily activities requiring locomotion, balance, and muscular endurance, and its decline is a well-established precursor to disability, loss of independence, and mortality in older adults [51, 52]. Importantly, low PF is a risk factor for CVD [53, 54] and predicts incident CVD events after adjustment for traditional risk factors [55, 56]. There is also compelling evidence that PA is central to preserving PF and mobility in later life [57, 58, 59, 60, 61, 62, 63]. Given its dual role as a key pathway through which PA influences aging-related outcomes and as an independent predictor of incident CVD, PF served as the primary outcome of this study.

PF was assessed using RAND-36 Physical Functioning subscale, which quantifies health-related limitations in daily physical activities on a 0–100 scale, with higher values indicating better functioning [64]. Because participants’ RAND-36 assessments and accelerometer wear wave did not always occur simultaneously, the first day of accelerometer wear at each of W1 and W2 was matched to the nearest RAND-36 assessment within a $\pm$190-day window. This window was chosen to minimize the time discrepancy between two different instruments while maximizing participant retention.

Baseline covariates were selected based on prior literatures and clinical relevance [7, 5]. With the exception of race/ethnicity and education, which were obtained at the WHI baseline, all others were collected at the OPACH baseline. Questionnaires assessed age, race/ethnicity (Black, White, Hispanic/Latina), education (high school equivalent or less, some college, college graduate or higher), alcohol consumption (nondrinker, $<$ 1 drink/week, $\geq$1 drink/week, unknown), current smoking status (no, yes), and self-rated general health (excellent/very good, good, poor/very poor). Moderate-to-vigorous PA (MVPA) was defined using the OPACH Calibration Study VM count cut-point ($\geq$519 counts per 15 seconds) [20]. Baseline average daily MVPA (minutes/day) was computed over awake wear time by summing MVPA across valid days and averaging over total number of valid days.

Multimorbidity was categorized as 0, 1, 2, or $\geq$ 3 of 12 prevalent conditions: history of CVD, cancer, diabetes, hip fracture, osteoarthritis, depression, chronic obstructive pulmonary disease, cognitive impairment, self-reported moderate/severe vision or hearing loss, $\geq$2 falls in past 12 months, and urinary incontinence [65]. Study staff measured systolic and diastolic blood pressure (SBP and DBP; mmHg) by auscultation using an aneroid sphygmomanometer, with participants seated following 5 minutes of quiet rest; 2 measurements were averaged. Height and weight were assessed by tape measure and bathroom scale to calculate BMI (kg/m²). Fasting (12-hour) blood samples were collected close in time ($\leq$6 months) to the accelerometer wear interval. Serum glucose, insulin, total cholesterol, high-density lipo-protein cholesterol (HDL-C), and high-sensitivity C-reactive protein (hsCRP) were measured using standardized Clinical Laboratory Improvement Act–approved methods at the University of Minnesota Fairview Advanced Research and Diagnostics Laboratory [66]. The Reynolds Risk Score, a summary CVD risk score representing the predicted 10-Year probability of having a clinical CVD event, was calculated using age, SBP, hsCRP, total and HDL cholesterol, diabetes, smoking, and family history of myocardial infarction [67].

We modeled longitudinal change in diurnal PA patterns using the Riemannian diffeomorphic framework of Zou et al. [32]. Participant’s daily PA pattern was represented as a smooth curve at each visit, the longitudinal change from two consecutive visits (e.g, baseline – W1, W1 – W2) were modeled as a Riemann Manifold diffeomorphism, a smooth transformation that deforms the source curve into the target curve.

Formally, let $C_{i}^{(q)}$ denote the $i^{th}$ participant’s smoothed PA curve at visit $q$. The change from visit $q{-}1$ to visit $q$ is modeled as

where $\gamma(\boldsymbol{v}_{i}^{(q)},\cdot)$ is a subject-specific diffeomorphism that maps the source curve $C_{i}^{(q-1)}$ to the target curve $C_{i}^{(q)}$, and $\circ$ denotes function composition.

The diffeomorphism is governed by a time-varying vector field $\boldsymbol{v}_{i}^{(q)}$ embedded in a reproducing kernel Hilbert space (RKHS, [68]) equipped with a Gaussian isotropic kernel [35, 32]. The deformation evolves over an auxiliary parameter $\tau\in[0,1]$, where $\tau=0$ corresponds to the source curve and $\tau=1$ to the target. Intuitively, the vector field $\boldsymbol{v}_{\tau}$ can be viewed as a series of forces that continuously deform the PA curve from one visit to the next at each infinitesimal step $\tau$. Crucially, these forces operate simultaneously in both the amplitude ($y$) and phase ($x$) domains, enabling the model to capture changes in PA magnitude (e.g., an overall increase or decrease in activity) and temporal/phase shifts (e.g., redistribution of PA bouts to earlier or later hours) within a single unified framework.

Subject to the constraint on deformation energy, $\int_{0}^{1}\|\boldsymbol{v}_{\tau}\|^{2}\,d\tau$, the deformation in Eq (1) is unique and fully determined by the initial momenta: the values of the vector field at $\tau=0$ evaluated at a set of preselected control points on the curve. In our analysis, the control points are selected as the $1{,}080$ minutes from 6:00 a.m. to midnight. At each minute $p$, the initial momentum is a two-dimensional vector,

where the $x$-component encodes the local temporal shift and the $y$-component encodes the local amplitude change; $\eta\in\{0,1\}$ corresponds to transition periods (baseline – W1 and W1 – W2). The optimal initial momenta are calculated by minimizing a penalized objective that balances the similarity of the deformed source curve to the observed target against total deformation energy, solved via Hamiltonian optimization using the fshapesTK package [35]. For participant $i$ and period $\eta$, we summarize the overall amount of deformation by the initial deformation energy, $d_{i}^{(\eta)}=\sum_{p=1}^{1080}\left[\bigl(m^{(\eta)}_{i,p,x}\bigr)^{2}+\bigl(m^{(\eta)}_{i,p,y}\bigr)^{2}\right]$. Details of initial momentum estimation, deformation energy and optimization can be found in previous publications [35, 32]. As an example, Figure 2 showed the estimated deformations and initial momenta of one participant from the OPACH-WHISH sub-cohort.

To extract the dominant modes of longitudinal PA change and to preserve the distinction between these two clinically interpretable $x$ (Temporal) and $y$ (Amplitude) domains, we applied multivariate functional principal component analysis (MFPCA) using the approach of Happ and Greven [38].

In this framework, univariate eigenfunctions and scores were first estimated separately on each domain, and the multivariate eigenfunctions and scores were obtained as linear combinations of their univariate counterparts with a weight factor. The estimation steps were described in more details in Appendix (A) and [38]. MFPCA was applied separately to each transition period with identical parameter settings.

An alternative, simpler strategy is to concatenate the $x$- and $y$-components of the initial momenta into a single vector function and apply standard univariate FPCA (UFPCA) [69]. If the measurement on both domains had identical grids (time-points), comparable variance, and the principal modes were not highly correlated across domains, concatenated UFPCA can approximate MFPCA reasonably. Otherwise, naive concatenation risks producing blended components in which temporal-shift modes are absorbed into amplitude modes (or vice versa), obscuring domain-specific interpretability. We therefore report results from both approaches and quantify their differences.

Three sets of predictors were derived in earlier sections to characterize longitudinal changes in diurnal PA patterns: (i) MFPCA principal component scores (Section 3.2), which capture the dominant modes of change in both temporal shift and amplitude; (ii) $\Delta\text{net-AUC}^{(\eta)}$ (Section 2.2), which summarizes net directional change (positive versus negative) in PA without information on temporal shift, and (iii) initial deformation energy $d_{i}^{(\eta)}$ (Section 3.1), which summarizes overall amount of PA change without directionality in timing and amplitude.

Because PC1 and $\Delta$net-AUC were highly collinear in both periods (Pearson $r=0.9$), they were not included in the same model. In contrast, deformation energy showed weak correlations with PC1 ($|r|\leq 0.15$) and with $\Delta$net-AUC ($|r|\leq 0.30$) in both periods, so deformation energy was modeled jointly with PC1 and, separately, with $\Delta$net-AUC. Two LMMs were carried out with participant-specific random effects to account for repeated PF outcome at W1 and W2: Model 1 included the stacked deformation-based change features (PC scores) and deformation energy, together with baseline PF scores and other covariates; Model 2 replaced the PC scores with the scalar change summary, $\Delta$net-AUC, while retaining deformation energy, baseline PF and covariates. Inclusion of interaction terms between period and the predictors were determined using likelihood ratio tests (LRTs) comparing nested models. To reduce overfitting and improve interpretability, we applied LASSO [70] to select covariates and PCs in the LMM via R package glmmLasso[71]. The LASSO penalty parameter was determined by Bayesian Information Criterion (BIC), and the significance threshold for two LMMs was adjusted to $\alpha=0.025$ by Bonferroni method.

All statistical analyses were conducted in the R statistical programming environment (version 4.3.1) [72], except for the initial momenta estimation and optimization steps, which were implemented using the MATLAB package fshapesTK [35] and executed on the National Research Platform Nautilus Kubernetes Cluster [73]. R package nlme for LMM is available from [74] and package MFPCA is available from [75].

At both periods, baseline–W1 and W1–W2, 15 PCs are required in MFPCA to explain at least 90% of the variability, with PC1 accounting for the largest share among these components (22.4% and 20.8%, respectively), followed by PC2 (10.9% and 11.5%). Detailed variance explained by each PC in both periods is reported in Appendix (B).

Figure 3 presents the MFPCA-estimated mean functions (top row) and corresponding deformations (bottom row) for the $x$ (temporal shift) and $y$ (magnitude shift) domains. In both periods, the mean $x$-domain function cross zero near 10:20 a.m., transitioning from rightward (positive) to leftward (negative) temporal shifts (3(a)). The mean $y$-domain function remains negative throughout, indicating persistent downward vertical shifts, with the deepest drop occurring at 7:52 p.m. in baseline–W1 (3(b)).

The bottom panels illustrate how the estimated initial momenta deform one visit’s smoothed PA curves toward the follow-up visit’s curves. In the baseline–W1 panel (3(c)), arrows point rightward and downward before 10:20 a.m. and leftward and downward thereafter, dragging the observed baseline mean (blue) toward W1 (green); the momentum-derived estimate mean PA curve (dashed red) closely follows the observed W1 mean. The W1–W2 panel (3(d)) exhibits a similar deformation pattern but with attenuated temporal and vertical shifts, consistent with the shorter average inter-visit interval (one year for W1–W2 versus three years for baseline–W1), and may also reflect smaller change in PA patterns at older ages.

Figure 4 presents the PC1 eigenfunction (top row) and its corresponding deformation (bottom row). Temporal shifts of PA under PC1 are similar across both periods (4(a)): predominantly leftward before 10:00 a.m. (negative eigenfunction $x$-coordinates) and rightward afterward (positive $x$-coordinates). Magnitude changes in PA are also consistent across periods and primarily upward (positive $y$-coordinates) (4(b)). The $y$-coordinate change is more pronounced than the $x$-coordinate change, indicating that PC1 captures more variability in PA magnitude than in timing. As with any PC, this mode is bidirectional: the plotted eigenfunction represents one direction of PC1, so a positive PC1 score corresponds to change in the displayed direction and a negative score to the opposite. These changes are visualized more directly in the bottom deformation panels (4(c), 4(d)): blue curves (identical to those in Figure 3(c), 3(d)) represent the observed sample mean curves at the earlier visit of each period (baseline and W1, respectively), and dashed red curves show the estimated curves deformed from the blue curves following PC1’s initial momenta. Because PC1 explains the largest portion of variability after removal of the mean function, it captures the residual primary mode of diurnal PA change, with positive scores indicating increase and negative indicating decrease in PA, especially between 10 a.m. and 7 p.m.

Deformation panels of PC2 are shown in Figure (5). Both temporal shifts and magnitude changes are consistent across the two periods. PC2 characterizes a leftward shift and localized increase in morning activity for positive scores (or a decrease for negative scores), illustrated by the up-and-leftward-pointing arrows from 6 a.m. to around 11 a.m. (5(a), 5(b)), followed by a slight downward shift thereafter.

Additional eigenfunctions and deformation figures are provided in Appendix (C), revealing distinct patterns of PA change during morning, afternoon, and evening phases. Beyond the top four PCs, each additional component explains only a small proportion of variability and PA changes associated with those PCs are minimal.

Importantly, in this OPACH–WHISH subcohort, the eigenfunctions and associated deformations of PA patterns for the top four PCs show similar patterns across the two periods. This cross-period similarity makes it reasonable to treat, for example, PC1 from baseline–W1 and PC1 from W1–W2 as representing the same underlying mode of PA change. Accordingly, the corresponding PC scores can be stacked across the two periods and included as a common predictor in the subsequent LMM association analyses.

Figure 6 shows comparisons between the MFPCA and concatenated UFPCA mean functions (top) and first PC eigenfunctions (bottom). The two approaches yield generally similar mean structures and PC1 patterns across both domains and both periods. However, the solid curves (UFPCA) exhibit a “smooth-over effect” at the transition between the end of the $x$-domain and the beginning of the $y$-domain (solid lines can be connected between left and right panels) , a feature that is absent in the dashed curves (MFPCA). This occurs because concatenation forced the two domains to join smoothly, thereby imposing continuity at the boundary and inducing contamination into the data structure between two domains.

Appendix (D) provides additional comparisons for PC2-PC5. Overall, the two methods produced consistent patterns, particularly for the four leading PCs. This consistency is likely a benefit from the centering and scaling applied to the $x$- and $y$-domain functions prior to decomposition, which place the two domains on comparable scales and reduced the chance that concatenation would artificially overweight one domain or obscure domain-specific structure.

For association analyses, two LMMs were described in Section (3.3) and LRT revealed significant deformation-energy-by-period interaction in both models ($p<0.01$). In Model 1, LASSO removed PC2, PC3, PC4, baseline smoking status, and DBP; in Model 2, baseline smoking status and DBP were removed as well.

Because PF and accelerometer assessments were not collected concurrently, PF was missing for 260 participants at W1 and 338 at W2 even after allowing a $\pm$190-day matching window. Several baseline covariates also had substantial missingness (e.g., 423 missing Reynolds Risk Score and 329 missing cholesterol/glucose/blood pressure measures). Furthermore, 106 participants with previous CVD diagnosis at baseline were excluded. After restricting to participants with complete data on PF and all retained covariates, the analytic sample size for both Model 1 and Model 2 was $N=1{,}157$. Significant terms of fixed-effect estimates are reported in Tables 2 and 3.

In Model 1, PC1 is positively associated with PF score ($\beta$ = 60.1, s.e = 14.51, $p<0.0001$) after adjusting for Lasso-selected covariates including baseline PF score, where increase in PC1 corresponding to an upward deformation as shown in Figures 4c, 4d. That is, increasing the PC1 score from 25^th to 75^th percentile, which corresponds to a 0.02-unit increase, is associated with an increase of 1.3-units in the PF score (s.e = 0.33).

Figure 7(a) visualizes the interaction over the observed deformation-energy range [0, 0.04]. The deformation-energy-per-period interaction is significant ($\beta$ = 387.4, s.e = 129.84, $p=0.003$), indicating that an increase of 0.04 in deformation energy has a 15.5 greater effect on PF in W1–W2 than in baseline–W1. This pattern suggests that PF may have been more sensitive to overall PA-pattern change in the later period, possibly because participants were older during W1–W2 and may have accumulated more prolonged sedentary time, making PF more responsive even to modest perturbations in PA patterns.

Among baseline covariates, age and BMI are negatively associated with outcome: one year increase in age is associated with 0.5 unit decrease in PF score (s.e = 0.10, $p<0.0001$) whereas 1 kg/m² increase in BMI is associated with 0.3 unit decrease in PF score (s.e = 0.09, $p=0.0004$). Self-rated good (vs excellent) health is associated with lower PF (2.4 units, s.e = 0.99, $p=0.02$). Higher baseline MVPA is associated with higher PF: one hour increase in baseline daily MVPA time is associated with 1.8 unit higher in PF score (s.e = 0.6, $p=0.02$) after fixing other predictors and covariates.

Results are similar in Model 2 (Table 3). $\Delta$net-AUC, used as a scalar summary of total activity change between consecutive visits, is positively associated with PF: a one-unit increase in $\Delta$net-AUC is associated with a 0.01 unit increase in PF score (s.e = 0.003, $p=0.0001$). Note that $\Delta$net-AUC was derived based on scaled net-AUC values, a one-unit change in $\Delta$net-AUC is small relative to scale of PF (0 – 100). The deformation-energy-per-period interaction remains significant ($p=0.005$), indicating an increase of 0.04 over observed deformation energy range has a 14.6 greater effect on PF in W1–W2 than in baseline–W1 (Figure 7(b)). Associations between significant baseline characteristics and PF in Model 2 are similar to those from Model 1 (Table 3).

All participants provided written informed consent to the WHI study, which was approved by the Institutional Review Board at the Fred Hutchinson Cancer Center. Participants to the WHISH study were WHI participants who were consented using a passive “opt-out” consent, approved by the Institutional Review Board at the Fred Hutchinson Cancer Center. University of California at San Diego exempted Institutional Review Board approval for the use of existing data for the present study.

Not applicable.

The authors do not have permission to share data.

The authors declare no potential conflict of interests.

JZ’s, LD’s, and AL’s work was supported by NIH/NHLBI R01HL166802. RZ’s and LN’s work was supported by NIH/NHLBI R01HL166802 and NIH/NHLBI 1R01HL168535. CD’s work was supported by NIH/NHLBI R01HL130483. RZ’s and LN’s work was partially supported by NIH/NHLBI R01HL130483.

RZ and JZ drafted the manuscript. JZ conceived the statistical design and its application to the OPACH-WHISH study with support from AL. RZ and JZ led the statistical analyses, with support from LN, CD, and CK. SN curated the OPACH baseline clinical data and contributed to their interpretation. YW supported implementation in the cluster computing environment for high-dimensional data. LD assisted with accelerometer data acquisition and management. ML, PR, RC, MS and SH contributed to the critical review and interpretation of the data. All authors contributed to review and edits of the manuscript.