Occlusion Coding

library(tidyverse)
library(readxl)
library(gt)

dat <- read_xlsx(
  "/mnt/datasets/SEEDLingS/1_RAWDATA/FACS/FACS_code_inventory.xlsx",
  sheet = 2,
  na = c("", "NA", "nan"),
  .name_repair = janitor::make_clean_names
)

# Reshape wide per-wave columns to long, then back to one row per subject-wave.
recode_var <- function(name) {
  case_when(
    str_detect(name, "^fullfaceavail_") ~ "ffa",
    str_detect(name, "^no_facs_vis")    ~ "vis",   # no underscore before wave
    str_detect(name, "^no_facs_hat")    ~ "hat",   # no underscore before wave
    str_detect(name, "^no_fac_sn_")     ~ "neu",
    str_detect(name, "^facs_tot_")      ~ "auc",
    .default = NA_character_
  )
}

dat_long <-
  dat |>
  select(
    subject = id,
    starts_with("fullfaceavail_"),
    starts_with("no_facs_vis"),
    starts_with("no_facs_hat"),
    starts_with("no_fac_sn_"),
    starts_with("facs_tot_")
  ) |>
  pivot_longer(-subject) |>
  mutate(
    var  = recode_var(name),
    wave = as.integer(str_extract(name, "\\d+(?=m?$)")),  # digits before optional "m"
    value = as.numeric(value)
  )

unmatched <- dat_long |> filter(is.na(var)) |> distinct(name) |> pull(name)
if (length(unmatched) > 0) {
  stop("Unrecognized occlusion columns: ", paste(unmatched, collapse = ", "))
}

dat2 <-
  dat_long |>
  select(-name) |>
  pivot_wider(names_from = var, values_from = value)

n_bad_neu <- sum(dat2$neu < 0, na.rm = TRUE)
cat("Sessions with negative neutral time (set to NA):", n_bad_neu, "\n")

Sessions with negative neutral time (set to NA): 0 

dat2 <- dat2 |> mutate(neu = if_else(neu < 0, NA_real_, neu))

check <-
  dat2 |>
  drop_na(ffa, vis, hat, neu, auc) |>
  mutate(diff = (neu + auc) - (ffa - vis - hat))

cat("Rows checked:", nrow(check),
    "| max |difference| (s):", format(max(abs(check$diff)), digits = 3), "\n")

Rows checked: 136 | max |difference| (s): 1.42e-14 

stopifnot(max(abs(check$diff)) < 1)   # seconds; should be ~0

stopifnot(all(dat2$wave %in% c(6, 9, 12, 15)))

write_rds(dat2, "facs_times.rds")

fig1df <-
  dat2 |>
  pivot_longer(vis:auc, names_to = "type", values_to = "sec") |>
  drop_na() |>
  summarize(.by = c(wave, type), sec = sum(sec)) |>
  mutate(
    wave = factor(wave),
    type = case_when(
      type == "auc" ~ "Any Other AU",
      type == "neu" ~ "No Coded AU",
      type == "hat" ~ "Hat Occlusion",
      type == "vis" ~ "Other Occlusion"
    ),
    type = factor(type, levels = c("Other Occlusion", "Hat Occlusion", "No Coded AU", "Any Other AU"))
  )

fig1df |>
  mutate(.by = wave, pct = sec / sum(sec)) |>
  gt(groupname_col = "wave") |>
  fmt_percent(columns = pct) |>
  tab_options(row_group.as_column = TRUE)

	type	sec	pct
6	Other Occlusion	552.901	10.14%
	Hat Occlusion	0.000	0.00%
	No Coded AU	2138.597	39.23%
	Any Other AU	2759.610	50.62%
9	Other Occlusion	872.787	15.29%
	Hat Occlusion	177.454	3.11%
	No Coded AU	1594.357	27.93%
	Any Other AU	3064.487	53.68%
12	Other Occlusion	275.052	5.14%
	Hat Occlusion	3319.348	62.06%
	No Coded AU	694.450	12.98%
	Any Other AU	1059.890	19.82%
15	Other Occlusion	293.760	8.43%
	Hat Occlusion	2430.587	69.72%
	No Coded AU	202.515	5.81%
	Any Other AU	559.495	16.05%

fig1df |>
  ggplot(aes(x = wave, y = sec, fill = type)) +
  geom_col(position = "fill", color = "black") +
  scale_fill_manual(values = c("grey30", "grey60", "white", "steelblue")) +
  scale_x_discrete(
    labels = c("6-8 months", "9-11 months", "12-14 months", "15-17 months")
  ) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Age Window",
    y = "Proportion of full-face time",
    fill = "Type"
  ) +
  theme_bw(base_size = 10)

ggsave(filename = "fig1.png", width = 6.5, height = 4.0, units = "in", dpi = 300)