2  SRI Analysis

2.1 Overview

2.2 Set the Environment

2.2.1 Load Packages

Code
library(actverse) # github.com/danielvartan/actverse
library(askpass)
library(broom)
library(checkmate)
library(cli)
library(colorspace)
library(effectsize)
library(fs)
library(geepack)
library(ggplot2)
library(glmtoolbox)
library(googlesheets4)
library(here)
library(hms)
library(infer)
library(janitor)
library(lockr) # github.com/danielvartan/lockr
library(lubridate)
library(lubritime) # github.com/danielvartan/lubritime
library(magrittr)
library(orbis) # github.com/danielvartan/orbis
library(patchwork)
library(parsnip)
library(performance)
library(plotr) # github.com/danielvartan/plotr
library(prettycheck) # github.com/danielvartan/prettycheck
library(psychometric)
library(pwr)
library(pwrss)
library(readr)
library(rutils) # github.com/danielvartan/rutils
library(scaler) # github.com/danielvartan/scaler
library(scales)
library(sodium)
library(stats)
library(stringr)
library(summarytools)
library(tidyr)
library(tsibble)

2.2.2 Load Custom Functions

Code
source(here("R", "anonymize_id.R"))
source(here("R", "ga.R"))

2.2.3 Set Data Directories

Code
data_dir <- here("data")
processed_data_dir <- here(data_dir, "processed")
actigraphy_processed_data_dir <- here(data_dir, "processed", "actigraphy")
sri_data_dir <- here(data_dir, "sri")
sleep_prop_dir <- here(data_dir, "sleep-prop")
lux_data_dir <- here(data_dir, "lux")
Code
for (i in ls(pattern = "_dir$")) {
  if (!dir_exists(get(i))) dir_create(get(i), recurse = TRUE)
}

2.2.4 Set Keys

Code
gs4_auth(
  cache = here(".secrets"),
  email = TRUE,
  use_oob = FALSE
)
Code
salt <- Sys.getenv("PREGNANCY_SALT") # askpass()

2.3 Set Exclusion Criteria

For SRI Data:

Last updated: 2025-09-19

  • n = 44

  • 5 days (drops 6) + 75% (drops 11): 27

  • 5 days (drops 6) + 62.5% (drops 2): 36

  • 4 days (drops 3) + 62.5% (drops 2): 39 # Ok

  • 4 days (drops 3) + 50% (drops 1): 40 # Ok

1 of the 40 was dropped when filtering the gestacional age range.

Code
min_days <- 4
min_valid_data <- 0.5 # By `sri` and gestational age (they are related)

2.4 Download Actigraphy Processed Data from OSF

See the osf-routines.qmd Quarto notebook to download and decrypt the actigraphy processed data from OSF.

2.5 Download and Import the Control Form Data

Code
control_form_data <-
  "13wtDr4fRD1wJSM-qdLOdGSchF8oXU1bG0Jtccu24GhY" |>
  read_sheet(
    sheet = "Dataset",
    col_types = "c",
    na = c("", "NA")
  )

2.6 Download and Import the Field Form Data

Code
field_form_data <-
  "1tY_TT0nPXFsErYQS2VED0-hqTzgHudA9BfhRhYG19fw" |>
  read_sheet(
    sheet = "Dataset",
    col_types = "c",
    na = c("", "NA")
  )

2.7 Filter and Tidy the Control Form Data

Code
control_form_data <-
  control_form_data |>
  `names<-`(paste0("X", seq_along(control_form_data))) |>
  dplyr::select(X5, X6, X12, X39, X40) |>
  rename(
    cpf = X5,
    id = X6,
    birth_center_delivery = X12,
    ultrasound_ga = X39,
    ultrasound_date = X40
  ) |>
  mutate(
    birth_center_delivery = case_when(
      tolower(trimws(birth_center_delivery)) == "sim" ~ TRUE,
      tolower(trimws(birth_center_delivery)) == "não" ~ FALSE,
      TRUE ~ NA
    ),
    ultrasound_date = ymd(ultrasound_date),
    ga_start = if_else(
      is.na(ultrasound_ga),
      as.Date(NA),
      ga_start(
        ultrasound = ultrasound_date,
        ga =
          str_extract(ultrasound_ga, "^\\d{1,2}") |>
          as.integer() |>
          weeks() |>
          add(
            str_extract(ultrasound_ga, "\\d(?= di)") |>
            as.integer() |>
            days()
          )
      )
    )
  ) |>
  dplyr::select(id, birth_center_delivery, ga_start)
control_form_data |> glimpse()
#> Rows: 54
#> Columns: 3
#> $ id                    <chr> "9e858cd8", "4b44f0ff", "a053ce64", "02250df5…
#> $ birth_center_delivery <lgl> FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, NA, T…
#> $ ga_start              <date> 2022-09-03, 2022-08-16, 2022-09-03, 2022-10-…

2.8 Filter and Tidy the Field Form Data

Code
field_form_data <-
  field_form_data |>
  `names<-`(paste0("X", seq_along(field_form_data))) |>
  dplyr::select(
    X1, X4, X5, X11, X91, X92, X95, X96, X97, X99, X101, X102, X106, X119,
    X118, X125, X126
  ) |>
  rename(
    timestamp = X1,
    birth_date = X4,
    cpf = X5,
    weight_before = X96,
    weight_current = X95,
    height = X97,
    ethnicity = X106,
    education = X118,
    family_income = X99,
    health_plan = X101,
    solo = X102,
    exercise = X119,
    work = X92,
    study = X91,
    gestations = X11,
    deliveries = X125,
    abortions = X126
  ) |>
  mutate(
    id = anonymize_id(cpf, salt),
    timestamp = mdy_hms(timestamp, tz = "America/Sao_Paulo"),
    birth_date = dmy(birth_date),
    age = scaler:::age(birth_date, timestamp),
    ethnicity = factor(
      ethnicity,
      levels = c(
        "Indígena",
        "Preta",
        "Parda",
        "Amarela",
        "Branca"
      ),
      ordered = TRUE
    ),
    education = factor(
      education,
      levels = c(
        "Não frequentou a escola",
        "Fundamental incompleto",
        "Fundamental completo",
        "Ensino médio incompleto",
        "Ensino médio completo",
        "Ensino superior incompleto",
        "Ensino superior completo",
        "Mestrado incompleto",
        "Mestrado completo",
        "Doutorado incompleto",
        "Doutorado completo",
        "Pós-doutorado incompleto",
        "Pós-doutorado completo"
      ),
      ordered = TRUE
    ),
    gestations = case_match(
      gestations,
      c("Sim", "0") ~ "1",
      "Não" ~ "0",
      .default = gestations
    ),
    deliveries = case_when(
      is.na(deliveries) ~ "0",
      TRUE ~ deliveries
    ),
    abortions = case_when(
      is.na(abortions) ~ "0",
      TRUE ~ abortions
    )
  ) |>
  mutate(
    across(
      .cols = all_of(
        c(
          "weight_before", "weight_current", "height",
          "family_income", "gestations", "deliveries",
          "abortions"
        )
      ),
      .fns = as.numeric
    )
  ) |>
  mutate(
    across(
      .cols = all_of(
        c("health_plan", "solo", "exercise", "work", "study")
      ),
      .fns = function(x) {
        case_match(
          x,
          "Sim" ~ TRUE,
          "Não" ~ FALSE
        )
      }
    )
  ) |>
  mutate(
    gestations = case_when(
      id == "ae5fc5cd" ~ 3,
      TRUE ~ gestations
    ),
  ) |>
  mutate(
    bmi_before = weight_before / ((height / 100)^2),
    bmi_current = weight_current / ((height / 100)^2)
  ) |>
  dplyr::select(-cpf) |>
  relocate(
    id, timestamp, birth_date, age, weight_before, weight_current, height, bmi_before, bmi_current, family_income
  )
field_form_data |> glimpse()
#> Rows: 74
#> Columns: 20
#> $ id             <chr> "9e858cd8", "a053ce64", "b37bd425", "db850249", "022…
#> $ timestamp      <dttm> 2023-05-21 11:18:49, 2023-05-26 17:55:52, 2023-06-3…
#> $ birth_date     <date> 1991-03-22, 2001-07-25, 1989-02-23, 1998-01-12, 199…
#> $ age            <dbl> 32.16388889, 21.83611111, 34.35277778, 25.47222222, …
#> $ weight_before  <dbl> 70, 50, 87, 70, 64, 84, 49, 60, 75, 66, 64, 72, 92, …
#> $ weight_current <dbl> 86, 75, 94, 76, 78, 87, 58, 91, 89, 79, 76, 89, 110,…
#> $ height         <dbl> 160, 156, 174, 169, 153, 168, 157, 164, 164, 159, 16…
#> $ bmi_before     <dbl> 27.34375000, 20.54569362, 28.73563218, 24.50894577, …
#> $ bmi_current    <dbl> 33.59375000, 30.81854043, 31.04769454, 26.60971255, …
#> $ family_income  <dbl> 4000, 2800, 6000, 5000, 3000, 4000, 3000, 9000, 6000…
#> $ gestations     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1…
#> $ study          <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALS…
#> $ work           <lgl> TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TR…
#> $ health_plan    <lgl> FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE…
#> $ solo           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
#> $ ethnicity      <ord> Branca, Parda, Branca, Branca, Parda, Parda, Branca,…
#> $ exercise       <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ education      <ord> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ deliveries     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0…
#> $ abortions      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0…

2.9 Import, Tidy and Transform the Processed Actigraphy Data

2.9.1 Get Paths to Processed Actigraphy Data Files

Code
files <-
  dir_ls(
    path = actigraphy_processed_data_dir,
    type = "file",
    regexp = ".txt$"
  )

2.9.2 Import, Tidy, Transform, and Write Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_actigraphy-processed-data.txt")

  i_ga_start <-
    control_form_data |>
    filter(id == i_id) |>
    pull(ga_start)

  i_data <-
    i |>
    read_acttrust(tz = "America/Sao_Paulo") |>
    mutate(
      ga =
        i_ga_start |>
        ga_point(timestamp) |>
        as.numeric() |>
        as.duration() |>
        add(
          as_hms(timestamp) |>
            as.numeric() |>
            as.duration()
        ),
      ga_week =
        ga |>
        as.numeric() |>
        divide_by(
          ddays(7) |>
            as.numeric()
        ) |>
        as.integer()
    ) |>
    suppressMessages() |>
    suppressWarnings()

  i_data |>
    write_rds(
      path(
        actigraphy_processed_data_dir,
        paste0(i_id, "_actigraphy-processed-data.rds")
      )
    )

  cli_progress_update()
}

cli_progress_done()

2.10 Compute SRI

2.10.1 Delete Old SRI Data Files

Code
old_files <-
  dir_ls(
    path = sri_data_dir,
    type = "file",
    regexp = "sri"
  )

if (length(old_files) > 0) file_delete(old_files)

2.10.2 Get Paths to Tidy and Transformed Actigraphy Data Files

Code
files <-
  dir_ls(
    path = actigraphy_processed_data_dir,
    type = "file",
    regexp = "_actigraphy-processed-data.rds$"
  )

2.10.3 Import Tidy and Transformed Actigraphy Data, Compute SRI and Write SRI Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_actigraphy-processed-data.rds")

  i_data <-
    i |>
    read_rds() |>
    as_tibble() |>
    group_by(ga_week) |>
    mutate(
      ga_week_available_days = (n() * 60) / as.numeric(ddays(1))
    ) |>
    ungroup(ga_week) |>
    filter(ga_week_available_days >= min_days) |>
    as_tsibble(index = timestamp)

  if (nrow(i_data) == 0) {
    cli_alert_info(
      paste0(
        "{.strong {col_red(i_id)}} was dropped because it has no gestational ",
        "week with at least {min_days} days of data."
      )
    )

    cli_progress_update()

    next
  }

  i_sri_data <- tibble()

  for (j in unique(i_data$ga_week)) {
    j_data <- i_data |> filter(ga_week == j)
    j_data_first_index <- which(i_data$ga_week == j)[1]

    if (j_data_first_index > 1) {
      j_data <-
        i_data |>
        slice(seq(j_data_first_index - 1440, j_data_first_index - 1)) |>
        as_tibble() |>
        bind_rows(j_data |> as_tibble()) |>
        as_tsibble(index = timestamp)
    }

    i_sri_data <-
      i_sri_data |>
      bind_rows(
        j_data |>
          sri(
            sleeping_states = c(1, 2),
            awake_states = 0,
            min_data = min_valid_data
          ) |>
          mutate(
            id = i_id,
            ga_week = j
          ) |>
          as_tibble()
      )
  }

  i_sri_data <-
    i_sri_data |>
    drop_na(sri) |>
    group_by(ga_week) |>
    mutate(dummy = n() >= (1440 * min_valid_data)) |>
    ungroup() |>
    filter(dummy) |>
    dplyr::select(-dummy)

  if (nrow(i_sri_data) == 0) {
    cli_alert_info(
      paste0(
        "{.strong {col_red(i_id)}} was dropped because it has no gestational ",
        "week with at least {round(min_valid_data * 100, 2)}% of SRI data."
      )
    )

    cli_progress_update()

    next
  }

  i_sri_data |>
    write_rds(
      path(
        sri_data_dir,
        paste0(i_id, "_sri-data.rds"))
      )

  cli_progress_update()
}

cli_progress_done()

2.10.4 Group and Transform SRI Data

2.10.4.1 Get Paths to SRI Data Files

Code
files <-
  dir_ls(
    path = sri_data_dir,
    type = "file",
    regexp = "_sri-data.rds$"
  )

2.10.4.2 Group and Transform SRI Data

Code
sri_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  sri_data <-
    i |>
    read_rds() |>
    mutate(
      time =
        time |>
        as.numeric() |>
        divide_by(60) |>
        as.integer() |>
        dminutes() |>
        as.numeric() |>
        as_hms()
    ) |>
    bind_rows(sri_data)

  cli_progress_update()
}

cli_progress_done()
Code
sri_data <-
  sri_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
sri_data |> glimpse()
#> Rows: 161,181
#> Columns: 8
#> $ time           <time> 00:00:00, 00:01:00, 00:02:00, 00:03:00, 00:04:00, 0…
#> $ state          <list> <Sleeping, Awake, Sleeping, Awake, Sleeping, Sleepi…
#> $ previous_state <list> <NA, Sleeping, Awake, Sleeping, Awake, Sleeping, Sl…
#> $ agreement      <list> <NA, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE>, <NA,…
#> $ sri            <dbl> -66.66666667, -66.66666667, -66.66666667, -66.666666…
#> $ valid_data     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ id             <chr> "e06521cc", "e06521cc", "e06521cc", "e06521cc", "e06…
#> $ ga_week        <fct> 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, …

2.10.5 Compute SRI Mean Data by Gestational Week

2.10.5.1 Get Paths to SRI Data Files

Code
files <-
  dir_ls(
    path = sri_data_dir,
    type = "file",
    regexp = "_sri-data.rds$"
  )

2.10.5.2 Summarize SRI Means by Gestational Week and Write SRI Mean Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_sri-data.rds")

  i_data <-
    i |>
    read_rds() |>
    as_tibble() |>
    group_by(ga_week) |>
    summarize(sri = mean(sri, na.rm = TRUE)) |>
    mutate(id = i_id) |>
    dplyr::select(id, ga_week, sri)

  i_data |>
    write_rds(
      path(
        sri_data_dir,
        paste0(i_id, "_sri-mean-data.rds")
      )
    )

  cli_progress_update()
}

cli_progress_done()

2.10.6 Group and Transform SRI Mean Data

2.10.6.1 Get Paths to SRI Mean Data Files

Code
files <-
  dir_ls(
    path = sri_data_dir,
    type = "file",
    regexp = "_sri-mean-data.rds$"
  )

2.10.6.2 Import and Group SRI Mean Data

Code
sri_mean_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  sri_mean_data <-
    i |>
    read_rds(i) |>
    bind_rows(sri_mean_data)

  cli_progress_update()
}

cli_progress_done()

2.10.6.3 Transform SRI Mean Data

Code
sri_mean_data <-
  sri_mean_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
sri_mean_data |> glimpse()
#> Rows: 118
#> Columns: 3
#> $ id      <chr> "e06521cc", "4cda994c", "e06521cc", "a8194247", "9e43829a",…
#> $ ga_week <fct> 33, 33, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35,…
#> $ sri     <dbl> 73.75293015, 50.77314815, 68.65350746, 75.85648148, 53.3755…

2.11 Get SRI Data Unique IDs

Code
sri_data_ids <- sri_data |> pull(id) |> unique()
sri_data_ids
#>  [1] "e06521cc" "4cda994c" "a8194247" "9e43829a" "5f28e20f" "5024b37e"
#>  [7] "3c2a132a" "ff8105e4" "f207fd6b" "b1c46113" "6ff93997" "1193992a"
#> [13] "f46f88c5" "b5b1d89e" "aeab920e" "a4899cae" "a027a958" "8369d5e8"
#> [19] "619b1d43" "2dec1b47" "1c559275" "d3350e10" "c826aad6" "b21511da"
#> [25] "9e858cd8" "6f89b17f" "5af4817a" "4c84ce02" "36566db2" "3196960f"
#> [31] "1d215e39" "18d5bdc1" "16791cc5" "f1e30a89" "db850249" "a053ce64"
#> [37] "7cbf7851" "4b44f0ff" "225f2a39" "0d483b7f"

2.12 Compute Sleep Proportions

2.12.1 Delete Old Sleep Proportion Files

Code
old_files <-
  dir_ls(
    path = sleep_prop_dir,
    type = "file",
    regexp = "sleep-prop"
  )

if (length(old_files) > 0) file_delete(old_files)

2.12.2 Get Paths to Tidy and Transformed Actigraphy Data Files

Code
files <-
  dir_ls(
    path = actigraphy_processed_data_dir,
    type = "file",
    regexp = "_actigraphy-processed-data.rds$"
  )

2.12.3 Import Tidy and Transformed Actigraphy Data, Compute Sleep Proportions and Write Sleep Proportion Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_actigraphy-processed-data.rds")

  if (!i_id %in% sri_data_ids) next

  i_data <-
    i |>
    read_rds() |>
    as_tibble() |>
    group_by(ga_week) |>
    mutate(
      ga_week_available_days = (n() * 60) / as.numeric(ddays(1))
    ) |>
    ungroup(ga_week) |>
    filter(ga_week_available_days >= min_days) |>
    as_tsibble(index = timestamp)

  if (nrow(i_data) == 0) {
    cli_alert_info(
      paste0(
        "{.strong {col_red(i_id)}} was dropped because it has no gestational ",
        "week with at least {min_days} days of data."
      )
    )

    cli_progress_update()

    next
  }

  i_sleep_prop_data <- tibble()

  for (j in unique(i_data$ga_week)) {
    j_data <- i_data |> filter(ga_week == j)

    i_sleep_prop_data <-
      i_sleep_prop_data |>
      bind_rows(
        j_data |>
          state_prop(
            state_values = 1 # Just "Sleeping", not "Resting" (2)
          ) |>
          mutate(
            id = i_id,
            ga_week = j
          ) |>
          as_tibble()
      )
  }

  i_sleep_prop_data |>
    write_rds(
      path(
        sleep_prop_dir,
        paste0(i_id, "_sleep-prop-data.rds"))
      )

  cli_progress_update()
}

cli_progress_done()

2.12.4 Group and Transform Sleep Proportions Data

2.12.4.1 Get Paths to Sleep Proportions Data Files

Code
files <-
  dir_ls(
    path = sleep_prop_dir,
    type = "file",
    regexp = "_sleep-prop-data.rds$"
  )

2.12.4.2 Group and Transform Sleep Proportions Data

Code
sleep_prop_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  sleep_prop_data <-
    i |>
    read_rds() |>
    mutate(
      time =
        time |>
        as.numeric() |>
        divide_by(60) |>
        as.integer() |>
        dminutes() |>
        as.numeric() |>
        as_hms()
    ) |>
    bind_rows(sleep_prop_data)

  cli_progress_update()
}

cli_progress_done()
Code
sleep_prop_data <-
  sleep_prop_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
sleep_prop_data |> glimpse()
#> Rows: 185,760
#> Columns: 5
#> $ time    <time> 00:00:00, 00:01:00, 00:02:00, 00:03:00, 00:04:00, 00:05:00…
#> $ state   <list> <Other, Other, State, Other, State, State, Other>, <Other,…
#> $ prop    <dbl> 0.4285714286, 0.4285714286, 0.4285714286, 0.4285714286, 0.5…
#> $ id      <chr> "e06521cc", "e06521cc", "e06521cc", "e06521cc", "e06521cc",…
#> $ ga_week <fct> 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33,…

2.12.5 Compute Sleep Proportions Mean Data by Gestational Week

2.12.5.1 Get Paths to Sleep Proportions Data Files

Code
files <-
  dir_ls(
    path = sleep_prop_dir,
    type = "file",
    regexp = "_sleep-prop-data.rds$"
  )

2.12.5.2 Summarize Sleep Proportions by Gestational Week and Write Sleep Proportions Mean Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_sleep-prop-data.rds")

  i_data <-
    i |>
    read_rds() |>
    as_tibble() |>
    group_by(ga_week) |>
    summarize(prop = mean(prop, na.rm = TRUE)) |>
    mutate(id = i_id) |>
    dplyr::select(id, ga_week, prop)

  i_data |>
    write_rds(
      path(
        sleep_prop_dir,
        paste0(i_id, "_sleep-prop-mean-data.rds")
      )
    )

  cli_progress_update()
}

cli_progress_done()

2.12.6 Group and Transform Sleep Proportions Mean Data

2.12.6.1 Get Paths to Sleep Proportions Mean Data Files

Code
files <-
  dir_ls(
    path = sleep_prop_dir,
    type = "file",
    regexp = "_sleep-prop-mean-data.rds$"
  )

2.12.6.2 Import and Group Sleep Proportions Mean Data

Code
sleep_prop_mean_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  sleep_prop_mean_data <-
    i |>
    read_rds(i) |>
    bind_rows(sleep_prop_mean_data)

  cli_progress_update()
}

cli_progress_done()

2.12.6.3 Transform Sleep Proportions Mean Data

Code
sleep_prop_mean_data <-
  sleep_prop_mean_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
sleep_prop_mean_data |> glimpse()
#> Rows: 129
#> Columns: 3
#> $ id      <chr> "e06521cc", "4cda994c", "e06521cc", "a8194247", "9e43829a",…
#> $ ga_week <fct> 33, 33, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35,…
#> $ prop    <dbl> 0.3103174603, 0.2721891534, 0.3164682540, 0.2696428571, 0.3…

2.13 Compute Illuminance Levels

2.13.1 Delete Old SRI Data Files

Code
old_files <-
  dir_ls(
    path = lux_data_dir,
    type = "file",
    regexp = "sri"
  )

if (length(old_files) > 0) file_delete(old_files)

2.13.2 Get Paths to Tidy and Transformed Actigraphy Data Files

Code
files <-
  dir_ls(
    path = actigraphy_processed_data_dir,
    type = "file",
    regexp = "_actigraphy-processed-data.rds$"
  )

2.13.3 Import Tidy and Transformed Actigraphy Data, Compute Illuminance Levels and Write Illuminance Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_actigraphy-processed-data.rds")

  if (!i_id %in% sri_data_ids) next

  i_data <-
    i |>
    read_rds() |>
    as_tibble()

  i_lux_data <- tibble()

  for (j in unique(i_data$ga_week)) {
    j_lux_data <-
      i_data |>
      filter(ga_week == j) |>
      dplyr::select(timestamp, light) |>
      mutate(time = as_hms(timestamp)) |>
      group_by(time) |>
      summarize(lux = mean(light, na.rm = TRUE)) |>
      ungroup() |>
      mutate(
        id = i_id,
        ga_week = j
      )

    i_lux_data <-
      i_lux_data |>
      bind_rows(j_lux_data)
  }

  i_lux_data |>
    write_rds(
      path(
        lux_data_dir,
        paste0(i_id, "_lux-data.rds"))
    )

  cli_progress_update()
}

cli_progress_done()

2.13.4 Group and Transform Illuminance Data

2.13.4.1 Get Paths to Illuminance Data Files

Code
files <-
  dir_ls(
    path = lux_data_dir,
    type = "file",
    regexp = "_lux-data.rds$"
  )

2.13.4.2 Group and Transform Illuminance Data

Code
lux_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  lux_data <-
    i |>
    read_rds() |>
    mutate(
      time =
        time |>
        as.numeric() |>
        divide_by(60) |>
        as.integer() |>
        dminutes() |>
        as.numeric() |>
        as_hms()
    ) |>
    bind_rows(lux_data)

  cli_progress_update()
}

cli_progress_done()
Code
lux_data <-
  lux_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
lux_data |> glimpse()
#> Rows: 237,329
#> Columns: 4
#> $ time    <time> 00:00:00, 00:01:00, 00:02:00, 00:03:00, 00:04:00, 00:05:00…
#> $ lux     <dbl> 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,…
#> $ id      <chr> "e06521cc", "e06521cc", "e06521cc", "e06521cc", "e06521cc",…
#> $ ga_week <fct> 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32,…

2.13.5 Compute Illuminance Mean Data by Gestational Week

2.13.5.1 Get Paths to Illuminance Data Files

Code
files <-
  dir_ls(
    path = lux_data_dir,
    type = "file",
    regexp = "_lux-data.rds$"
  )

2.13.5.2 Summarize Illuminance by Gestational Week and Write Illuminance Mean Data Files

Code
cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  i_id <- str_remove(basename(i), "_lux-data.rds")

  i_data <-
    i |>
    read_rds() |>
    as_tibble() |>
    group_by(ga_week) |>
    summarize(lux = mean(lux, na.rm = TRUE)) |>
    mutate(id = i_id) |>
    dplyr::select(id, ga_week, lux)

  i_data |>
    write_rds(
      path(
        lux_data_dir,
        paste0(i_id, "_lux-mean-data.rds")
      )
    )

  cli_progress_update()
}

cli_progress_done()

2.13.6 Group and Transform Illuminance Mean Data

2.13.6.1 Get Paths to Illuminance Mean Data Files

Code
files <-
  dir_ls(
    path = lux_data_dir,
    type = "file",
    regexp = "_lux-mean-data.rds$"
  )

2.13.6.2 Import and Group Illuminance Mean Data

Code
lux_mean_data <- tibble()

cli_progress_bar(total = length(files), clear = FALSE)

for (i in files) {
  lux_mean_data <-
    i |>
    read_rds(i) |>
    bind_rows(lux_mean_data)

  cli_progress_update()
}

cli_progress_done()

2.13.6.3 Transform Illuminance Mean Data

Code
lux_mean_data <-
  lux_mean_data |>
  mutate(ga_week = as.factor(ga_week)) |>
  arrange(ga_week)
lux_mean_data |> glimpse()
#> Rows: 170
#> Columns: 3
#> $ id      <chr> "e06521cc", "e06521cc", "a8194247", "5024b37e", "4cda994c",…
#> $ ga_week <fct> 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 34, 34, 34, 34,…
#> $ lux     <dbl> 381.44080671, 211.18389400, 145.59865162, 54.87876126, 123.…

2.14 Categorize Participantes by SRI Levels

Code
sri_data_by_id <-
  sri_data |>
  summarize(
    sri = mean(sri, na.rm = TRUE),
    .by = id
  )
Code
sri_quintiles <-
  sri_data_by_id |>
  pull(sri) |>
  quantile(probs = c(0, 0.2, 0.4, 0.6, 0.8, 1))
sri_quintiles
#>          0%         20%         40%         60%         80%        100% 
#> 44.07514451 53.52458860 61.75108230 69.86563612 74.73463204 83.22344322
Code
sri_top_quintile <- sri_quintiles[5]
Code
sri_bottom_quintile <- sri_quintiles[2]
Code
regular_sri_ids <-
  sri_data_by_id |>
  filter(sri >= sri_top_quintile) |>
  pull(id)
regular_sri_ids
#> [1] "3c2a132a" "6ff93997" "aeab920e" "2dec1b47" "c826aad6" "b21511da"
#> [7] "5af4817a" "3196960f"
Code
others_sri_ids <-
  sri_data_by_id |>
  filter(
    sri > sri_bottom_quintile &
    sri < sri_top_quintile
  ) |>
  pull(id)
others_sri_ids
#>  [1] "e06521cc" "4cda994c" "a8194247" "9e43829a" "5f28e20f" "5024b37e"
#>  [7] "ff8105e4" "f207fd6b" "b1c46113" "f46f88c5" "b5b1d89e" "a4899cae"
#> [13] "a027a958" "8369d5e8" "1c559275" "d3350e10" "9e858cd8" "6f89b17f"
#> [19] "36566db2" "1d215e39" "18d5bdc1" "7cbf7851" "225f2a39" "0d483b7f"
Code
irregular_sri_ids <-
  sri_data_by_id |>
  filter(sri <= sri_bottom_quintile) |>
  pull(id)
irregular_sri_ids
#> [1] "1193992a" "619b1d43" "4c84ce02" "16791cc5" "f1e30a89" "db850249"
#> [7] "a053ce64" "4b44f0ff"

2.15 Filter Datasets by Gestational Age Range

2.15.1 Define Gestational Age Range

Code
sri_mean_data |>
  summarize(
    n = unique(id) |> length(),
    .by = ga_week
  ) |>
  print(n = Inf)
#> # A tibble: 8 × 2
#>   ga_week     n
#>   <fct>   <int>
#> 1 33          2
#> 2 34          7
#> 3 35         11
#> 4 36         20
#> 5 37         29
#> 6 38         26
#> 7 39         18
#> 8 40          5
Code
ga_range <- c(36, 39) # Gestational age range

2.15.2 Filter Datasets

Code
datasets <- c(
  "sri_data",
  "sri_mean_data",
  "sleep_prop_data",
  "sleep_prop_mean_data",
  "lux_data",
  "lux_mean_data"
)

for (i in datasets) {
  assign(i, get(i) |>
    filter(
      between(
        ga_week |>
          as.character() |>
          as.numeric(),
        ga_range[1],
        ga_range[2]
      )
    )
  )
}

2.16 Explore Data

2.16.1 SRI

2.16.1.1 Mean SRI

2.16.1.1.1 By Pregnant Women
Code
sri_mean_data_by_id <-
  sri_data |>
  summarize(
    sri = mean(sri, na.rm = TRUE),
    .by = id
  )

sri_mean_data_by_id
sri_mean_data_by_id |> glimpse()
#> Rows: 39
#> Columns: 2
#> $ id  <chr> "ff8105e4", "f46f88c5", "f207fd6b", "e06521cc", "b5b1d89e", "ae…
#> $ sri <dbl> 74.96453373, 53.61662257, 80.27673720, 72.84796592, 73.72551478…
2.16.1.1.2 Global SRI
2.16.1.1.2.1 Using All the Data

Mean SRI considering all the SRI data points across all pregnant women.

Code
sri_data |> pull(sri) |> mean(na.rm = TRUE)
#> [1] 65.88659463
2.16.1.1.2.2 Using the Means by Pregnant Women
Code
sri_data |>
  summarize(
    sri = mean(sri, na.rm = TRUE),
    .by = id
  ) |>
  pull(sri) |>
  mean(na.rm = TRUE)
#> [1] 65.77606715
2.16.1.1.3 Mean SRI by Gestational Age Week
Code
sri_data |>
  summarize(
    sri = mean(sri, na.rm = TRUE),
    .by = ga_week
  )
2.16.1.1.4 Histogram
Code
bin_width <- 5
Code
sri_mean_data_by_id |>
  ggplot(aes(x = sri)) +
  geom_histogram(
    binwidth = bin_width,
    color = "black",
    fill = "white"
  ) +
  geom_density(
    aes(y = ..count.. * bin_width),
    color = "red",
    size = 1
  ) +
  labs(
    x = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
    y = "Frequência", # Frequency
  )
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(count)` instead.

2.16.1.2 By Gestational Age Week

2.16.1.2.1 Scatter Plot with Linear Trend
Code
sri_mean_data |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  ggplot(aes(x = ga_week, y = sri)) +
  geom_point() +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    color = "red"
  ) +
  labs(
    x = "Semana Gestacional", # Gestational Age Week
    y = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
  )

Code
  theme(legend.position = "none")
#> <theme> List of 1
#>  $ legend.position: chr "none"
#>  @ complete: logi FALSE
#>  @ validate: logi TRUE
2.16.1.2.2 Boxplots
Code
sri_mean_data |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  mutate(ga_week = factor(ga_week)) |>
  ggplot(aes(x = sri, y = ga_week)) +
  geom_boxplot() +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 4,
    color = "red",
    shape = 18
  ) +
  geom_point(
    position = position_jitter(height = 0.1, width = 0),
    size = 2,
    alpha = 0.25,
    color = "blue"
  ) +
  coord_flip() +
  labs(
    x = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
    y = "Semana Gestacional", # Gestational Age Week
  ) +
  scale_color_viridis_d() +
  theme(legend.position = "none")

2.16.1.3 By Gestational Age Week & Time of Day

2.16.1.3.1 Global Trend
Code
sri_plot_global <-
  sri_data |>
  ggplot(aes(x = time, y = sri, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = sri),
    color = "#000000",
    fill = lighten("#000000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sri_plot_global

2.16.1.3.2 Regular Trend
Code
sri_plot_regular <-
  sri_data |>
  filter(id %in% regular_sri_ids) |>
  ggplot(aes(x = time, y = sri, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = sri),
    color = "#0000FF",
    fill = lighten("#0000FF", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sri_plot_regular

2.16.1.3.3 “Others” Trend
Code
sri_plot_others <-
  sri_data |>
  filter(id %in% others_sri_ids) |>
  ggplot(aes(x = time, y = sri, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = sri),
    color = "#5E5E5E",
    fill = lighten("#5E5E5E", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sri_plot_others

2.16.1.3.4 Irregular Trend
Code
sri_plot_irregular <-
  sri_data |>
  filter(id %in% irregular_sri_ids) |>
  ggplot(aes(x = time, y = sri, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = sri),
    color = "#FF0000",
    fill = lighten("#FF0000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Regularidade do Sono (SRI)", # Sleep Regularity Index (SRI)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sri_plot_irregular

2.16.2 Sleep Proportions

2.16.2.1 Mean Sleep Proportions by Pregnant Women

Code
sleep_prop_mean_by_id_data <-
  sleep_prop_data |>
  summarize(
    prop = mean(prop, na.rm = TRUE),
    .by = id
  ) |>
  mutate(
    prop_hours =
      prop |>
      multiply_by(24) |>
      multiply_by(60) |>
      multiply_by(60) |>
      as_hms()
  )

sleep_prop_mean_by_id_data

2.16.2.2 Mean Global Sleep Proportion

2.16.2.2.1 Using All Data

Mean sleep proportion considering all the proportion data points across all pregnant women.

Code
sleep_prop_data |>
  pull(prop) |>
  mean(na.rm = TRUE) |>
  multiply_by(24) |>
  multiply_by(60) |>
  multiply_by(60) |>
  as_hms()
#> 06:51:45.177531
2.16.2.2.2 Using the Means by Pregnant Women
Code
sleep_prop_data |>
  summarize(
    prop = mean(prop, na.rm = TRUE),
    .by = id
  ) |>
  pull(prop) |>
  mean(na.rm = TRUE) |>
  multiply_by(24) |>
  multiply_by(60) |>
  multiply_by(60) |>
  as_hms()
#> 06:52:19.960714

2.16.2.3 Mean Sleep Proportion by Gestational Age Week

Code
sleep_prop_data |>
  summarize(
    prop = mean(prop, na.rm = TRUE),
    .by = ga_week
  ) |>
  mutate(
    prop_hours =
      prop |>
      multiply_by(24) |>
      multiply_by(60) |>
      multiply_by(60) |>
      as_hms()
  )

2.16.2.4 By Gestational Age Week

2.16.2.4.1 Scatter Plot with Linear Trend
Code
sleep_prop_mean_data |>
  mutate(
    prop_hours =
      prop |>
      multiply_by(24) |>
      multiply_by(60) |>
      multiply_by(60) |>
      as_hms()
  ) |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  ggplot(aes(x = ga_week, y = prop_hours)) +
  geom_point() +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    color = "red"
  ) +
  labs(
    x = "Semana Gestacional", # Gestational Age Week
    y = "Tempo Total de Sono (TTS)", # Total Sleep Time
  ) +
  theme(legend.position = "none")

2.16.2.4.2 Boxplots
Code
sleep_prop_mean_data |>
  mutate(
    prop_hours =
      prop |>
      multiply_by(24) |>
      multiply_by(60) |>
      multiply_by(60) |>
      as_hms()
  ) |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  mutate(ga_week = factor(ga_week)) |>
  ggplot(aes(x = prop_hours, y = ga_week)) +
  geom_boxplot() +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 4,
    color = "red",
    shape = 18
  ) +
  geom_point(
    position = position_jitter(height = 0.1, width = 0),
    size = 2,
    alpha = 0.25,
    color = "blue"
  ) +
  coord_flip() +
  labs(
    x = "Tempo Total de Sono (TTS)", # Total Sleep Time
    y = "Semana Gestacional", # Gestational Age Week
  ) +
  scale_color_viridis_d() +
  theme(legend.position = "none")

2.16.2.5 By Gestational Age Week & Time of Day

2.16.2.5.1 Global Trend
Code
sleep_prop_plot_global <-
  sleep_prop_data |>
  mutate(per = prop * 100) |>
  ggplot(aes(x = time, y = per, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = per),
    color = "#000000",
    fill = lighten("#000000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Percentual de Tempo Dormindo (%)",
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sleep_prop_plot_global

2.16.2.5.2 Regular Trend
Code
sleep_prop_plot_regular <-
  sleep_prop_data |>
  mutate(per = prop * 100) |>
  filter(id %in% regular_sri_ids) |>
  ggplot(aes(x = time, y = per, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = per),
    color = "#0000FF",
    fill = lighten("#0000FF", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Percentual de Tempo Dormindo (%)", # Sleep Proportion (%)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sleep_prop_plot_regular

2.16.2.5.3 “Others” Trend
Code
sleep_prop_plot_others <-
  sleep_prop_data |>
  mutate(per = prop * 100) |>
  filter(id %in% others_sri_ids) |>
  ggplot(aes(x = time, y = per, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = per),
    color = "#5E5E5E",
    fill = lighten("#5E5E5E", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Percentual de Tempo Dormindo (%)", # Sleep Proportion (%)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sleep_prop_plot_others

2.16.2.5.4 Irregular Trend
Code
sleep_prop_plot_irregular <-
  sleep_prop_data |>
  mutate(per = prop * 100) |>
  filter(id %in% irregular_sri_ids) |>
  ggplot(aes(x = time, y = per, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = per),
    color = "#FF0000",
    fill = lighten("#FF0000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Percentual de Tempo Dormindo (%)", # Sleep Proportion (%)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(legend.position = "none")

sleep_prop_plot_irregular

2.16.3 Illuminance Levels

2.16.3.1 Statistics

2.16.3.1.1 Mean Illuminance by Pregnant Women
Code
lux_data |>
  summarize(
    lux = mean(lux, na.rm = TRUE),
    .by = id
  )
2.16.3.1.2 Mean Global Illuminance
2.16.3.1.2.1 Using All Data

Mean SRI considering all the SRI data points across all pregnant women.

Code
lux_data |> pull(lux) |> mean(na.rm = TRUE)
#> [1] 137.3573708
2.16.3.1.2.2 Using the Means by Pregnant Women
Code
lux_data |>
  summarize(
    lux = mean(lux, na.rm = TRUE),
    .by = id
  ) |>
  pull(lux) |>
  mean(na.rm = TRUE)
#> [1] 135.9001966
2.16.3.1.3 Mean Illuminance by Gestational Age Week
Code
lux_data |>
  summarize(
    lux = mean(lux, na.rm = TRUE),
    .by = ga_week
  )

2.16.3.2 By Gestational Age Week

2.16.3.2.1 Scatter Plot with Linear Trend
Code
lux_mean_data |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  ggplot(aes(x = ga_week, y = lux)) +
  geom_point() +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    color = "red"
  ) +
  labs(
    x = "Semana Gestacional", # Gestational Age Week
    y = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
  ) +
  theme(legend.position = "none")

2.16.3.2.2 Boxplots
Code
lux_mean_data |>
  mutate(ga_week = ga_week |> as.character() |> as.numeric()) |>
  mutate(ga_week = factor(ga_week)) |>
  ggplot(aes(x = lux, y = ga_week)) +
    geom_boxplot() +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 4,
    color = "red",
    shape = 18
  ) +
  geom_point(
    position = position_jitter(height = 0.1, width = 0),
    size = 2,
    alpha = 0.25,
    color = "blue"
  ) +
  coord_flip() +
  # scale_y_discrete(limits = rev) +
  labs(
    x = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
    y = "Semana Gestacional", # Gestational Age Week
  ) +
  scale_color_viridis_d() +
  theme(legend.position = "none")

2.16.3.3 By Gestational Age Week & Time of Day

2.16.3.3.1 Global Trend
Code
lux_plot_global <-
  lux_data |>
  ggplot(aes(x = time, y = lux, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = lux),
    color = "#000000",
    fill = lighten("#000000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(
    limits = c(0, 1200),
    breaks = seq(0, 1200, by = 300)
  ) +
  theme(legend.position = "none")

lux_plot_global

2.16.3.3.2 Regular Trend
Code
lux_plot_regular <-
  lux_data |>
  filter(id %in% regular_sri_ids) |>
  ggplot(aes(x = time, y = lux, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = lux),
    color = "#0000FF",
    fill = lighten("#0000FF", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(
    limits = c(0, 1200),
    breaks = seq(0, 1200, by = 300)
  ) +
  theme(legend.position = "none")

lux_plot_regular

2.16.3.3.3 “Others” Trend
Code
lux_plot_others <-
  lux_data |>
  filter(id %in% others_sri_ids) |>
  ggplot(aes(x = time, y = lux, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = lux),
    color = "#5E5E5E",
    fill = lighten("#5E5E5E", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(
    limits = c(0, 1200),
    breaks = seq(0, 1200, by = 300)
  ) +
  theme(legend.position = "none")

lux_plot_others

2.16.3.3.4 Irregular Trend
Code
lux_plot_irregular <-
  lux_data |>
  filter(id %in% irregular_sri_ids) |>
  ggplot(aes(x = time, y = lux, color = id)) +
  geom_smooth(
    se = FALSE
  ) +
  geom_smooth(
    aes(x = time, y = lux),
    color = "#FF0000",
    fill = lighten("#FF0000", 0.5),
    linewidth = 2,
    se = TRUE
  ) +
  scale_color_grey(
    start = 0.7,
    end = 0.9
  ) +
  labs(
    x = "Hora do Dia", # Time of day (Hour)
    y = "Índice de Iluminância (Lux)", # Illuminance Index (Lux)
  ) +
  scale_x_time(
    breaks = breaks_width("6 hours"),
    labels = label_time("%-H") # Use "%#H" for Windows
  ) +
  scale_y_continuous(
    limits = c(0, 1200),
    breaks = seq(0, 1200, by = 300)
  ) +
  theme(legend.position = "none")

lux_plot_irregular

2.16.4 SRI versus Illuminance

Code
wrap_plots(
  A = sri_plot_regular,
  X = ggplot() + theme_void(),
  C = lux_plot_regular,
  B = sri_plot_irregular,
  Y = ggplot() + theme_void(),
  D = lux_plot_irregular,
  ncol = 3,
  nrow = 2,
  widths = c(0.4875, 0.025, 0.4875),
  axes = "collect",
  tag_level = "new"
)

2.16.5 Sleep Proportion versus Illuminance

Code
wrap_plots(
  A = sleep_prop_plot_regular,
  X = ggplot() + theme_void(),
  C = lux_plot_regular,
  B = sleep_prop_plot_irregular,
  Y = ggplot() + theme_void(),
  D = lux_plot_irregular,
  ncol = 3,
  nrow = 2,
  widths = c(0.4875, 0.025, 0.4875),
  axes = "collect",
  tag_level = "new"
)

2.16.6 SRI versus Sleep Duration

Code
plot_data <-
  sri_mean_data_by_id |>
  left_join(
    sleep_prop_mean_by_id_data,
    by = "id"
  ) |>
  mutate(
    group = case_when(
      id %in% regular_sri_ids ~ "Regular",
      id %in% irregular_sri_ids ~ "Irregular",
      TRUE ~ "Outros"
    )
  )
Code
plot_data |> glimpse()
#> Rows: 39
#> Columns: 5
#> $ id         <chr> "ff8105e4", "f46f88c5", "f207fd6b", "e06521cc", "b5b1d89…
#> $ sri        <dbl> 74.96453373, 53.61662257, 80.27673720, 72.84796592, 73.7…
#> $ prop       <dbl> 0.2855902778, 0.2251157407, 0.3055720899, 0.3017609127, …
#> $ prop_hours <time> 06:51:15.000000, 05:24:10.000000, 07:20:01.428571, 07:1…
#> $ group      <chr> "Outros", "Outros", "Outros", "Outros", "Outros", "Regul…
Code
plot_data |>
  summarize(
    n = n(),
    x = mean(sri, na.rm = TRUE),
    y =
      prop_hours |>
      as.numeric() |>
      mean(na.rm = TRUE) |>
      as_hms(),
    x_se =
      sri |>
      sd(na.rm = TRUE) |>
      divide_by(sqrt(n())),
    y_se =
      prop_hours |>
      as.numeric() |>
      sd(na.rm = TRUE) |>
      divide_by(sqrt(n())) |>
      as_hms()
  )
Code
plot_data_se <-
  plot_data |>
  summarize(
    n = n(),
    x = mean(sri, na.rm = TRUE),
    y =
      prop_hours |>
      as.numeric() |>
      mean(na.rm = TRUE) |>
      as_hms(),
    x_se =
      sri |>
      sd(na.rm = TRUE) |>
      divide_by(sqrt(n())),
    y_se =
      prop_hours |>
      as.numeric() |>
      sd(na.rm = TRUE) |>
      divide_by(sqrt(n())) |>
      as_hms(),
    .by = group
  )
Code
plot_data_se
Code
plot_data_se <- plot_data_se |> filter(!group == "Outros")
Code
plot_data |>
  ggplot(
    aes(x = sri, y = prop_hours, color = group, shape = group)
  ) +
  geom_point(
    size = 4
  ) +
  geom_crossbar(
    data = plot_data_se,
    aes(x = x, y = y, ymin = y - y_se, ymax = y + y_se, color = group),
    width = 0,
    inherit.aes = FALSE
  ) +
  geom_errorbar(
    data = plot_data_se,
    aes(x = x, y = y, xmin = x - x_se, xmax = x + x_se, color = group),
    height = 0,
    inherit.aes = FALSE
  ) +
  scale_shape_manual(
    values = c(
      "Regular" = 15,
      "Irregular" = 16,
      "Outros" = 6
    )
  ) +
  scale_color_manual(
    values = c(
      "Regular" = "blue",
      "Irregular" = "red",
      "Outros" = "black"
    )
  ) +
  scale_y_time(
    limits = c(as_hms("04:30:00"), as_hms("08:45:00")),
    breaks = "1 hour",
    labels = function(x) {
      x |>
        as.POSIXct() |>
        format("%H:%M")
    }
  ) +
  scale_x_continuous(
    limits = c(42.5, 85),
    breaks = seq(40, 100, by = 10)
  ) +
  labs(
    # Sleep Regularity Index (SRI)
    x = "Índice de Regularidade do Sono (SRI)",
    # Average Daily Total Sleep Time (TST)
    y = "Média do Tempo Total de Sono (TTS)",
    color = NULL,
    shape = NULL
  ) +
  guides(
    color = guide_legend(override.aes = list(linetype = 0))
  )

2.17 Differences Test

2.18 Compute an A Priori Power Analysis

Code
set.seed(2025)
Code
n_clusters <- 66 # Number of participants
cluster_size <- 2 # Average number of weeks per participant
alpha <- 0.3 # Intra-class correlation
beta <- - 2 # Effect size (per week)
sigma <- 5 # Standard deviation of the outcome
type_one_error <- 0.05
n_sim <- 1000
rejects <- 0
Code
cli_progress_bar(total = n_sim, clear = FALSE)

for (i in seq_len(n_sim)){
  sim_data <- tibble(
    id = seq_len(n_clusters) |> rep(each = cluster_size),
    x = seq_len(cluster_size) |> rep(n_clusters)
  )

  correlated_errors <-
    rep(
      rnorm(n_clusters, 0, sigma * sqrt(alpha)),
      each = cluster_size
    ) |>
    add(
      rnorm(n_clusters * cluster_size, 0, sigma * sqrt(1 - alpha))
    )

  sim_data <- sim_data |> mutate(y = beta * x + correlated_errors)

  fit <- geeglm(
    y ~ x,
    family = gaussian,
    id = id,
    data = sim_data,
    corstr = "exchangeable"
  )

  p_value <- summary(fit)$coefficients[2,4]

  if (p_value < type_one_error) rejects <- rejects + 1

  cli_progress_update()
}
#> ■■■■■■■                           19% | ETA:  4s
#> ■■■■■■■■■■■■■                     39% | ETA:  3s
#> ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  2s
#> ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
#> 

cli_progress_done()

power_estimate <- rejects / n_sim
power_estimate
#> [1] 0.803

2.19 Model Data

2.19.1 Prepare Model Data

Code
model_data <-
  sri_mean_data |>
  mutate(
    ga_week =
      ga_week |>
      as.character() |>
      as.numeric()
  ) |>
  mutate(
    ga_week_centered =
      ga_week |>
      scale(center = TRUE, scale = FALSE) # Mean deviations
  ) |>
  group_by(ga_week) |>
  mutate(n = n()) |>
  ungroup() |>
  left_join(
    field_form_data |> dplyr::select(-timestamp),
    by = "id"
  )

2.19.2 Inspect Model Data

model_data |> nrow()
#> [1] 93
model_data |>
  pull(id) |>
  unique() |>
  length()
#> [1] 39
model_data |>
  pull(ga_week) |>
  unique() |>
  length()
#> [1] 4
model_data |>
  summarize(
    n = n(),
    .by = id
  ) |>
    arrange(n) |>
    print(n = Inf)
#> # A tibble: 39 × 2
#>    id           n
#>    <chr>    <int>
#>  1 9e43829a     1
#>  2 3c2a132a     1
#>  3 1193992a     1
#>  4 c826aad6     1
#>  5 b21511da     1
#>  6 36566db2     1
#>  7 db850249     1
#>  8 7cbf7851     1
#>  9 4b44f0ff     1
#> 10 225f2a39     1
#> 11 f207fd6b     2
#> 12 aeab920e     2
#> 13 a8194247     2
#> 14 8369d5e8     2
#> 15 5f28e20f     2
#> 16 d3350e10     2
#> 17 16791cc5     2
#> 18 f1e30a89     2
#> 19 a053ce64     2
#> 20 0d483b7f     2
#> 21 f46f88c5     3
#> 22 b5b1d89e     3
#> 23 6ff93997     3
#> 24 5024b37e     3
#> 25 4cda994c     3
#> 26 2dec1b47     3
#> 27 9e858cd8     3
#> 28 6f89b17f     3
#> 29 5af4817a     3
#> 30 4c84ce02     3
#> 31 3196960f     3
#> 32 1d215e39     3
#> 33 18d5bdc1     3
#> 34 ff8105e4     4
#> 35 e06521cc     4
#> 36 a4899cae     4
#> 37 a027a958     4
#> 38 619b1d43     4
#> 39 1c559275     4
model_data |> glimpse()
#> Rows: 93
#> Columns: 23
#> $ id               <chr> "ff8105e4", "f46f88c5", "f207fd6b", "e06521cc", "b…
#> $ ga_week          <dbl> 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36…
#> $ sri              <dbl> 74.36507937, 73.05555556, 83.63642535, 67.72067154…
#> $ ga_week_centered <dbl[,1]> <matrix[25 x 1]>
#> $ n                <int> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20…
#> $ birth_date       <date> 1992-06-23, 2000-12-22, 1990-04-12, 1986-08-16, 20…
#> $ age              <dbl> 32.85277778, 24.19444444, 35.03055556, 38.7055555…
#> $ weight_before    <dbl> 68, 62, 63, 65, 85, 74, 67, 52, 64, 89, 77, 55, 80…
#> $ weight_current   <dbl> 80.0, 72.0, 76.0, 77.0, 78.0, 87.0, 81.7, 65.0, 71…
#> $ height           <dbl> 161, 156, 158, 170, 169, 170, 168, 165, 150, 169, …
#> $ bmi_before       <dbl> 26.23355580, 25.47666009, 25.23634033, 22.49134948…
#> $ bmi_current      <dbl> 30.86300683, 29.58579882, 30.44383913, 26.64359862…
#> $ family_income    <dbl> 13000.0, 2600.0, 4000.0, 4000.0, NA, 2500.0, 6.8, …
#> $ gestations       <dbl> 2, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 3, 1,…
#> $ study            <lgl> TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
#> $ work             <lgl> TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, T…
#> $ health_plan      <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
#> $ solo             <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FA…
#> $ ethnicity        <ord> Preta, Parda, Parda, Branca, Parda, Parda, Branca,…
#> $ exercise         <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, …
#> $ education        <ord> Ensino superior completo, Ensino superior completo…
#> $ deliveries       <dbl> 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,…
#> $ abortions        <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,…

2.19.3 Model Data with Generalized Estimating Equations (GEE)

2.19.3.1 Prepare Model Data for GEE

Code
model_data_gee <-
  model_data |>
  mutate(
    id =
      id |>
      as.factor() |>
      as.numeric()
  ) |>
  arrange(id, ga_week)
model_data_gee |>
  summarize(
    n = n(),
    .by = id
  ) |>
    arrange(n) |>
    print(n = Inf)
#> # A tibble: 39 × 2
#>       id     n
#>    <dbl> <int>
#>  1     2     1
#>  2     7     1
#>  3    10     1
#>  4    11     1
#>  5    12     1
#>  6    21     1
#>  7    23     1
#>  8    30     1
#>  9    32     1
#> 10    34     1
#> 11     1     2
#> 12     3     2
#> 13    17     2
#> 14    22     2
#> 15    26     2
#> 16    28     2
#> 17    29     2
#> 18    33     2
#> 19    36     2
#> 20    37     2
#> 21     4     3
#> 22     6     3
#> 23     8     3
#> 24     9     3
#> 25    13     3
#> 26    14     3
#> 27    15     3
#> 28    16     3
#> 29    19     3
#> 30    20     3
#> 31    24     3
#> 32    31     3
#> 33    38     3
#> 34     5     4
#> 35    18     4
#> 36    25     4
#> 37    27     4
#> 38    35     4
#> 39    39     4

2.19.3.2 Model Data

model_gee <- geeglm(
  sri ~ ga_week_centered,
  family = gaussian,
  id = id,
  weights = n,
  data = model_data_gee,
  corstr = "exchangeable"
)

2.19.3.3 Inspect Model Results

model_gee
#> 
#> Call:
#> geeglm(formula = sri ~ ga_week_centered, family = gaussian, data = model_data_gee, 
#>     weights = n, id = id, corstr = "exchangeable")
#> 
#> Coefficients:
#>      (Intercept) ga_week_centered 
#>     64.980416414     -2.386303277 
#> 
#> Degrees of Freedom: 93 Total (i.e. Null);  91 Residual
#> 
#> Scale Link:                   identity
#> Estimated Scale Parameters:  [1] 168.6295269
#> 
#> Correlation:  Structure = exchangeable    Link = identity 
#> Estimated Correlation Parameters:
#>        alpha 
#> 0.4109994934 
#> 
#> Number of clusters:   39   Maximum cluster size: 4
summary(model_gee)
#> 
#> Call:
#> geeglm(formula = sri ~ ga_week_centered, family = gaussian, data = model_data_gee, 
#>     weights = n, id = id, corstr = "exchangeable")
#> 
#>  Coefficients:
#>                   Estimate   Std.err       Wald Pr(>|W|)    
#> (Intercept)      64.980416  1.789900 1317.97521  < 2e-16 ***
#> ga_week_centered -2.386303  1.254093    3.62069 0.057065 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Correlation structure = exchangeable 
#> Estimated Scale Parameters:
#> 
#>             Estimate  Std.err
#> (Intercept) 168.6295 24.31045
#>   Link = identity 
#> 
#> Estimated Correlation Parameters:
#>        Estimate   Std.err
#> alpha 0.4109995 0.1265896
#> Number of clusters:   39  Maximum cluster size: 4

2.19.3.4 Transform p-value for an One-Tailed Hypothesis

Code
coefs <- summary(model_gee)$coefficients
beta  <- coefs["ga_week_centered", "Estimate"]
se    <- coefs["ga_week_centered", "Std.err"]

p_one_tailed <- pnorm(beta / se)
p_one_tailed
#> [1] 0.02853256

2.19.4 Model Diagnostics

2.19.4.1 Check Normality of Residuals

Code
tibble(x = residuals(model_gee)[, 1]) |>
  ggplot(aes(sample = x)) +
  stat_qq() +
  stat_qq_line() +
  labs(
    x = "Quantis Teóricos (Normal Padrão)",
    y = "Quantis Amostrais"
  )

tibble(x = residuals(model_gee)[, 1]) |>
  rutils:::normality_summary("x") |>
  print(n = Inf)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
#> # A tibble: 11 × 4
#>    test                     statistic_1 statistic_2  p_value
#>    <chr>                          <dbl>       <dbl>    <dbl>
#>  1 Anderson-Darling               1.71        NA    0.000205
#>  2 Bonett-Seier                  52.9         -1.50 0.134   
#>  3 Cramer-von Mises               0.320       NA    0.000188
#>  4 D'Agostino Omnibus Test        5.72        NA    0.0574  
#>  5 D'Agostino Skewness Test      -2.25        NA    0.0247  
#>  6 D'Agostino Kurtosis Test      -0.818       NA    0.413   
#>  7 Jarque-Bera                    5.69         2    0.0581  
#>  8 Lilliefors (K-S)               0.124       NA    0.00117 
#>  9 Pearson chi-square            29.3         NA    0.00111 
#> 10 Shapiro-Francia                0.954       NA    0.00259 
#> 11 Shapiro-Wilk                   0.957       NA    0.00547

2.19.4.2 Use performance Package to Check the Model

Code
check_model(model_gee)
#> Converting missing values (`NA`) into regular values currently not
#>   possible for variables of class `NULL`.

2.19.5 Check the Standardized Effect by Week

\[ d_{Total} = \cfrac{\beta}{\text{SD}(Y)} \]

Code
beta <- coef(model_gee)["ga_week_centered"]
scale_parameter <- summary(model_gee)$dispersion[[1]]
sd_residuals <- sqrt(scale_parameter)

d_by_week <- beta / sd_residuals
d_by_week
#> ga_week_centered 
#>       -0.1837633
Code
se_beta <- summary(model_gee)$coefficients["ga_week_centered", "Std.err"]

# 1.645 == One-tailed 95% CI z-score
ci_lower <- (beta - (1.645 * se_beta)) / sd_residuals
ci_upper <- (beta + (1.645 * se_beta)) / sd_residuals
c(ci_lower, ci_upper)
#> ga_week_centered ga_week_centered 
#>      -0.34262857      -0.02489806
interpret_cohens_d(d_by_week)
#> ga_week_centered 
#>     "very small" 
#> (Rules: cohen1988)

2.19.6 Check Standardized Effect for 5 Weeks

\[ d_{Total} = \cfrac{\beta \times \Delta\text{Weeks}}{\text{SD}(Y)} \]

Code
d_5_weeks <- d_by_week * 5
d_5_weeks
#> ga_week_centered 
#>       -0.9188166
interpret_cohens_d(d_5_weeks)
#> ga_week_centered 
#>          "large" 
#> (Rules: cohen1988)

2.20 Compute an A Posteriori Power Analysis

Code
set.seed(2025)
Code
n_clusters <- summary(model_gee)$clusz |> length()
cluster_size <- summary(model_gee)$clusz |> mean() |> floor()
alpha <- summary(model_gee)$corr[[1]]
beta <- coefficients(model_gee)[2]
sigma <- summary(model_gee)$dispersion[[1]] |> sqrt()
n_sim <- 1000
rejects <- 0
cli_progress_bar(total = n_sim, clear = FALSE)

for (i in seq_len(n_sim)){
  sim_data <- tibble(
    id = seq_len(n_clusters) |> rep(each = cluster_size),
    x = seq_len(cluster_size) |> rep(n_clusters)
  )

  correlated_errors <-
    rnorm(n = n_clusters, mean = 0, sd = sigma * sqrt(alpha)) |>
    rep(each = cluster_size) |>
    add(
      rnorm(n_clusters * cluster_size, 0, sigma * sqrt(1 - alpha))
    )

  sim_data <-
    sim_data |>
    mutate(
      y = beta * x + correlated_errors
    )

  fit <- geeglm(
    y ~ x,
    family = gaussian,
    id = id,
    data = sim_data,
    corstr = "exchangeable"
  )

  p_value <- summary(fit)$coefficients[2,4]

  if (p_value < 0.05) rejects <- rejects + 1

  cli_progress_update()
}
#> ■■■■■■                            15% | ETA:  6s
#> ■■■■■■■■■■■■■■■■■■                56% | ETA:  3s
#> ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     92% | ETA:  1s
#> ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
#> 

cli_progress_done()

power_estimate <- rejects / n_sim
power_estimate
#> [1] 0.191