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

## ----remotes, eval = FALSE----------------------------------------------------
# install.packages("remotes")

## ----loading-remotes, warning = FALSE, eval = FALSE---------------------------
# library(remotes)

## ----menstrualcycleR-package, eval = FALSE------------------------------------
# remotes::install_github("eisenlohrmoullab/menstrualcycleR")

## ----loading menstrualcycleR-package------------------------------------------
library(menstrualcycleR)

## ----tidyverse, eval=FALSE----------------------------------------------------
# install.packages(c("dplyr", "ggplot2"))

## ----loading-tidyverse, warning=FALSE-----------------------------------------
library(dplyr)
library(ggplot2)

## ----loading-pcks, warning=FALSE, eval = FALSE--------------------------------
# library(tidyverse)

## ----gamm-packages-install, eval = FALSE--------------------------------------
# install.packages("mgcv")
# install.packages("marginaleffects")

## ----gamm-packages-libraries--------------------------------------------------
library(mgcv)
library(marginaleffects)

## ----demo dataset-------------------------------------------------------------
cycle_df = cycledata
dim(cycle_df)
head(cycle_df)

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

## ----pacts-scaling------------------------------------------------------------
cycle_df_scaled = pacts_scaling(data = cycle_df, id = id, date = daterated, menses = menses, ovtoday = ovtoday, lower_cyclength_bound = 21, upper_cyclength_bound = 35)

## ----view-data----------------------------------------------------------------
names(cycle_df_scaled)

## ----cycledata_check----------------------------------------------------------
checkdata = cycledata_check(cycle_df_scaled, symptom_columns = c("symptom"))

## ----cycledata_check_by_id----------------------------------------------------
checkdata$by_id

## ----cycledata_check_overall--------------------------------------------------
checkdata$overall

## ----cycledata_check_plot,fig.width=7, fig.height=6---------------------------
checkdata$data_symptom_plots

## ----ovulation_summary--------------------------------------------------------
ov_summary = summary_ovulation(cycle_df_scaled)

## ----ovulation_summary_total--------------------------------------------------
ov_summary$ovstatus_total

## ----ovulation_summary_id-----------------------------------------------------
ov_summary$ovstatus_id

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

## ----eval = FALSE-------------------------------------------------------------
# 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
# )

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


## ----cycle_plot_menses_raw_data-----------------------------------------------
cycle_plot_df_menses$data

## ----cycle_plot_menses_raw_summary--------------------------------------------
cycle_plot_df_menses$summary

## ----cycle_plot_menses_raw_plot, fig.width=7, fig.height=6--------------------
cycle_plot_df_menses$plot

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


## ----cycle_plot_ov_roll_plot_access, fig.width=7, fig.height=6----------------
cycle_plot_df_ov$plot

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

## ----cycle_plot_menses_id_2_raw_plot, fig.width=7, fig.height=5---------------
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

## ----cycle_plot_menses_id_2_roll_plot, fig.width=7, fig.height=5--------------
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

## ----eval = FALSE-------------------------------------------------------------
# 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
# )

## ----cycle_plot_menses_id_2_summary-------------------------------------------
cycle_plot_menses_id_2$symptom$Cycle_1$summary

## ----cycle_plot_menses_id_2_plot_access, fig.width=7, fig.height=5------------
cycle_plot_menses_id_2$symptom$Cycle_1$plot

## ----cycle_plot_ov_id_2_summary, fig.width=7, fig.height=5--------------------
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

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

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

## -----------------------------------------------------------------------------
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
)

## -----------------------------------------------------------------------------
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

## ----gam_plot, fig.width=7, fig.height=5--------------------------------------
# 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()

# Print the plot
print(gamplot)