Getting Started with menstrualcycleR and Phase-Aligned Cycle Time Scaling (PACTS)

1.5.1 Optional: imputing the next menses from a confirmed ovulation

The imputation above works in one direction: given an observed menses onset, it estimates the ovulation 15 days earlier. The reverse case also happens: a cycle has a biomarker-confirmed ovulation but the next menses onset was never recorded (the participant stopped surveying, or the closing period was missed). By default such a cycle cannot be scaled, because there is no anchor to close it.

Setting impute_next_menses = TRUE in pacts_scaling() opts in to closing that cycle by imputing a menses onset forward from the confirmed ovulation, placed at ovulation + next_menses_luteal_days (default 14 days, so ovulation + 14 = the last follicular day, i.e. the “LH+15” convention). Imputation is skipped whenever an observed menses already closes the cycle, at any distance, so a recorded period is never overridden – the next_menses_max_window argument from earlier versions no longer affects this (as of 0.1.7; see ?pacts_scaling), because bounding that check by a day count was itself a bug: a real closing menses just beyond the old window was invisible to it and got overwritten. Imputed onsets are flagged in a new menses_impute column, exactly as imputed ovulations are flagged in ovtoday_impute, so the two can be reported separately. This option is off by default, so existing analyses are unchanged; only the general rule is applied, and any study-specific gating (for example, not imputing across a treatment phase or a documented off-study break) is left to the caller.

4.1.1 1. Cycle Scaling Function: pacts_scaling()

This function adds standardized cycle time variables to your dataset.

• scaled_cycleday:
  A continuous cycle time variable centered on menses onset (menses == 1 → 0), ranging from –1 (start of luteal phase) to +1 (ovulation). This maps -1 and +1 as adjacent days. This is what the supplementary analyses of Nagpal et al. (2025) include. Includes only biomarker-confirmed ovulation cycles.

• scaled_cycleday_impute:
  Same as above, but also includes cycles with imputed ovulation (assigned to day –15 before next menses). Increases dataset coverage, but with lower precision. This maps -1 and +1 as adjacent days. This is what the supplementary analyses of Nagpal et al. (2025) include.

• scaled_cycleday_ov:
  A cycle time variable centered on ovulation (ovtoday == 1 → 0), ranging from –1 (start of follicular phase) to +1 (end of luteal phase). Biomarker-confirmed ovulation only. This maps -1 and +1 as adjacent days. This is what the supplementary analyses of Nagpal et al. (2025) include.

• scaled_cycleday_imp_ov:
  Same as above, but includes imputed ovulation (ovtoday_impute == 1) for cycles lacking confirmation. Centered on either confirmed or imputed ovulation. This maps -1 and +1 as adjacent days. This is what the supplementary analyses of Nagpal et al. (2025) include.

• cyclic_time:
  A continuous cycle time variable centered on menses onset (menses == 1 → 0), ranging from –1 (start of luteal phase) to +1 (ovulation). Includes only biomarker-confirmed ovulation cycles. This ensures that -1 and +1 are hormonally equivalent, rather than adjacent. This is what the main analyses of Nagpal et al. (2025) include. Other iterations of cycle time (scaled_cycleday) are being phased out.

• cyclic_time_impute:
  Same as above, but also includes cycles with imputed ovulation (assigned to day –15 before next menses). Increases dataset coverage, but with lower precision. This ensures that -1 and +1 are hormonally equivalent, rather than adjacent. This is what the main analyses of Nagpal et al. (2025) include.

• cyclic_time_ov:
  A cycle time variable centered on ovulation (ovtoday == 1 → 0), ranging from –1 (start of follicular phase) to +1 (end of luteal phase). Biomarker-confirmed ovulation only. This ensures that -1 and +1 are hormonally equivalent, rather than adjacent. This is what the main analyses of Nagpal et al. (2025) include.

• cyclic_time_imp_ov:
  Same as above, but includes imputed ovulation (ovtoday_impute == 1) for cycles lacking confirmation. Centered on either confirmed or imputed ovulation. This ensures that -1 and +1 are hormonally equivalent, rather than adjacent. This is what the main analyses of Nagpal et al. (2025) include.

• ovtoday_impute:
  A binary indicator identifying imputed ovulation days (value = 1), estimated as 15 days before the next menses onset.

• menses_impute: Present only when impute_next_menses = TRUE. A binary indicator identifying imputed menses onsets (value = 1) placed forward from a confirmed ovulation (at ovulation + 14 days by default), distinguishing them from observed onsets.

4.1.2 2. Data Checking Functions

These functions assess the completeness and coverage of your dataset:

• cycledata_check():
  Summarizes non-missing data by cycle phase (follicular, luteal) and overall. Also produces visualizations by participant and cycle phase.

• summary_ovulation():
  Reports ovulation detection status:

  • How many cycles had biomarker-confirmed ovulation
  • How many cycles required imputed ovulation using the day –15 rule (ovtoday_impute == 1), based on standard luteal phase length

4.1.3 3. Data Visualization Functions

Use these functions to plot outcomes across the menstrual cycle:

• cycle_plot():
  Creates a sample-level plot of symptom or outcome trajectories across the standardized cycle. Supports centering on either menses or ovulation, with flexible y-axis scaling.

• cycle_plot_individual():
  Provides participant-level plots across cycles for selected symptoms. Useful for individual inspection or diagnostic review.

Both functions: - Calculate per-cycle means and deviations - Allow centering on menses or ovulation - Are compatible with scaled time variables

4.1.4 4. Launching the Shiny App: launch_app()

This function launches an interactive Shiny application for uploading, checking, and analyzing menstrual cycle data using the package functions.

The app includes tools for: - Scaled cycle time computation - Data completeness checks - Visualization by person and symptom - PMDD, MRMD, and PME diagnosis using the Carolina Premenstrual Assessment Scoring System (C-PASS) - Please note that cycles utilized for diagnosis in the CPASS are defined by counting days forward (+10) and backward (-7) from menses onset.

C-PASS is implemented via the cpass R package. Please cite this package if you utilize the CPASS tool in the shiny app:
Symul L, Eisenlohr-Moul T (2025). cpass: PMDD and MRMD Diagnoses Following The Carolina Premenstrual Assessment Scoring System (C-PASS). R package version 0.1.0.

For more on C-PASS as a diagnostic tool, see:
Eisenlohr-Moul, T. A., Girdler, S. S., Schmalenberger, K. M., Dawson, D. N., Surana, P., Johnson, J. L., & Rubinow, D. R. (2017).
Toward the Reliable Diagnosis of DSM-5 Premenstrual Dysphoric Disorder: The Carolina Premenstrual Assessment Scoring System (C-PASS).
American Journal of Psychiatry, 174(1), 51–59. https://doi.org/10.1176/appi.ajp.2016.15121510 Please cite this paper if you utilize the CPASS tool in the shiny app.

The shiny app is also accessible online here:
👉 https://menstrualcycledata.shinyapps.io/shiny/

4.2.1 A second dataset, for awkward cases

cycledata is an ordinary diary: every cycle falls inside 21-35 days, every participant’s first recorded period comes before their first rating, and no cycle is left open at the end. That makes it a clean demonstration of PACTS, but it also means the package’s optional settings – imputing an ovulation before a first recorded period, or a next menses onset from a confirmed ovulation, or widening the phase-length bounds – change nothing at all when run on it.

A second demo dataset exists for those. It holds 15 participants, each one a different awkward case, named in words in a case column.

cases = unique(cycledata_special[, c("id", "case")])
writeLines(sprintf("%2d  %s", cases$id, cases$case))
#>  1  ordinary: three complete cycles, ovulation confirmed
#>  2  8 rated days before the first recorded menses onset
#>  3  30 rated days before the first onset: how many scale is a setting, not a cap
#>  4  confirmed ovulation among the leading days: scales with no setting
#>  5  confirmed ovulation at the end, no closing menses onset
#>  6  an 18-day cycle, ovulation not confirmed: under lower_cyclength_bound
#>  7  a 42-day cycle, ovulation not confirmed: over upper_cyclength_bound
#>  8  confirmed ovulation 22 days before the next onset: past luteal_phase_max_days
#>  9  confirmed ovulation 26 days after the onset: past follicular_phase_max_days
#> 10  confirmed ovulation 5 days before the next onset: under luteal_phase_min_days
#> 11  an 11-day hole in the diary: rows come back with no rating
#> 12  ovulation confirmed on no day: imputed columns only
#> 13  no menses onset and no confirmed ovulation: nothing can scale
#> 14  confirmed ovulation but no onset ever: one day scales, or a whole luteal phase
#> 15  both opt-in rules at once, at opposite ends of one diary

Use it to see what an optional setting actually does before applying it to your own data. ?cycledata_special documents each case, including which are recoverable and by which setting, and the “Preparing Your Data” vignette works through several of them. Like cycledata, it is simulated and represents no particular outcome or population.

4.3.1 Before and After: Filling Missing Dates Within a Cycle

Below is an example showing how pacts_scaling() fills in rows for missing dates in a real participant’s cycle.

Original data (cycle_df)

• This participant skipped reporting on 2024-01-23.
• The time series is incomplete for cycle scaling until the gap is filled.

After scaling (cycle_df_scaled)

• The function inserts a placeholder row for 2024-01-23 with NA in the symptom column.
• This ensures complete daily continuity from menses to the end of the cycle.

4.3.2 Summary of Row Counts

By ensuring complete daily coverage, pacts_scaling() allows downstream modeling and visualization tools to interpret menstrual cycle time as a continuous, standardized variable across individuals and cycles.

4.3.3 Cycle time variables are added in cycle_df_scaled

names(cycle_df_scaled)
#>  [1] "id"                                "date"                             
#>  [3] "menses"                            "ovtoday"                          
#>  [5] "symptom"                           "daterated"                        
#>  [7] "m2mcount"                          "mcyclength"                       
#>  [9] "cycle_incomplete"                  "cyclenum"                         
#> [11] "mcyclength_complete"               "ovtoday_impute"                   
#> [13] "scaled_cycleday"                   "scaled_cycleday_ov"               
#> [15] "scaled_cycleday_impute"            "scaled_cycleday_imp_ov"           
#> [17] "cyclic_time"                       "cyclic_time_impute"               
#> [19] "cyclic_time_impute_extended_phase" "cyclic_time_ov"                   
#> [21] "cyclic_time_imp_ov"                "cyclic_time_imp_ov_extended_phase"
#> [23] "luteal_length"

• scaled_cycleday:
  A continuous cycle time variable centered on menses onset (menses == 1 → 0), ranging from –1 (start of luteal phase) to +1 (ovulation). Includes only biomarker-confirmed ovulation cycles.

• scaled_cycleday_impute:
  Same as above, but also includes cycles with imputed ovulation (assigned to day –15 before next menses). Increases dataset coverage, but with lower precision.

• scaled_cycleday_ov:
  A cycle time variable centered on ovulation (ovtoday == 1 → 0), ranging from –1 (start of follicular phase) to +1 (end of luteal phase). Biomarker-confirmed ovulation only.

• scaled_cycleday_imp_ov:
  Same as above, but includes imputed ovulation (ovtoday_impute == 1) for cycles lacking confirmation. Centered on either confirmed or imputed ovulation.

• ovtoday_impute:
  A binary indicator identifying imputed ovulation days (value = 1), estimated as 15 days before the next menses onset.

• menses_impute: Present only when impute_next_menses = TRUE. A binary indicator identifying imputed menses onsets (value = 1) placed forward from a confirmed ovulation (ovulation + 14 days by default).

• m2mcount: A forward day count starting from each menses onset (day of onset = 0), reset at every new cycle.

• cyclenum: Which complete menses-to-menses cycle, in sequence, this row belongs to, per id. NA on cycle_incomplete == 1 rows (a still-open trailing cycle isn’t a “complete cycle” to number yet).

• mcyclength: Indicates the length of a menses-to-menses cycle an observation belongs to. On cycle_incomplete == 1 rows this is not a cycle length – it is the number of days observed so far in a still-open trailing cycle. Use mcyclength_complete instead if you want to threshold cycle length directly (see Section 3, “Cycle Length Inclusion Criteria”, above).

• mcyclength_complete: Same as mcyclength, but NA whenever cycle_incomplete == 1. Prefer this column for any length restriction – an incomplete cycle’s NA here can never satisfy a numeric comparison, so it drops out automatically without a separate cycle_incomplete == 0 filter.

• luteal_length: Indicates the length of the luteal phase an observation belongs to, when ovulation was not imputed. This is computed only for cycles that meet the criteria to be scaled.

• cycle_incomplete A binary indicator identifying if the observation is part of an incomplete menses-to-menses cycle. On these rows mcyclength is days-observed-so-far, not a cycle length – mcyclength_complete is NA instead.

• cyclic_time A continuous cycle time variable centered on menses onset (menses == 1 → 0), ranging from –1 to +1. Unlike scaled_cycleday, -1 and + 1 are parameterized to have the same hormonal meaning of ovulation. Includes only biomarker-confirmed ovulation cycles.

• cyclic_time_impute Same as above, but also includes cycles with imputed ovulation (assigned to day –15 before next menses). Increases dataset coverage, but with lower precision.

• cyclic_time_impute_extended_phase: Binary flag (0/1) marking rows where cyclic_time_impute was filled by the last-resort phase-cap fallback described in Section 3 above (value 1), rather than by a confirmed or imputed-ovulation anchor (0). Report this rate alongside the imputed-ovulation rate for full transparency about coverage.

• cyclic_time_ov A cycle time variable centered on ovulation (ovtoday == 1 → 0), ranging from –1 to +1. Unlike scaled_cycleday_ov, -1 and + 1 are parameterized to have the same hormonal meaning of menses_onset. Includes only biomarker-confirmed ovulation.

• cyclic_time_imp_ov Same as above, but includes imputed ovulation (ovtoday_impute == 1) for cycles lacking confirmation. Centered on either confirmed or imputed ovulation.

• cyclic_time_imp_ov_extended_phase: The same flag as cyclic_time_impute_extended_phase, for the ovulation-centered cyclic_time_imp_ov column.

4.6.1 cycle_plot()

The cycle_plot() function generates a group-level visualization of symptom or outcome data across the standardized menstrual cycle. This function supports flexible centering (on menses or ovulation), person-centered or raw scaling, and optional smoothing via rolling averages. Rolling averages are supported by the package zoo.

Let’s take a look at the parameters of cycle_plot:

cycle_plot(
  data,
  symptom,
  centering = "menses",
  include_impute = TRUE,
  y_scale = "person-centered_roll",
  rollingavg = 5,
  align_val = "center",
  se = FALSE
)

• data:
  A dataframe containing the cycle-aligned data. Must include the relevant scaled cycle time variables (scaled_cycleday, scaled_cycleday_ov, cyclic_time, etc.) and the symptom column. This function will only work if your dataset has been run through pacts_scaling() first.

• symptom:
  A string specifying the name of the symptom or outcome variable to visualize.

• centering:
  Whether to center the x-axis on "menses" (default) or "ovulation". This determines which scaled cycle day variable is used.

• include_impute:
  Logical. If TRUE, uses both confirmed and imputed ovulation cycles. Set to FALSE to restrict the analysis to biomarker-confirmed cycles only.

• y_scale:
  Specifies how to scale the y-axis:

  • "person-centered": Mean person-centered values
  • "person-centered_roll": Smoothed person-centered values (default) that applies a rolling average to data, using zoo::rollapply
  • "means": Raw mean values of participants

• rollingavg:
  Number of days used for rolling average smoothing (default = 5). Computed using zoo::rollapply.

• align_val:
  Controls alignment of the rolling average window ("center", "left", or "right"). Passed to zoo::rollapply.

• se:
  Logical. If TRUE, includes a standard error ribbon around the group mean.

4.6.2 Note on Rolling Averages

If y_scale = "person-centered_roll" is selected, rolling averages are calculated using the zoo::rollapply() function, which applies a moving average over the time series, with the settings as:

zoo::rollapply(
  variable,         # vector of values
  rollingavg,       # size of the moving window (default = 5 days)
  FUN = function(x) mean(x, na.rm = TRUE), # specifies a mean function to apply to each rolling window of values, ignoring missing values
  align = "center", # align the output to the center of the window
  fill = NA,        # fill edges with NA where full window is not available
  partial = TRUE    # allow smaller windows at the beginning and end
)

If you use the "person-centered_roll" option for the y_scale, please cite the zoo package:
Zeileis A, Grothendieck G (2005). “zoo: S3 Infrastructure for Regular and Irregular Time Series.”
Journal of Statistical Software, 14(6), 1–27.
https://doi.org/10.18637/jss.v014.i06

4.6.3 Symptom, Menses-centered

cycle_plot_df_menses <- cycle_plot(
  cycle_df_scaled,
  "symptom",
  centering = "menses",
  include_impute = TRUE,
  y_scale = "person-centered_roll", 
  se = TRUE
)

cycle_plot() returns a list with three components:

1. data: The original data augmented with person-mean (.m) and deviation (.d) values for the symptom.

A “person-mean” is the average of all of a single participant’s repeated scores or measurements over time. It represents a way to summarize what is “typical” for a given person, and is used for “between-person” or “trait-level” analyses.

• By extension, “deviations from person-mean” represent the differences between a person’s score on a specific day (or other repeated measure interval) and their usual average (their person-mean). It tells us how much higher or lower they are compared to their normal, and is used in “within-person” or “state-level” analyses.

2. summary: A summary of the outcome variable at each 5% increment of the scaled cycle timeline. It also provides standard errors if se = TRUE is included in the function parameters.

3. plot: A ggplot2 object visualizing the trajectory of the symptom variable across the cycle.

cycle_plot_df_menses$data
#> # A tibble: 744 × 28
#>       id date       menses ovtoday symptom daterated  m2mcount mcyclength
#>    <int> <date>      <dbl>   <dbl>   <dbl> <date>        <dbl>      <dbl>
#>  1     1 2024-01-20      1       0       5 2024-01-20        1         24
#>  2     1 2024-01-21      0       0       5 2024-01-21        2         24
#>  3     1 2024-01-22      0       0       3 2024-01-22        3         24
#>  4     1 2024-01-23      0       0      NA 2024-01-23        4         24
#>  5     1 2024-01-24      0       0       2 2024-01-24        5         24
#>  6     1 2024-01-25      0       0       1 2024-01-25        6         24
#>  7     1 2024-01-26      0       0       1 2024-01-26        7         24
#>  8     1 2024-01-27      0       0       1 2024-01-27        8         24
#>  9     1 2024-01-28      0       0       3 2024-01-28        9         24
#> 10     1 2024-01-29      0       0       4 2024-01-29       10         24
#> # ℹ 734 more rows
#> # ℹ 20 more variables: cycle_incomplete <dbl>, cyclenum <int>,
#> #   mcyclength_complete <dbl>, ovtoday_impute <int>, scaled_cycleday <dbl>,
#> #   scaled_cycleday_ov <dbl>, scaled_cycleday_impute <dbl>,
#> #   scaled_cycleday_imp_ov <dbl>, cyclic_time <dbl>, cyclic_time_impute <dbl>,
#> #   cyclic_time_impute_extended_phase <int>, cyclic_time_ov <dbl>,
#> #   cyclic_time_imp_ov <dbl>, cyclic_time_imp_ov_extended_phase <int>, …

cycle_plot_df_menses$summary
#> # A tibble: 22 × 3
#>    cycleday_5perc mean_dev_roll     se
#>             <dbl>         <dbl>  <dbl>
#>  1           0          1.03    0.178 
#>  2           0.05       0.677   0.125 
#>  3           0.1        0.0248  0.182 
#>  4           0.15      -0.398   0.198 
#>  5           0.2       -0.907   0.132 
#>  6           0.25      -1.13    0.109 
#>  7           0.3       -1.14    0.0702
#>  8           0.35      -0.663   0.107 
#>  9           0.4       -0.00551 0.146 
#> 10           0.45       0.613   0.122 
#> # ℹ 12 more rows

cycle_plot_df_menses$plot
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_line()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_ribbon()`).

This plot is a ggplot object and can be edited using the package ggplot2 and saved using the function ggsave().

4.6.4 Symptom, Ovulation-centered

Here is the same data, centered on ovulation.

cycle_plot_df_ov <- cycle_plot(
  cycle_df_scaled,
  "symptom",
  centering = "ovulation",
  include_impute = TRUE,
  y_scale = "person-centered_roll", 
  se = TRUE
)

cycle_plot_df_ov$plot
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_line()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_ribbon()`).

4.7.1 cycle_plot_individual

The cycle_plot_individual() function generates cycle-specific plots and summaries for a given participant ID. This function is especially useful for quality control, or individualized data visualization. It produces one plot per cycle, centered on either menses or ovulation, and includes summary data for each. Rolling averages are supported by the package zoo.

Let’s take a look at the parameters of cycle_plot_individual:

cycle_plot_individual(
  data,
  id,
  symptoms,
  centering = "menses",
  y_scale = "person-centered",
  include_impute = TRUE,
  rollingavg = 5
)

• data:
  A dataframe containing scaled cycle data with at least id, cyclenum, and the relevant symptom variables. The dataframe must already include the appropriate scaled cycle time columns.

• id:
  A numeric or character value specifying the participant ID to visualize.

• symptoms:
  A character vector of one or more symptom variables to plot.

• centering:
  Center the plot on "menses" (default) or "ovulation". Determines which scaled time variable is used.

• y_scale:
  Specifies the type of y-axis transformation:

  • "person-centered": Deviation from participant mean
  • "person-centered_roll": Smoothed person-centered values (default) that applies a rolling average to data, using zoo::rollapply
  • "raw": Uncentered raw values
  • "roll": Rolling average of raw values

• include_impute:
  Logical. If TRUE, includes cycles with imputed ovulation.

• rollingavg:
  Number of days to use in the rolling average window (default = 5). Applies only when using "roll" or "person-centered_roll" scales. Computed using zoo::rollapply.

  Note: If you use the "person-centered_roll" or "roll" option for the y_scale, please cite the zoo package:
  Zeileis A, Grothendieck G (2005). “zoo: S3 Infrastructure for Regular and Irregular Time Series.”
  Journal of Statistical Software, 14(6), 1–27.
  https://doi.org/10.18637/jss.v014.i06

4.7.2 Symptom, Menses-centered

This is what the raw data of symptom looks like across the first (and only cycle) for ID = 2:

cycle_plot_menses_id_2 <- cycle_plot_individual(
  cycle_df_scaled,
  id = 2, 
  "symptom",
  centering = "menses",
  y_scale = "raw",
  include_impute = TRUE
  
)

cycle_plot_menses_id_2$symptom$Cycle_1$plot

As you can see, it is difficult to examine any trends. Smoothing using y_scale = "roll" or y_scale = "person-centered_roll" will help us understand what the trajectory looks like for participant ID = 2.

cycle_plot_menses_id_2 <- cycle_plot_individual(
  cycle_df_scaled,
  id = 2, 
  "symptom",
  centering = "menses",
  y_scale = "roll",
  include_impute = TRUE,
  rollingavg = 3
)

cycle_plot_menses_id_2$symptom$Cycle_1$plot

Note on Rolling Averages

If y_scale = "person-centered_roll" or y_scale = "roll" is selected, rolling averages are calculated using the zoo::rollapply() function, which applies a moving average over the time series, with the settings as:

zoo::rollapply(
  variable,         # vector of values
  rollingavg,       # size of the moving window (default = 5 days)
  FUN = function(x) mean(x, na.rm = TRUE), # specifies a mean function to apply to each rolling window of values, ignoring missing values
  align = "center", # align the output to the center of the window
  fill = NA,        # fill edges with NA where full window is not available
  partial = TRUE    # allow smaller windows at the beginning and end
)

4.7.3 cycle_plot_individual Output

cycle_plot_individual() returns a list with two components accessible for every menses-to-menses cycle that exists for an individual:

1. summary: A summary of the outcome variable at each 5% increment of the scaled cycle timeline with computed deviations/ rolling deviations.

“Deviations from person-mean” represent the differences between a person’s score on a specific day (or other repeated measure interval) and their usual average (their person-mean). It tells us how much higher or lower they are compared to their normal, and is used in “within-person” or “state-level” analyses.

2. plot: A ggplot2 object visualizing the trajectory of the symptom variable across the cycle.

summary and plot are accessible for each menses-to-menses cycle a participant ID has in the dataset.

4.7.4 Cycle_1 summary for ID = 2

cycle_plot_menses_id_2$symptom$Cycle_1$summary
#> # A tibble: 21 × 7
#>    cycleday_5perc mean_dev mean_dev_roll raw_sx sx_roll cycleday mcyclength
#>             <dbl>    <dbl>         <dbl>  <dbl>   <dbl>    <dbl>      <dbl>
#>  1           0       2.29          1.62     5      4.33     15           26
#>  2           0.05    0.286         1.29     3      4        16           26
#>  3           0.1    -0.214        -0.381    2.5    2.33     17.5         26
#>  4           0.15   -1.71         -1.71     1      1        19           26
#>  5           0.2    -1.71         -1.71     1      1        20           26
#>  6           0.25   -1.71         -1.71     1      1        21           26
#>  7           0.3    -1.71         -1.38     1      1.33     22           26
#>  8           0.35   -0.714        -1.05     2      1.67     23           26
#>  9           0.4     0.786         0.786    3.5    3.5      24.5         26
#> 10           0.45    2.29          2.29     5      5        26           26
#> # ℹ 11 more rows

4.7.5 Cycle_1 plot for ID = 2

cycle_plot_menses_id_2$symptom$Cycle_1$plot

Cycles can also be examined ovulation-centered.

cycle_plot_ov_id_2 <- cycle_plot_individual(
  cycle_df_scaled,
  id = 2, 
  "symptom",
  centering = "ovulation",
  y_scale = "roll",
  include_impute = TRUE,
  rollingavg = 3
)

cycle_plot_ov_id_2$symptom$Cycle_1$plot

All of these plots are ggplot objects and can be edited using the package ggplot2 and saved using the function ggsave().

5.2.1 Data Preparation

Here, we log-transform our outcome variable symptom, because oftentimes many biological and behavioral outcomes (e.g. hormone levels, reaction times, symptom scores) are right-skewed. Log-transforming these outcomes often results in a distribution that is more symmetric and approximately normal. However, you should always assess your model for heteroscedasticity and examine the distribution of residuals.

cycle_df_scaled$symptom_log = log(cycle_df_scaled$symptom + 1) #log-transforming our outcome variable symptom

We now create a dataframe called datSX that includes only complete cases for the variables cyclic_time_impute and symptom_log.

GAMMs in the mgcv package use listwise deletion, meaning the model will only include rows with no missing values across any variables in the formula. By pre-filtering the dataset using complete.cases(), we can determine exactly how many observations will be included in the model. This ensures transparency and avoids confusion when interpreting model results.

selected_vars <- c("cyclic_time_impute", "symptom_log" )
datSX <- cycle_df_scaled[complete.cases(cycle_df_scaled[selected_vars]), ]

Let’s run the GAMM!

datSX$id = as.factor(datSX$id) # ALWAYS factor id before putting it in a gam formula

gamm1 <- mgcv::gam(
  symptom_log ~ # outcome: log-transformed symptom score
    s(cyclic_time_impute, bs = "cc") + # fixed effect: population-average smooth of cycle time ("cc" = cyclic cubic spline)
    s(id, bs = 're') + # random intercept: per-person differences in average symptom level
    s(cyclic_time_impute, id, bs = "re"), # random slope: per-person linear tilt on the cycle trajectory
  knots= list(cyclic_time_impute = c(-1,1)), # tie -1 and +1 to the same hormonal timepoint (cyclic boundary)
  data = datSX, # data: complete cases only (listwise deletion)
  method = 'REML' # estimate smoothing parameters via REML
)

This model includes:

• symptom_log
  The outcome variable (log-transformed symptom).

• s(cyclic_time_impute)
  A smooth term modeling the fixed effect of cycle time (scaled, centered on menses) across the sample, where -1 and +1 both represent the same hormonal timepoint across the cycle (ovulation). This captures nonlinear, average trends in the symptom over the cycle.

• s(id, bs = 're')
  A random intercept for each participant. This accounts for individual differences in average symptom level.

• s(cyclic_time_impute, id, bs = 're') A random slope: allows each participant’s trajectory to be tilted (steeper, flatter, or reversed) linearly relative to the population-level smooth. Because the basis is "re" on both margins, this term is linear in cyclic_time_impute per person — it shifts the slope of the line, not the shape of the curve. It does not give each participant their own nonlinear cyclic trajectory; every participant shares the same population-level curve shape (s(cyclic_time_impute, bs = "cc")) and only its intercept and linear slope vary by person. (A true per-person nonlinear cyclic trajectory would require a factor-smooth interaction, e.g. s(cyclic_time_impute, id, bs = "fs", xt = list(bs = "cc")), which is a substantially larger model not used here.)

• data = datSX
  The dataset used for model fitting. It should contain only complete cases for all model variables.

• method = 'REML'
  Fits the model using Restricted Maximum Likelihood, which improves estimation of smoothness and variance components, especially in the presence of random effects.

5.3.1 Model Family & Link

Family: gaussian 
Link function: identity 

• The outcome is modeled using a Gaussian (normal) distribution with an identity link (i.e., predicted values are on the same scale as the outcome). This is appropriate for continuous, unbounded outcome variables.

5.3.2 Model Formula

symptom_log ~ s(cyclic_time_impute) + s(id, bs = "re") +
              s(cyclic_time_impute, id, bs = "re")

• This model includes:
  • A smooth term for population-level cyclic trends.
  • A random intercept smooth function for each participant.
  • A random slope giving each participant a linear tilt relative to the population-level cyclic curve (not a separate nonlinear shape per person).

5.3.3 Parametric Coefficients

(Intercept)  1.30963    0.01801   72.72   <2e-16 ***

• The intercept represents the average value of the log-transformed symptom when all smooth terms are at their reference level.

5.3.4 Smooth Terms

s(cyclic_time_impute)     edf = 7.47   F = 81.60   p < 2e-16 ***
s(id)                     edf = 13.71  F = 1.52    p = 0.00022 ***
s(cyclic_time_impute,id)  edf = 14.75  F = 1.51    p = 0.000375 ***

Each smooth term contributes to explaining nonlinear or hierarchical effects:

• s(cyclic_time_impute): Captures the overall nonlinear trend in symptom expression across the standardized cycle.
  • The high edf (≈7.5) reflects flexibility in the shape of the curve. Higher edf indicates more “wiggles” and complexity in the nonlinear pattern.
  • The very small p-value suggests this population-level pattern is strongly significant for nonlinearity
• s(id, bs = "re"): Models random intercept to account for baseline differences across individuals.
  • p = 0.0002 indicates that between-person variation in baseline symptom level is detectable at this sample size. Whether it is large enough to matter is a separate question, answered by the size of the variance component rather than by the p-value.
• s(cyclic_time_impute, id, bs = "re"): Captures each participant’s linear tilt (steeper, flatter, or reversed) relative to the population-level curve — not a separate nonlinear shape per person.
  • p = 0.000375 indicates detectable heterogeneity in the linear slope of symptom change across individuals, on top of the shared population-level nonlinear curve. As above, detectable is not the same as large.

5.3.5 Model Fit

R-sq.(adj) =  0.528
Deviance explained = 55.5%
-REML = 125.69
Scale est. = 0.077082
n = 611

• Adjusted R² of 0.528: The model explains over 52% of the variance in the log-transformed symptom outcome. It is roughly the squared correlation between observed and fitted values, adjusted for model complexity. It is used to assess overall model fit when working with Gaussian models.

• Deviance explained: Another metric showing over 55% of outcome variability is accounted for. Deviance is a likelihood-based measure, more general than R², but comes from comparing the model’s deviance (a measure of error) to that of a null model. It should be used when comparing models or working with non-normal outcomes.

• REML: Indicates the restricted maximum likelihood used in model fitting—important for selecting smoothing parameters in GAMs.

• n = 611: The number of complete observations included in the model after listwise deletion.

These results suggest that on average, symptom varies significantly across the cycle in a nonlinear way (shared population-level curve), and that individuals also differ significantly in their linear slope relative to that curve – i.e. some participants’ symptom trajectories tilt more steeply up or down across the cycle than others, on top of the same underlying nonlinear shape.

5.4.1 First, we need to compute predicted or model-implied values

plotdat <- expand.grid(cyclic_time_impute = seq(-1, 1, by = 0.05),
                      id = 0) # setting id = 0 suppresses random effects, to model just the fixed effect (sample-wide) of your outcome across the cycle

# Predict using the model for each dataset and add predictions
pred <- marginaleffects::predictions(gamm1, newdata = plotdat, type = "response", transform = function(x) exp(x) - 1) # applying a transform function, to undo the log transformation on symptom. The transform can be removed if your outcome was not log-transformed 
plotdat$estimate = pred$estimate
plotdat$conf.low = pred$conf.low
plotdat$conf.high = pred$conf.high

5.4.2 Plotting model-implied values

# Plotting
gamplot <- ggplot(plotdat, aes(x = cyclic_time_impute, y = estimate)) +
  scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, by = 0.50), 
                     labels = c("Ovulation", "50%L", "Menses Onset", "50%F", "Ovulation")) +
  labs(x = "", y = "Symptom") + # You can change the y-axis label to reflect your outcome
  
  geom_rect(xmin =-0, xmax = 0.08, ymin = -Inf, ymax = Inf,
            fill = "grey70", alpha = 0.2, color = "white") +
  geom_rect(xmin = 0.92, xmax = 1, ymin = -Inf, ymax = Inf,
            fill = "grey87", alpha = 0.2, color = "white") +
  geom_rect(xmin = -1, xmax = -0.92, ymin = -Inf, ymax = Inf,
            fill = "grey87", alpha = 0.2, color = "white") +
  geom_line(size = 1, show.legend = TRUE) +
  # Adding CI ribbon with translucent light grey color
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "lightgrey", alpha = 0.3) +
  theme_minimal()
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

# Print the plot
print(gamplot)