## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = FALSE
)

## ----libs, message = FALSE, warning = FALSE-----------------------------------
library(menstrualcycleR)
library(dplyr)

## ----shape--------------------------------------------------------------------
head(cycledata, 4)

## ----dupes--------------------------------------------------------------------
dups <- cycledata %>%
  count(id, daterated) %>%
  filter(n > 1)
nrow(dups)

## ----coverage-----------------------------------------------------------------
scaled <- pacts_scaling(cycledata, id = id, date = daterated,
                        menses = menses, ovtoday = ovtoday)

c(all_rows             = nrow(scaled),
  rows_with_an_outcome = sum(!is.na(scaled$symptom)))

## ----sources------------------------------------------------------------------
period_log <- data.frame(
  id           = c(1, 1, 1, 2, 2, 2),
  period_start = as.Date(c("2024-01-03", "2024-01-31", "2024-02-27",
                           "2024-01-08", "2024-02-10", "2024-03-09"))
)

set.seed(1)
wearable <- data.frame(
  id   = rep(c(1, 2), each = 75),
  date = c(seq(as.Date("2024-01-01"), by = "day", length.out = 75),
           seq(as.Date("2024-01-05"), by = "day", length.out = 75)),
  resting_hr = round(c(rnorm(75, 58, 2), rnorm(75, 62, 2)), 1)
)

head(period_log, 3)
head(wearable, 3)

## ----join---------------------------------------------------------------------
daily <- wearable %>%
  left_join(mutate(period_log, menses = 1L),
            by = c("id", "date" = "period_start")) %>%
  mutate(
    menses  = coalesce(menses, 0L),
    ovtoday = 0L          # no ovulation biomarker in this example
  )

daily %>% filter(menses == 1) %>% head(3)

## ----scale-example------------------------------------------------------------
daily_scaled <- pacts_scaling(daily, id = id, date = date,
                              menses = menses, ovtoday = ovtoday)

daily_scaled %>%
  select(id, date, resting_hr, cyclic_time_impute, ovtoday_impute) %>%
  filter(!is.na(cyclic_time_impute)) %>%
  head(4)

## ----checks-------------------------------------------------------------------
checks <- c(
  date_is_Date        = inherits(daily$date, "Date"),
  no_duplicate_rows   = !any(duplicated(daily[c("id", "date")])),
  menses_is_binary    = all(daily$menses %in% c(0, 1)),
  has_ovtoday_column  = "ovtoday" %in% names(daily),
  every_id_has_onset  = all(tapply(daily$menses, daily$id, sum) > 0)
)
checks

## ----coverage-by-variable-----------------------------------------------------
vars <- c("cyclic_time", "cyclic_time_impute", "cyclic_time_ov", "cyclic_time_imp_ov")
sapply(vars, function(v) sum(!is.na(scaled[[v]])))

## ----sensitivity--------------------------------------------------------------
# The confirmed-only subset is the rows cyclic_time covers.
c(rows_confirmed_only = sum(!is.na(scaled$cyclic_time)),
  rows_incl_imputed   = sum(!is.na(scaled$cyclic_time_impute)))

## ----extended-filter, eval = FALSE--------------------------------------------
# strict <- subset(scaled, cyclic_time_impute_extended_phase == 0)

## ----special-cases------------------------------------------------------------
cases <- unique(cycledata_special[, c("id", "case")])
writeLines(sprintf("%2d  %s", cases$id, cases$case))

## ----special-leading----------------------------------------------------------
off <- pacts_scaling(cycledata_special, id = id, date = daterated,
                     menses = menses, ovtoday = ovtoday)

on <- pacts_scaling(cycledata_special, id = id, date = daterated,
                    menses = menses, ovtoday = ovtoday,
                    impute_leading_ovulation = TRUE)

covered <- function(x) !is.na(x$cyclic_time_impute)

c(days_gained     = sum(covered(on)) - sum(covered(off)),
  with_a_rating   = sum(covered(on) & !is.na(on$symptom)) -
                    sum(covered(off) & !is.na(off$symptom)),
  rows_fabricated = nrow(on) - nrow(off))

## ----special-placement--------------------------------------------------------
sapply(c(5, 15, 25, 40), function(k) {
  x <- pacts_scaling(cycledata_special, id = id, date = daterated,
                     menses = menses, ovtoday = ovtoday,
                     impute_leading_ovulation = TRUE,
                     leading_ovulation_luteal_days = k)
  sum(covered(x)) - sum(covered(off))
})

## ----special-floor------------------------------------------------------------
most <- pacts_scaling(cycledata_special, id = id, date = daterated,
                      menses = menses, ovtoday = ovtoday,
                      impute_leading_ovulation = TRUE, impute_next_menses = TRUE,
                      lower_cyclength_bound = 20, upper_cyclength_bound = 43,
                      luteal_phase_min_days = 4)

unique(most$case[!is.na(most$id) & is.na(most$cyclic_time_impute) &
                   !most$id %in% most$id[covered(most)]])

## ----special-fill, eval = FALSE-----------------------------------------------
# library(dplyr)
# library(tidyr)
# most |> group_by(id) |> fill(case, .direction = "downup") |> ungroup()

## ----ovsummary----------------------------------------------------------------
summary_ovulation(scaled)$ovstatus_total

## ----extended-----------------------------------------------------------------
mean(scaled$cyclic_time_impute_extended_phase == 1, na.rm = TRUE)