---
title: >
  "Objective Assessment of Wear Time During the Initial Phase of Orthodontic Aligner Therapy Using Microsensors - 
  a Randomized Controlled Trial"
subtitle: "Data Supplement - Manuscript (ref. AJODO-D-25-00406R2) accepted 2025/10/03"
author: 
  - name: Uwe Baumert
    affiliation:
      - id: LMU
      - name: Department of Orthodontics and Dentofacial Orthopedics, Dental School, LMU Medical Center, LMU Munich, Germany
  - name: Hisham Sabbagh
    affiliation:
      - ref: LMU
date: 2025-04-21
date-modified: 2025-10-08
editor: source
bibliography: grateful-refs.bib
format: 
  html:
    theme: cosmo
    fontsize: 1.1em
    linestretch: 1.7
    toc: true
    toc-depth: 4
    toc-title: Contents
    toc-location: left
    number-sections: true
    number-depth: 4
    fig-align: center
    fig-format: png
    fig-cap-location: bottom
    tbl-cap-location: top
    page-layout: full
    embed-resources: true
    link-external-newwindow: true
    html-math-method: katex
    html-table-processing: none
    code-fold: true
    code-tools: true
    code-summary: "Show the code"
    code-links:
      - text: Data file
        icon: download
        href: Theramon_data_long2.sav
    code-block-bg: true
    code-block-border-left: "#31BAE9"
    css: style.css
    crossref:
      chapters: true
execute: 
  echo: true
  cache: true
---

```{r}
#| label: r-setup
#| echo: true
#| error: false
#| warning: false
#| message: false
#| include: false

# knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(broom)
library(ggpubr)
library(sjlabelled)
library(labelled)
library(ggpmisc)
library(rstatix)
library(gghighlight)
library(patchwork)
library(knitr)
library(kableExtra)
library(scales)
library(gtsummary)
library(flextable)
library(qreport)
library(grid)
library(ggsci)
library(ggmosaic)
library(grateful)
library(lme4)
library(lmerTest)
library(broom.mixed)
library(texreg)
library(easystats)
library(report)
library(effects)
library(sjPlot)
library(ggrepel)
library(emmeans)

set.seed(2)
options(ggplot2.discrete.colour = "npg") # ggpubr - palette definition

## formatter for flextable() output; from URL https://github.com/davidgohel/flextable/issues/118 (2025-07-07)
p_val_format <- function(x){
  z <- scales::pvalue_format()(x)
  z[!is.finite(x)] <- ""
  z}

options("digits" = 2)

setwd("~/Documents/01_Manuskripte/Sabbagh_Theramon_AJODO/R")

# load("./Theramon_data-supplement_accepted.RData")
```


# Data preparation and code book


```{r}
#| label: data-import
#| echo: true
#| error: false
#| warning: false
#| message: false
#| include: true

theramon_long <- haven::read_sav("Theramon_data_long2.sav")

theramon_long <- theramon_long |>
  mutate(Id = factor(Id), 
         age = as.numeric(sprintf(age_yr, fmt = "%#.1f")), # 'age' rounded to 1 significant digit
         study_arm = factor(study_arm, 
                            levels = c("0", "1"), 
                            labels = c("Group 1 (unaware)", "Group 2 (aware)")),
         sex = factor(sex, 
                      levels = c("0", "1"), 
                      labels = c("Male", "Female")),
         visit_n = as.integer(visit),
         visit = factor(visit, levels = c("1", "2", "3", "4", "5", "6"),
                              labels = c("T0-T1", "T1-T2", "T2-T3", "T3-T4", "T4-T5", "T5-T6")),
         row = row_number()) |>
  rename(wear_time.visit = wear_appointmentx_h,
         mean_dApp_hours = mean_dApp_days,
         delta_visit.days = diff_Appointment_days,
         delta_visit.hours = dApp_h) |>
  select(Id, study_arm, age, sex, visit, visit_n, wear_time.visit, delta_visit.days, delta_visit.hours, mean_wear_h, mean_diff_days, mean_dApp_hours, row)
# [1] "Id"                "study_arm"         "age"               "sex"               "visit"             "visit_n"          
# [7] "wear_time.visit"   "delta_visit.days"  "delta_visit.hours" "mean_wear_h"       "mean_diff_days"    "mean_dApp_hours"  

theramon_long <- theramon_long |>
  group_by(Id) |>
  mutate(cum_days = cumsum(delta_visit.days)) |>
  ungroup() |>
  var_labels(
    Id = "Id",
    age = "Age (years)",
    study_arm = "Study arm",
    visit_n = "Period number",
    visit = "Period",
    sex = "Sex",
    wear_time.visit = "WT at visit (h)",
    delta_visit.days = "Interval between visits (days)",
    delta_visit.hours = "Change in WT between visits (h)",
    cum_days = "Cumulative treatment time (days)",
    mean_wear_h = "Mean WT (h)",
    mean_dApp_hours = "Mean WT change (h)",
    row = "Observation number")

#names(theramon_long)
#  [1] "Id"                "study_arm"         "age"               "sex"               "visit"             "visit_n"          
#  [7] "wear_time.visit"   "delta_visit.days"  "delta_visit.hours" "mean_wear_h"       "mean_diff_days"    "mean_dApp_hours"  
# [13] "row"               "cum_days"    
```

Originally, the data was collected in a Microsoft Excel spreadsheet in wide format. 
In accordance with the study protocol, the following data was entered anonymized for each participant after completion of the 6th wear period:

- `Id` (factor; 1...40), 
- `sex` (factor; male = 0, female = 1), 
- `age` (continuous; years),
- `study_arm` (factor; Group 1 (unaware) = 0, Group 2 (aware) = 1), 
- and data exported from the *Theramon* program:
  - for each wear period the duration of the specific period (days) (`delta_visit.days`) and 
  - the wear time (WT) in hours (`wear_time.visit`).

This data was used to calculate other variables for each patient:

- change in mean wear time between visits (h) (`delta_visit.hours`) for each study subject
- and mean wear time (h) (`mean_wear_h`) for each study subject. 

Do aid further analysis, the `row` number corresponding to an individual observation was added to the `theramon_long` data frame ("Observation number").

For analysis, this table was reshaped to "long" format in Excel and then imported in SPSS 29. 
Each variable was labelled accordingly and if necessary, factor levels were defined.
The statistical analysis and figure preparation was done using this data file (`Theramon_data_long.sav`). 

The final statistical analysis was done using R [@base] using RStudio version 2025.5.1.513. All analysis and output was documented in this Quarto file.

In R, the variable `cum_days` representing cumulative treatment time (days) was additionally calculated using the R base function `cumsum`. 

Data is available in the data frame `theramon_long`, which is in "long" format. Variable definitions are summarized in @tbl-codebook1.


```{r}
#| label: tbl-codebook1
#| tbl-cap: "Variable definitions (aka 'codebook'). WT = wear time"
#| echo: true
#| error: false
#| warning: false
#| message: false

labelled::generate_dictionary(theramon_long) |> 
  labelled::convert_list_columns_to_character() |> 
  select(pos, variable, label, col_type, levels) |> 
  rename(Pos = pos, 
         Variable = variable,
         Label = label, 
         "Col. type" = col_type, 
         Levels = levels)  |>
  flextable() |>
  bold(bold = TRUE, part = "header") |>
  autofit()
```


# Baseline data


## Patient's `age`


In the following table (@tbl-desc-age-study_arm), a summary of `age` stratified by study arm is given.


```{r}
#| label: tbl-desc-age-study_arm
#| tbl-cap: "Descriptive statistics for `age` by `study_arm`."
#| echo: true
#| error: false
#| warning: false
#| message: false

set.seed(2)

theramon_long |>
  filter(visit == "T0-T1") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  tbl_summary(include = c(age),
              by = study_arm,
              statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} [{p25}, {p75}]", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1) |>
  add_overall() |>
  add_n() |>
  as_flex_table() |>
  set_table_properties(width = 0.7)
```


## Variable: `sex`


Cross-tabulation was applied to investigate the distribution of `sex` by `study_arm` (@tbl-stat-baseline-sex).


```{r}
#| label: tbl-stat-baseline-sex
#| tbl-cap: "Fisher's exact test: `sex` by `study_arm`."
#| echo: true
#| error: false
#| warning: false
#| message: false

set.seed(2)

theramon_long |>
  filter(visit == "T0-T1") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(sex, study_arm) |>
  tbl_summary(by = study_arm) |>
  add_overall() |>
  as_flex_table() |>
  set_table_properties(layout = "autofit") 
```


A mosaic plot shows the properties (@fig-mosaic-baseline-sex).


```{r}
#| label: fig-mosaic-baseline-sex
#| fig-cap: "Mosaicplot: `study_arm` by `sex`."
#| echo: true
#| error: false
#| warning: false
#| message: false
#| fig-width: 2.5
#| fig-height: 2.5
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg")

mosaic.data <- theramon_long |>
  filter(visit == "T0-T1") |>
  select(sex, study_arm) |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  as_tibble()

mosaic.ggplot <- ggplot(mosaic.data) +
  geom_mosaic(aes(x = product(sex, study_arm), fill = study_arm)) +
  scale_fill_npg() +
  labs(x = "Study arm", y = "Sex") +
  theme_bw(base_size = 12) +
  theme(axis.title.x = element_text(size = 12, face = "bold"),
        axis.text.x = element_text(size = 10),
        axis.title.y = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 10)) |>
  theme(legend.position = "below")

mosaic.ggplot

# since ggplot >= 4.0.0 breaks the current available 'ggmosaic' 0.3.3
# the mosaic plot can also be generated using 'vcd' 
#vcd::mosaic(~ sex + study_arm, data = mosaic.data, highlighting = "study_arm", highlighting_fill = pal_npg(palette = c("nrc"))(2))

```


## Baseline summary


Baseline summary statistics is shown in @tbl-stat-baseline-summary.


```{r}
#| label: tbl-stat-baseline-summary
#| tbl-cap: "Summary statistics of baseline data."
#| echo: true
#| error: false
#| warning: false
#| message: false

set.seed(2)

theramon_long |>
  filter(visit == "T0-T1") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  tbl_summary(include = c(age, sex),
              by = study_arm,
              statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1) |>
  add_overall() |>
  add_n() |>
  as_flex_table() |>
  set_table_properties(layout = "autofit")
```


# Outcomes


## Time interval (days) between two consecutive appointments


Variables:

- `delta_visit.days`: duration of the individual wear periods T0-T1; 2: T1-T2; ...; T5-T6;
- `mean_diff_days`: covers duration of all individual wear periods (T0-T1 up to and including T5-T6, i.e. T0-T6) and is the mean of these for each individual case.

Patients were scheduled for six clinical appointments following a biweekly protocol.
However, due to factors such as professional obligations, illness, and other circumstances, strict adherence to the prescribed schedule was not always feasible.
In @tbl-desc-delta_visit_days the intervals between two consecutive appointments (i.e., wear periods) and the overall mean time of all wear periods covering 
periods T0-T1, ..., T5-T6 were summarized.


```{r}
#| label: tbl-desc-delta_visit_days
#| tbl-cap: "Number of days between appointments and overall time between appointments"
#| echo: true
#| error: false
#| warning: false
#| message: false

set.seed(2)

tbl_mean_diff_days <- theramon_long |>
  filter(visit == "T5-T6") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  tbl_summary(include = mean_diff_days,
                by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(mean_diff_days = "Period T0-T6")) |>
  add_overall()

tbl_delta_visit.1 <- theramon_long |>
  filter(visit == "T0-T1") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T0-T1")) |>
  add_overall()

tbl_delta_visit.2 <- theramon_long |>
  filter(visit == "T1-T2") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T1-T2")) |>
  add_overall()

tbl_delta_visit.3 <- theramon_long |>
  filter(visit == "T2-T3") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T2-T3")) |>
  add_overall()

tbl_delta_visit.4 <- theramon_long |>
  filter(visit == "T3-T4") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T3-T4")) |>
  add_overall()

tbl_delta_visit.5 <- theramon_long |>
  filter(visit == "T4-T5") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T4-T5")) |>
  add_overall()

tbl_delta_visit.6 <- theramon_long |>
  filter(visit == "T5-T6") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, delta_visit.days) |>
  tbl_summary(by = study_arm,
              type = list(all_continuous() ~ "continuous2",
                          all_categorical() ~ "continuous2"),
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(delta_visit.days = "Period T5-T6")) |>
  add_overall()

tbl_stack(list(tbl_delta_visit.1,
          tbl_delta_visit.2,
          tbl_delta_visit.3,
          tbl_delta_visit.4,
          tbl_delta_visit.5,
          tbl_delta_visit.6,
          tbl_mean_diff_days)) |>
  as_flex_table() |>
  set_table_properties(layout = "autofit")
```


The abacus plot (@fig-abacus-plot) shows the distribution of visits among the patients, illustrating both the clustering of appointment waves and 
the individuals' variability in time intervals between consecutive appointments.


```{r}
#| label: fig-abacus-plot
#| fig-cap: "Abacus plot: visit waves for each study subject. Grey shaded areas label the 14-days visit intervals as proposed."
#| echo: true
#| error: false
#| warning: false
#| message: false
#| fig-height: 10
#| fig-width: 10
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg")

data_cumdays <- theramon_long |>
  select(Id, study_arm, visit, cum_days) |>
  mutate(Id = as.numeric(Id),
         visit = as.numeric(visit))|>
  # add_row(Id = 1:40, visit = 0, cum_days = 0) |>
  add_row(Id = c(1:2, 4:7, 9:11, 13, 15, 18, 21, 25, 26, 29, 31, 34, 36, 37), visit = 0, cum_days = 0, study_arm = "Group 1 (unaware)") |>
  add_row(Id = c(3, 8, 12, 14, 16, 17, 19, 20, 22:24, 27, 28, 30, 32, 33, 35, 38:40), visit = 0, cum_days = 0, study_arm = "Group 2 (aware)") |>
  mutate(visit = factor(visit, levels = c("0", "1", "2", "3", "4", "5", "6"))) |>
  arrange(Id, visit)

fig.abacus_plot_cumdays_with0_high <- data_cumdays |>
  ggplot(aes(x = fct_rev(as.factor(Id)), 
             y = cum_days, 
             group = Id, colour = study_arm)) +
  geom_point(size = 1.2, aes(colour = study_arm)) +
  scale_y_continuous(name = "Visit waves (days)", limits = c(0, 120), breaks = seq(0, 120, by = 10), expand = c(0,1)) + # 10-tägig
  labs(x = "Case Id") +
  guides(color = guide_legend(title = "Study arm:")) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 13, ymax = 15, alpha = .2) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 27, ymax = 29, alpha = .2) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 41, ymax = 43, alpha = .2) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 55, ymax = 57, alpha = .2) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 69, ymax = 71, alpha = .2) +
  annotate("rect", xmin = 0, xmax = Inf, ymin = 83, ymax = 85, alpha = .2) +
  coord_flip() +
  theme_bw(base_size = 14) +
  theme(axis.title.x = element_text(size = 14, face = "bold"),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14, face = "bold"),
        axis.text.y = element_text(size = 12)) +
  theme(panel.grid.minor.x = element_blank()) +
  theme(legend.position = "bottom")

fig.abacus_plot_cumdays_with0_high
```


```{r}
#| label: fig_s1_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_S1 panel_visit-waves.pdf", plot = fig.abacus_plot_cumdays_with0_high, device = cairo_pdf, units = "in", width = 10, height = 8, dpi = 300)
ggsave("Fig_S1 panel_visit-waves.png", plot = fig.abacus_plot_cumdays_with0_high, device = png, units = "in", width = 10, height = 8, dpi = 300)
```


## Wear period specific wear time

Variables:

* `wear_time.visit`: wear time for each case/wear period
* `mean_wear_h`: mean wear time over all wear periods for each case

Wear times per wear period and the overall mean wear time (periods T0-T1 ... T5-T6) were summarized for both study arms
and the complete study group in @tbl-desc-mean-WT-between-visits.


```{r}
#| label: tbl-desc-mean-WT-between-visits
#| tbl-cap: "Descriptive statistics of `wear_time.visit` by `study_arm` and `visit` and `mean_wear_h` by `study_arm`."
#| echo: true
#| error: false
#| warning: false
#| message: false

tbl_mean_wear_h <- theramon_long |>
  filter(visit == "T5-T6") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  tbl_summary(include = mean_wear_h,
              by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(mean_wear_h = "Period T0-T6")) |>
  add_overall()

tbl_wear_time_visit.1 <- theramon_long |>
  filter(visit == "T0-T1") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T0-T1")) |>
  add_overall()

tbl_wear_time_visit.2 <- theramon_long |>
  filter(visit == "T1-T2") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T1-T2")) |>
  add_overall()

tbl_wear_time_visit.3 <- theramon_long |>
  filter(visit == "T2-T3") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T2-T3")) |>
  add_overall()

tbl_wear_time_visit.4 <- theramon_long |>
  filter(visit == "T3-T4") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T3-T4")) |>
  add_overall()

tbl_wear_time_visit.5 <- theramon_long |>
  filter(visit == "T4-T5") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T4-T5")) |>
  add_overall()

tbl_wear_time_visit.6 <- theramon_long |>
  filter(visit == "T5-T6") |>
  mutate(study_arm = recode(study_arm, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  select(study_arm, wear_time.visit) |>
  tbl_summary(by = study_arm,
                            statistic = list(all_continuous() ~ "{mean} ({sd}) \n {median} ({p25}, {p75})", all_categorical() ~ "{n} ({p}%)"),
              digits = all_continuous() ~1,
              label = list(wear_time.visit = "Period T5-T6")) |>
  add_overall()

tbl_stack(list(tbl_wear_time_visit.1,
          tbl_wear_time_visit.2,
          tbl_wear_time_visit.3,
          tbl_wear_time_visit.4,
          tbl_wear_time_visit.5,
          tbl_wear_time_visit.6,
          tbl_mean_wear_h)) |>
  as_flex_table() |>
  set_table_properties(layout = "autofit")
```


Trajectories of each participant overlaid with the mean WT (and its 95% CI) at each visit for both study arms are shown in @fig-indiv-mean-WT-studyarms.


```{r}
#| label: fig-indiv-mean-WT-studyarms
#| fig-cap: "Individual mean WT for the six wear periods and both study arms"
#| echo: true
#| error: false
#| warning: false
#| message: false
#| fig-width: 12
#| fig-height: 5
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg")

fig.study_indiv <- theramon_long |>
  select(Id, study_arm, visit, wear_time.visit) |>
  mutate(visit = as.numeric(visit)) |>
  arrange(Id, visit) |>
  ggplot(aes(x = visit, y = wear_time.visit, group = Id)) +
  geom_line(colour = "lightgrey") +
  geom_point(colour = "darkgrey") +
  guides(shape = "none", color = "none") +
  labs(x = "Wear period", y = "Mean wear time (h) with 95% CI") +
  scale_x_continuous(limit = c(.5, 6.5), breaks = seq(1, 6, by = 1), expand = c(0,0), labels = c("T0-T1", "T1-T2", "T2-T3", "T3-T4", "T4-T5", "T5-T6")) +
  scale_y_continuous(limit = c(0, 25), breaks = seq(0, 25, by = 10), expand = c(0,0)) +
  facet_grid(cols = vars(study_arm), labeller = labeller(study_arm = c("0" = "Temperature", "1" = "Wear time"))) +
  theme_bw(base_size = 14) +
  theme(axis.title.x = element_text(size = 14, face = "bold"),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14, face = "bold"),
        axis.text.y = element_text(size = 12))+
  theme(legend.position = "bottom") +
  theme(panel.grid.major.y = element_line(color = "darkgrey", linetype = "dashed"), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank()) +
  theme(strip.text.x = element_text(size = 12, face = "bold"))

fig.study_indiv2 <- fig.study_indiv + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0.2, color = "red", group = 1, linewidth = .8) + 
  stat_summary(fun = "mean", geom = "point", shape = 16, size = 2.5, color = "red", group = 1)

fig.study_indiv2
```


```{r}
#| label: fig_4_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_4 panel_individuelle_Tragezeit_Id-wise_CI.pdf", plot = fig.study_indiv2, device = cairo_pdf, units = "in", width = 10, height = 4, dpi = 300)
ggsave("Fig_4 panel_individuelle_Tragezeit_Id-wise_CI.png", plot = fig.study_indiv2, device = png, units = "in", width = 10, height = 4, dpi = 300)
```


The individual patters during the six periods of this study are shown in the next figure (@fig-indiv-mean-WT-across-visits). 
Ordinary least square regression was applied to visualize the individual differences as indicated by the different, individual slopes.


```{r}
#| label: fig-indiv-mean-WT-across-visits
#| fig-cap: "Wear time (WT) for each individual case across the six follow-up appointments. Slopes were calculated for each case separately."
#| echo: true
#| error: false
#| warning: false
#| message: false
#| fig-width: 12
#| fig-height: 14
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg") 

theramon_long.slopes <- theramon_long |>
  group_by(Id) |>
  do(tidy(lm(wear_time.visit ~ as.numeric(visit), data = .))) |>
  select(Id, term, estimate) |>
  pivot_wider(names_from = term, values_from = estimate) |>
  rename(intercept = '(Intercept)', 
         slope = 'as.numeric(visit)') |>
  ungroup() 

theramon_long.w_slopes <- theramon_long |>
  group_by(Id) |>
  do(tidy(lm(wear_time.visit ~ as.numeric(visit), data = .))) |>
  select(Id, term, estimate) |>
  pivot_wider(names_from = term, values_from = estimate) |>
  rename(intercept = '(Intercept)', 
         slope = 'as.numeric(visit)') |>
  ungroup()  |>
  merge(theramon_long, by = 'Id')

case_names <- as_labeller(function(x) paste0('Case ', x))

## https://github.com/aphalo/ggpmisc/issues/12
## https://stackoverflow.com/questions/7549694/add-regression-line-equation-and-r2-on-graph
## https://www.reddit.com/r/rstats/comments/bwp5sg/how_do_i_get_r2_values_and_line_equation_for_geom/

theramon_long.slopes_round <- theramon_long.slopes |>
  mutate_if(is.numeric, round, 3)

fig.wear_time.Id_lm <- theramon_long |>
  arrange(Id, visit) |>
  ggplot(aes(x = visit, y = wear_time.visit, group = 1, color = study_arm)) +
  geom_line() +
  geom_point(shape = 21, size = 2, fill = "white") +
  stat_poly_line(formula = y ~ x, colour = "black", linewidth = 0.5, se = FALSE) +
  geom_text(data = theramon_long.slopes_round, aes(x = 2, y = 25, label = paste0("slope: ", slope)), colour = "black", size = 3) +
  guides(shape = "none") + 
  labs(x = "Wear period", y = "Wear time (h)") +
  scale_y_continuous(limit = c(0, 30), breaks = seq(0, 30, by = 10), expand = c(0,0)) +
  scale_colour_discrete(name = "Study arm", labels = c("Group 1 (unaware)", "Group 2 (aware)")) +
  facet_wrap(facet = vars(Id), ncol = 5, labeller = case_names) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  theme(panel.grid.major.y = element_line(color = "lightgrey", linetype = "dashed"), 
        panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank()) +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 30, vjust = 0.8, hjust = 0.8))

fig.wear_time.Id_lm
```


```{r}
#| label: fig_3_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_3 individuelle_Tragezeit_Id-wise_lm.pdf", plot = fig.wear_time.Id_lm, device = cairo_pdf, units = "in", width = 12, height = 14, dpi = 300)
ggsave("Fig_3 individuelle_Tragezeit_Id-wise_lm.png", plot = fig.wear_time.Id_lm, device = png, units = "in", width = 12, height = 14, dpi = 300)
```


To get an idea on the different patters, the slopes calcuated for @fig-indiv-mean-WT-across-visits above, were used to classify the study subjects into those
with increasing slope over the six periods (*slope* >= 0.1), with constant slope over the six periods (-0.1 <= *slope* <= 0.1),
and with decreasing slope (*slope* < -0.1). The findings were summarized in @fig-indiv-trajectories.


```{r}
#| label: fig-indiv-trajectories
#| fig-cap: 
#|   - "Individual participant trajectories. (a) increasing slope (*slope* >= 0.1), \n (b) constant slope (-0.1 ≤ *slope* ≤ 0.1 ), and (c) decreasing slope (*slope* < -0.1)."
#| echo: true
#| error: false
#| warning: false
#| message: false
#| fig-height: 12
#| fig-width: 8
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg")

fig.theramon_slope <- theramon_long.w_slopes |>
  ggplot(aes(x = as.numeric(visit), y = wear_time.visit, group = Id)) +
  geom_line(aes(colour = Id), linewidth = 1.5) +
  labs(x = "Wear period", y = "Wear time (h)") +
  scale_x_continuous(limit = c(.5, 6.5), breaks = seq(1, 6, by = 1), expand = c(0,0), labels = c("T0-T1", "T1-T2", "T2-T3", "T3-T4", "T4-T5", "T5-T6")) +
  scale_y_continuous(limit = c(0, 25), breaks = seq(0, 25, by = 10), expand = c(0,0)) +
  facet_grid(cols = vars(study_arm)) +
  theme_bw(base_size = 14) +
  theme(axis.title.x = element_text(size = 14, face = "bold"),
        axis.title.y = element_text(size = 14, face = "bold"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12)) +
  theme(panel.grid.major.y = element_line(color = "darkgrey", linetype = "dashed"), 
        panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank()) +
  theme(strip.text.x = element_text(size = 12, face = "bold")) +
  theme(legend.position = "none")

# fig.theramon_slope

## Increasing slop from visit 1 to 6
# theramon_long.w_slopes |>
#   filter(slope > 0.1) |>
#   filter(as.numeric(visit) == 1 | as.numeric(visit) == 6) |>
#   group_by(Id, visit) |>
#   select(Id, visit, slope, sex, wear_time.visit, study_arm) |>
#   pivot_wider(names_from = visit, values_from = c(wear_time.visit)) |>
#   kable(digits = 3)

fig.theramon_slope_inc <- fig.theramon_slope + 
  gghighlight(slope  > 0.1, calculate_per_facet = TRUE, keep_scales = TRUE, label_key = Id, 
              use_direct_label = FALSE, unhighlighted_params = list(linewidth = .8)) 

## more or less stable slope from visit 1 to 6
# theramon_long.w_slopes |>
#   filter(slope <= 0.1 & slope >= -0.1) |>
#   filter(as.numeric(visit) == 1 | as.numeric(visit) == 6) |>
#   group_by(Id, visit) |>
#   select(Id, visit, slope, sex, wear_time.visit, study_arm) |>
#   pivot_wider(names_from = visit, values_from = wear_time.visit) |>
#   kable(digits = 3)

fig.theramon_slope_stable <- fig.theramon_slope +
  gghighlight(max(slope) <=0.1, min(slope) >= -0.1, calculate_per_facet = TRUE, keep_scales = TRUE, 
              label_key = Id, use_direct_label = FALSE, unhighlighted_params = list(linewidth = .8))

## Decreasing slope from visit 1 to 6
# theramon_long.w_slopes |>
#   filter(slope < -0.1) |>
#   filter(as.numeric(visit) == 1 | as.numeric(visit) == 6) |>
#   group_by(Id, visit) |>
#   arrange(Id) |>
#   select(Id, visit, slope, sex, wear_time.visit, study_arm) |>
#   pivot_wider(names_from = visit, values_from = wear_time.visit) |>
#   kable(digits = 3)

fig.theramon_slope_dec <- fig.theramon_slope + 
  gghighlight(min(slope) < -0.1,  calculate_per_facet = TRUE, keep_scales = TRUE, label_key = Id, 
              use_direct_label = FALSE, unhighlighted_params = list(linewidth = .8)) 

grid3 <- grid::grid.draw(fig.theramon_slope_inc + 
                           fig.theramon_slope_stable + 
                           fig.theramon_slope_dec + 
                           plot_layout(ncol = 1) + 
                           plot_annotation(tag_levels = 'a'))
```


```{r}
#| label: fig_S2_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_S2 panel_indiv_Tragezeit_slopes.pdf", plot = grid3, device = cairo_pdf, units = "in", width = 10, height = 15, dpi = 300)
ggsave("Fig_S2 panel_indiv_Tragezeit_slopes.png", plot = grid3, device = png, units = "in", width = 10, height = 15, dpi = 300)
```


# Linear mixed model


Subject-specific trajectories (@fig-indiv-mean-WT-across-visits and @fig-indiv-trajectories) revealed
extensive variation in `wear_time.visit` between subjects. Within each subject `wear_time.visit` was also fluctuating over time, 
and either slightly increasing, nearly constant, or decreasing (@fig-indiv-trajectories).

Due to the nested structure of the data, a linear mixed model (LMM) was established to test the effect of *time* (variable `visit`; ordered), 
study arm* (variable `study_arm`; nominal; Group 1 (unaware) = 0, Group 2 (aware) = 1), and the fixed predictors (`sex`, `age`) on wear period specific
wear time (h) (variable `wear_time.visit`), and potential interactions between them after controlling for subject (`Id`) variation. 


## Model selection


Models of different complexity were evaluated in this order:

1. `fit_lmer_1ml`:   Intercept model: first and second level predictors (`visit` and `study_arm`) and the random intercept (`Id`).
2. `fit_lmer_2aml`:  Potential interaction between `visit` and `study_arm`.
3. `fit_lmer_3aml`:  Addition of `sex`, and `age`.
4. `fit_lmer_4ml`:   Potential interaction between `study_arm` and `sex`.
5. `fit_lmer_4bml`:  Potential interaction between `study_arm` and `age`.

In all of the following steps, REML was deactivated (*REML = FALSE*) to ease model comparisons as suggested by the `lme4` R package authors [@lme4; @lmerTest].


### Intercept model

The first model (`fit_lmer_1ml`) consisted of `wear_time.visit` as response variable, `visit` (time; first level) and `study_arm` (groups; second level), 
and the subject-contributed random effects (`Id`). The latter were described by:

1. the random intercept, which is an estimate for the variation of the intercept and as such wear time at visit (`wear_time.visit`) between subjects (`Id`).
2. the level 1-residues, describing the variation of `wear_time.visit`  within a study subject.


```{r}
#| label: mm_visit_studyarm
#| include: true
#| error: false
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)

fit_lmer_1ml <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + (1 | Id), data = theramon_long, REML = FALSE)
summary(fit_lmer_1ml)
```


### Interaction between `visit` and `study_arm`


Next, a potential interaction between time and study arm (`visit:study_arm`) was analyzed and compared this model to the previous one without interaction, i.e., `fit_lmer_1ml`.


```{r}
#| label: mm_visit-studyarm-interaction
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)

fit_lmer_2ml <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + visit:study_arm + (1 | Id), data = theramon_long, REML = FALSE)
summary(fit_lmer_2ml) |> print(width = 100)
```


To test, if the addition of `study_arm:visit` interaction was contributing further information to the `fit_lmer_1ml` model a likelihood ratio test was applied:


```{r}
#| label: mm_visit-studyarm-interaction_LR
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)

anova(fit_lmer_1ml, fit_lmer_2ml)
```


According to this results, the interaction between `visit` and `study_arm` was not significant. 
Therefore, this interaction was not further considered and we continued with the less complex model `fit_lmer_1ml`.


### Addition of the level-invariant predictors `sex` and `age`


Next, we added the level-invariant predictors `sex` and `age`to the model from the previous step, yielding model 

`fit_lmer_3aml`: `wear_time.visit ~ visit + study_arm + sex + age + (1|Id)`.

To avoid, that the reference of `age` is set to 0 years, the minimum age of this patient cohort (i.e., `{r} round(min(theramon_long$age), 1)` years) 
was subtracted from `age` to yield `age-min(age)` as reference.


The model's summary is shown next:


```{r}
#| label: mm_visit_studyarm_sex_age
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)
options(digits = 3)

fit_lmer_3aml <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + I(age-17.2) + (1 | Id), data = theramon_long, REML = FALSE)
summary(fit_lmer_3aml)
```


Next, we tested if either `age` (i.e., `I(age-17.2)`) or `sex` were significantly contributing to the model using a likelihood ratio test.


```{r}
#| label: mm_visit_studyarm_sex_age_LR
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)

anova(fit_lmer_3aml)
```

According to these results, `age` was not significant and thus removed from the model. 

We continued with model `fit_lmer_3ml` containing `visit`, `study_arm`, and `sex` as predictors, and `Id` as random intercept:


```{r}
#| label: mm_visit_studyarm_sex
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)
options(digits = 3)

fit_lmer_3ml <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + (1 | Id), data = theramon_long, REML = FALSE)
summary(fit_lmer_3ml)
```


### Interaction between `study_arm` and `sex`


Then, a potential interaction between `study_arm` and `sex` (i.e., `study_arm:sex`) was tested, yielding `fit_lmer_4ml`.


```{r}
#| label: mm_visit_age_studyarm-sex-interaction
#| include: true
#| error: false
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)

fit_lmer_4ml <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long, REML = FALSE)
summary(fit_lmer_4ml)
```


### Model comparison


Model comparison between the 'null' model and the models accepted so far (@tbl-mixed-model-comparison).


```{r}
#| label: tbl-mixed-model-comparison
#| tbl-cap: "Model comparision between the models `fit_lmer_1ml`, `fit_lmer_3ml`, \n and `fit_lmer_4aml`."
#| error: false
#| include: true
#| warning: false
#| message: false
#| results: asis

set.seed(2)
options(ggplot2.discrete.colour = "npg")

sjPlot::tab_model(fit_lmer_1ml, fit_lmer_3ml, fit_lmer_4ml, show.reflvl = TRUE, p.val = "satterthwaite", df.method = "satterthwaite",
                  show.aic = TRUE, show.loglik = TRUE,
                  dv.labels = c("fit_lmer_1ml", "fit_lmer_3ml", "fit_lmer_4ml"))
```


Additionally, the performance of these models using `compare_performance()` (the "performance" R package) was compared (@tbl-mixed-model-performance).


```{r}
#| label: tbl-mixed-model-performance
#| tbl-cap: "Ranked mixed model performance."
#| include: true
#| error: false
#| warning: false
#| message: false

set.seed(2)
options(ggplot2.discrete.colour = "npg")

performance::compare_performance(fit_lmer_1ml, fit_lmer_3ml, fit_lmer_4ml, rank = TRUE) |> 
  select("Name", "R2_conditional", "R2_marginal", "ICC", "RMSE", "Sigma", "AIC_wt", "AICc_wt", "BIC_wt", "Performance_Score") |>
  rename("R2 cond." = R2_conditional,
         "R2 marg." = R2_marginal,
         "Performance Score" = Performance_Score) |>
  as_flextable() |> 
  colformat_double(digits = 3) |>
  bold(bold = TRUE, part = "header") |> autofit()
```


Taken together, we accepted `fit_lmer_4ml`.


## Accepted model

### General description of the accepted model

The formula of the accepted model was: 

`wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id)`. 

This time, REML was enabled. Degrees of freedom and p-values were calculated with the Kenward-Roger approximation.


```{r}
#| label: mm_final-model-summary
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)

final <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long, REML = TRUE)

# Degrees of freedom and t-statistics were computed using Kenward-Roger's method.
summary(final, ddf = "Kenward-Roger")
```



The model's performance (`performance::performance(final)`) is summarized in @tbl-mm_final-model_performance.


```{r}
#| label: tbl-mm_final-model_performance
#| tbl-cap: "Performance of the accepted model `final`."
#| include: true
#| error: false
#| warning: false
#| message: false

set.seed(2)

performance::performance(final) |>
  rename("R2_cond" = R2_conditional,
         "R2_marg" = R2_marginal) |>
  flextable() |> 
  colformat_double(j = c("AIC", "AICc", "BIC"), digits = 0) |>
  colformat_double(j = c("R2_cond", "R2_marg", "ICC", "RMSE", "Sigma"), digits = 3) |>
  bold(bold = TRUE, part = "header") |> 
  set_table_properties(layout = "autofit")
```


### LMM model assumptions

The `final` model was then analyzed for deviations from the assumptions of linear mixed models using `performance::check_model()`
from the "performance" R package (@fig-model-assumptions).


```{r}
#| label: fig-model-assumptions
#| fig-cap: "Accepted model overview. Visual check of various model assumptions using `performance::check_model()`."
#| error: false
#| include: true
#| warning: false
#| message: false
#| fig-width: 10
#| fig-height: 12
#| fig-align: center

final_performance.fig <- performance::check_model(final)
final_performance.fig

# cairo_pdf("Fig_S4 final_performance.pdf", family = "sans", height = 15, width = 10)
# final_performance.fig
# dev.off()
```


#### Heteroscedasticity


Checking the assumption of heteroscedasticity using the Breusch-Pagan test was applied as implemented in the `check_heteroscedasticity()` function
implemented in the "performance" R package [@performance] and depicted (@fig-mm-assumptions-heteroscedasticity)


```{r}
#| label: fig-mm-assumptions-heteroscedasticity
#| fig-cap: "Check for homogeneity of variance"
#| error: false
#| include: true
#| warning: false
#| message: false
#| fig-width: 4
#| fig-height: 2.5
#| fig-align: center

plot(performance::check_heteroscedasticity(final))
```


```{r}
#| label: mm-assumptions-heteroscedasticity
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

performance::check_heteroscedasticity(final)
```


#### Outliers {#sec-lmm-outliers}


Check for potential outliers using `check_outliers()` from the "performance" R package using Cook's distance.


```{r}
#| label: mm-assumptions-outliers-cook
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

performance::check_outliers(final)
```


Observation "145" (Id = 3, visit = T3-T4) was identified as a potential outlier according to Cook's distance.
The related data is summarized in the following table (@tbl-assumptions-outliers_observations).


```{r}
#| label: tbl-assumptions-outliers_observations
#| tbl-cap: "Observations identified as potential outliers by observation #145."
#| include: true
#| error: false
#| warning: false
#| message: false

theramon_long |>
  mutate(Row = row_number()) |>
  filter(Row == 145) |>
  select(Row, Id, study_arm, sex, visit, wear_time.visit) |>
  flextable() |> 
  colformat_num(digits = 3) |>
  colformat_double(digits = 1) |>
  bold(bold = TRUE, part = "header") |> 
  set_table_properties(layout = "autofit")
```


A sensitivity analysis was done (see to section @sec-lmm-sensitivity for details).


#### multicollinearity of model terms


Multicollinearity of the model terms was assessed with `check_collinerarity()` from the "performance" R package [@performance]
by calculating the (generalized) variance inflation factor (VIF) (@tbl-mm-assumptions-collinearity).


```{r}
#| label: tbl-mm-assumptions-collinearity
#| tbl-cap: "Collinarity check"
#| include: true
#| error: false
#| warning: false
#| message: false

performance::check_collinearity(final) |>
  rename("VIF\n CI.l" = VIF_CI_low,
         "VIF\n CI.h" = VIF_CI_high,
         "SE\n factor" = SE_factor,
         "Tolerance\n CI.l" = Tolerance_CI_low,
         "Tolerance\n CI.h" = Tolerance_CI_high) |>
  flextable() |> 
  colformat_num(digits = 3) |>
  colformat_double(digits = 2) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()
  #set_table_properties(layout = "autofit")
```


For all predictors, the VIF was <5. This indicated a low correlation of that predictor with the other predictors.


#### Normality


The model check for (non-)normality of residuals was done using a Q-Q and density plots  (`check_normality` from R package "performance")
(@fig-mm-assumptions-normality). This command checks the studentized residuals of a mixed model for normal distribution. 

Visual inspection of this plot showed a deviation from normal distribution especially in the lower (i.e., left) tail of the distribution.


```{r}
#| label: fig-mm-assumptions-normality
#| fig-cap: "Check normality of residuals."
#| fig-subcap:
#|   - "Q-Q plot"
#|   - "Density plot"
#| layout-ncol: 2
#| fig-width: 5
#| fig-height: 3
#| fig-align: center
#| include: true
#| error: false
#| warning: false
#| message: false


set.seed(2)
options(ggplot2.discrete.colour = "npg") 

plot(performance::check_normality(final), type = "qq")

plot(performance::check_normality(final), type = "density")
```


#### Uniformity of ranodmized quantile residuals


To test the model for uniformity of randomized quantile residuals, the command `check_residuals()` (R package "performance") was applied to the `final` model. 
This function "*checks generalized linear (mixed) models for uniformity of randomized quantile residuals, 
which can be used to identify typical model misspecification problems, such as over/underdispersion, zero-inflation, 
and residual spatial and temporal autocorrelation.* ... [This is done by simulated quantile residuals.] ... 
*The command then "tests the distribution of the quantile residuals against the uniform distribution using a 
Kolmogorov-Smirnov test.*" (from the help page of this function)


```{r}
#| label: mm-assumptions-residuals
#| error: false
#| include: true
#| warning: false
#| message: false
#| classes: styled-output

final.res <- performance::simulate_residuals(final)

performance::check_residuals(final.res)
```


### Model description


#### Model parameters


Summary of the finally 'accepted' model is shown in @tbl-mm_final-model_sjplot. *P*-values (and degrees of freedom) were calculated using Kenward-Roger's method.


```{r}
#| label: tbl-mm_final-model_sjplot
#| tbl-cap: "Comparision of model summarys of the full `final` model."
#| include: true
#| error: false
#| warning: false
#| message: false
#| results: asis

set.seed(2)
options(ggplot2.discrete.colour = "npg")

sjPlot::tab_model(final, show.reflvl = TRUE, p.val = "kr", df.method = "kr", dv.labels = c("Accepted model"))

```


The model and its results were reported as follows using the `report()` ("reports" R package) applying Kenward-Roger's method for 95% CI calculation.


```{r}
#| label: mm_final-model-report
#| error: false
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 100)

# parameters::standardize_parameters
# standardize_parameters(m, method = "pseudo", ci_method = "satterthwaite") => ci_method: ci_kenward

report::report(final, ci_method = "kenward")
```


#### Model figures

@fig-final-random_effect depicts the plots from the `final` model showing the random estimates and marginal effects next to each other.


```{r}
#| label: fig-final-random_effect
#| fig-cap: 
#|   - "Plots from the `final` model (Kenward-Roger). (a) Forest plot of random effects by subject's `Id`. 
#|   (b) Plot of the predicted values (marginal effects) for `wear_time.visit` by `visit`, `study_arm` and `sex`. "
#| fig-width: 10
#| fig-height: 6
#| fig-align: center
#| error: false
#| warning: false
#| message: false
#| include: true

set.seed(2)
options(ggplot2.discrete.colour = "npg") 

fig.final_random <- sjPlot::plot_model(final, type = "re", colors = "bw", title = "", ci_method = "kenward") + 
  theme_bw() + 
  theme(legend.position = "bottom")
fig.LMM_interaction2 <- sjPlot::plot_model(final, type = "pred", terms = c("visit", "sex", "study_arm"), title = "", ci_method = "kenward") + 
  theme_bw() + 
  theme(legend.position = "bottom")

grid4 <- grid::grid.draw(fig.final_random + fig.LMM_interaction2  + plot_layout(nrow = 1, width = c(1,2))  + plot_annotation(tag_levels = 'a'))

ggsave("Fig_5 panel_LMM_intersection.pdf", plot = grid4, device = cairo_pdf, units = "in", width = 10, height = 5, dpi = 300)
ggsave("Fig_5 panel_LMM_intersection.png", plot = grid4, device = png, units = "in", width = 10, height = 5, dpi = 300)
```


```{r}
#| label: LMM-effect
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

# tbl-cap: "Model data for random effects intercept for each subject and predictions according to visit, sex, and study_arm (Kenward-Roger)."

set.seed(2)

sjPlot::get_model_data(final, type = "pred", terms = c("visit", "sex", "study_arm"), ci_method = "kenward") |>
  as_data_frame() |>
  select(facet, group, x, predicted, std.error, conf.low, conf.high) |>
  rename(Sex = group,
         Group = facet,
         Visit = x,
         Predicted = predicted,
         "Std. error" = std.error,
         Conf.l = conf.low,
         Conf.h = conf.high) |>
  mutate(Group = recode(Group, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  flextable() |> 
  colformat_double(j = c("Predicted", "Std. error", "Conf.l", "Conf.h"), digits = 1) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()
```


Fitting the `final` LMM to @fig-indiv-mean-WT-across-visits resulted in the following @fig-indiv-mean-WT-across-visits-finalmodel:


```{r}
#| label: fig-indiv-mean-WT-across-visits-finalmodel
#| fig-cap: "Mean WT for each individual case across the six follow-up appointments: prediction using fixed effects only and including random effects."
#| error: false
#| warning: false
#| message: false
#| fig-width: 12
#| fig-height: 14
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg") 

theramon_long.lmer <- theramon_long |>
  mutate(final.fixef = predict(final, re.form = NA),   # exclude all random effects (cf. lme4::predict.merMod)
         final.ranef = predict(final, re.form =~(1|Id))) # include all random effects

fig.wear_time.Id_lmm <- theramon_long.lmer |>
  ggplot(aes(x = visit, y = wear_time.visit, group = 1)) +
  geom_point(shape = 21, size = 2, fill = "white") +
  # geom_smooth(method = "lm", se = FALSE, aes(alpha = "OLS"), colour = "black", linewidth = 0.5) +
  geom_line(aes(y = final.fixef, alpha = "fixed effects only"), colour = "blue", linewidth = 0.7, linetype = 1) + 
  geom_line(aes(y = final.ranef, alpha = "fixed and random effects"), colour = "green", linewidth = 0.7, linetype = 1) +
  guides(shape = "none") +
  labs(x = "Wear period", y = "Wear time (h)") + 
  scale_y_continuous(limit = c(0, 30), breaks = seq(0, 30, by = 10), expand = c(0,0)) +
  scale_alpha_manual(name = NULL, values = c(1,1), breaks = c("fixed effects only", "fixed and random effects")) +
  facet_wrap(facet = vars(Id), ncol = 5, labeller = case_names) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  theme(panel.grid.major.y = element_line(colour = "lightgrey", linetype = "dashed"),
        panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank()) +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 30, vjust = 0.8, hjust = 0.8))

fig.wear_time.Id_lmm
```


```{r}
#| label: fig_s3_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_S3 individuelle_Tragezeit_Id-wise_LMM.pdf", plot = fig.wear_time.Id_lmm, device = cairo_pdf, units = "in", width = 12, height = 14, dpi = 300)
ggsave("Fig_S3 individuelle_Tragezeit_Id-wise_LMM.png", plot = fig.wear_time.Id_lmm, device = png, units = "in", width = 12, height = 14, dpi = 300)
```


#### Multiple comparisions {#sec-lmm-comparisons}


Using functions from the "emmeans" package [@emmeans], marginal means (aka Least-Squares Means) and contrasts were calculated 
for specified factors or factor combinations from the `final` model. The *p*-values are based on the t-distribution using 
degres of freedom based on Kenward-Roger's method. Adjustments for multiple comparisons were not applied, instead unadjusted *p*-values were reported.


```{r}
#| label: multitest-emmeams_full1
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

options(width = 100)
set.seed(2)

final.allsets <- emmeans::lsmeans(final, ".", lmer.df = "k")

cat("\n")
cat("*** contrast(final.allsets, revpairwise, lmer.df = k, adjust = none) ***\n")
cat("\n")
emmeans::contrast(final.allsets, "revpairwise", lmer.df = "k", adjust = "none")
cat("\n ... and confidence intervals ... ")
cat("\n")
confint(emmeans::contrast(final.allsets, "revpairwise", lmer.df = "k", adjust = "none"))
cat("******************** \n")
cat("\n")
```


The following summarized the contrasts extracted from the `final' model:


```{r}
#| label: multitest-emmeams_full2
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output # <- only for summary() or console output

set.seed(2)
options(width = 150)

# emmeans(final, ~ visit|study_arm*sex) |> contrast(interaction = c("consec", "consec", "consec"), simple = "each", combine = TRUE, adjust = "none")
# emmeans(final, ~ visit|study_arm*sex) |> contrast(interaction = c("consec", "consec", "consec"), simple = "each", combine = TRUE) |> confint(pairs, reverse = TRUE)

cat("\n")
cat("*** emmeans(final, ~ visit|study_arm*sex) |> contrast(consec, simple = each, combine = TRUE) NOT adjusted! ***\n")
cat("\n")
emmeans::lsmeans(final, ~ visit|study_arm*sex, lmer.df = "k") |> 
  emmeans::contrast("consec", simple = "each", combine = TRUE, adjust = "none")
cat("\n")
cat("\n ... and confidence intervals ... ")
cat("\n")
emmeans::lsmeans(final, ~ visit|study_arm*sex, lmer.df = "k") |> 
  emmeans::contrast("consec", simple = "each", combine = TRUE, adjust = "none") |>
  confint(pairs, reverse = TRUE)
cat("\n")
cat("******************** \n")
cat("\n")
```


### Testing outliers and sensitivity {#sec-lmm-sensitivity}


In @sec-lmm-outliers one observation (#145) was identified as potential outlier (@tbl-outliers_observations). See also @fig-outliers_observations below.

```{r}
#| label: tbl-outliers_observations
#| tbl-cap: "Observations identified as potential outliers or influental observations."
#| error: false
#| warning: false
#| message: false
#| include: true

tbl_outliers <- theramon_long |>
  filter(row == 145) |>
  select(row, Id, study_arm, sex, visit, wear_time.visit) 

tbl_case2 <- theramon_long |>
  filter(Id == "2") |>
  select(row, Id, study_arm, sex, visit, wear_time.visit)
  
tbl_outliers_case2 <- full_join(tbl_outliers, tbl_case2, by = c("row", "Id", "study_arm", "sex", "visit", "wear_time.visit"))

tbl_outliers_case2 |>
  flextable() |> 
  colformat_num(digits = 3) |>
  colformat_double(digits = 2) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()

case2.min <- theramon_long |> filter(Id == "2") |> select(wear_time.visit) |> min()
case2.max <- theramon_long |> filter(Id == "2") |> select(wear_time.visit) |> max()

```


Additionally, subject/case #2 displayed very low wear times at each visit, ranging from `{r} case2.min` hours to `{r} case2.max` hours (@tbl-outliers_observations). 
See also @fig-outliers_observations below.
This case might also be influencing the model.

In this section, we evaluated possible effects of these observations on the linear mixed model in three steps:

1. outlier observation removed from the data set
2. influential case #2 removed from the data set
3. outlier observation and potential influential case #2 removed from the data set

Each time, the model was run on the reduced dataset and the model parameter were then compared.


```{r}
#| label: fig-outliers_observations
#| fig-cap: "Potentially influenting case #2 and outlier observation from case #3."
#| error: false
#| warning: false
#| message: false
#| include: true
#| fig-width: 6
#| fig-height: 3
#| fig-align: center

fig.wear_time.Id_lm_outlier <- theramon_long |> 
  filter(Id == "2" | Id == "3") |>
  arrange(Id, visit) |>
  ggplot(aes(x = visit, y = wear_time.visit, group = 1, color = study_arm, label = row)) +
  geom_line() +
  geom_point(shape = 21, size = 2, fill = "white") +
  geom_label_repel(data = tbl_outliers_case2, aes(label = row), box.padding = 0.35, point.padding = 0.5, color = 'black', segment.color = 'black') +
  # geom_text_repel(aes(label = row), box.padding = 0.35, point.padding = 0.5, segment.color = 'black') +
  guides(shape = "none") + 
  labs(x = "Wear period", y = "Wear time (h)") +
  scale_y_continuous(limit = c(0, 30), breaks = seq(0, 30, by = 10), expand = c(0,0)) +
  scale_colour_discrete(name = "Study arm", labels = c("Group 1 (unaware)", "Group 2 (aware)")) +
  facet_wrap(facet = vars(Id), ncol = 5, labeller = case_names) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  theme(panel.grid.major.y = element_line(color = "lightgrey", linetype = "dashed"), 
        panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank()) +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 30, vjust = 0.8, hjust = 0.8))

fig.wear_time.Id_lm_outlier
```


#### Outlier

To evaluate, if potential outliers (@tbl-outliers_observations) had an impact on model performance, the `final` model was updated by removing 
observation #145 (Id = 3, Group 2, female, Visit 3 = T3-T4) yielding model `final.rowdel`.


```{r}
#| label: mm-outlier-effect
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

set.seed(2)
options(width = 100)

# final.rowdel <- update(final, subset = -c(145))

theramon_long.rowdel <- theramon_long |>
  slice(-145)

final.rowdel <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long.rowdel, REML = TRUE)

summary(final.rowdel, ddf = "Kenward-Roger")
```


#### Influential case

Additionally, the wear times of case #2 were very short in comparison to the other study participants (@tbl-outliers_observations).
What happens, if these were removed from the model? 
The `final` model was updated by deleting the observations related to case #2 from the model yielding model `final.case2`:


```{r}
#| label: mm-final-effect_case2
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

set.seed(2)
options(width = 100)

# final.case2 <- update(final, subset = -c(1, 4, 7, 10, 13, 16))

theramon_long.case2 <- theramon_long |>
  filter(Id != "2")

final.case2 <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long.case2, REML = TRUE)

summary(final.case2, ddf = "Kenward-Roger")
```


#### Combination of outlier and influential case

What happens, if both, the outlier observation and the influential subject, were removed from the data frame? 
This yielded model `final.reduced`.


```{r}
#| label: mm-final-effect_reduced
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

# final.reduced <- update(final, subset = -c(1, 4, 7, 10, 13, 16, 145))

set.seed(2)
options(width = 100)

theramon_long.reduced <- theramon_long.rowdel |>
  filter(Id != "2")

final.reduced <- lmerTest::lmer(wear_time.visit ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long.reduced, REML = TRUE)

theramon.reduced_eff <- Effect(c("visit", "sex", "study_arm"),
                           final.reduced,
                           xlevels = list(visit=seq(1, 6, by=1)))

theramon.full_eff <- Effect(c("visit", "sex", "study_arm"),
                           final,
                           xlevels = list(visit=seq(1, 6, by=1)))

summary(final.reduced, ddf = "Kenward-Roger")
```


Check for outliers:

```{r}
#| label: mm-reduced-outliers
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

check_outliers(final.reduced)
```


#### Comparison


In @tbl-mm-outlier-effect2 we compared the accepted model (`final`) with those after elemination of the potential outlier (`final.rowdel`), 
removal of the potentially influencing case 2 (`final.case2`), or both (`final.reduced`).

* `final.rowdel`: The intercept increased, the estimates for visit(T1-T2) up to visit(T5-T6) increased, the same for sex(Female). 
The estimates for `study_arm` = Group 2 and for the interaction term study_arm(Group2):sex(Female) decreased. 
* `final.case2`: Slight changes in the estimates, sex(Female) and the interaction term study_arm(Group2):sex(Female) 
were not significant anymore and should be probably removed from the analysis?!
* `final.reduced`: The number of observations was reduced from 240 to 233, and the number of groups dropped from 40 to 39.
In comparison to the `final` model (complete data frame) deleting the outlier and the influential subject led to a better 
performance of the `final.reduced` model as indicated by the lower AIC and BIC. 
After adjusting for the other variables in the model, the estimates for each `visit` period slightly increased, and the estimate for `study_arm` was reduced. 
After adjusting for the other variables in the model, the estimate for female `sex` was reduced by `{r} 5.33-3.62`, and was not significant anymore. 
After adjusting for the other variables in the model, the estimate for the interaction of `study_arm`:`sex` for female subjects increased by `{r} -5.54-(-3.72)`, 
and also was not significant anymore.


```{r}
#| label: tbl-mm-outlier-effect2
#| tbl-cap: "Comparision between the accpeted (`final`) and those after elimination of outlier and/or potentially influenting case #2."
#| error: false
#| include: true
#| warning: false
#| message: false
#| results: asis

sjPlot::tab_model(final, final.rowdel, final.case2, final.reduced, p.val = "kr", df.method = "kr",
                  show.reflvl = TRUE, collapse.ci = TRUE, # show.ngroups = TRUE, 
                  # show.se = TRUE, show.ngroups = TRUE, show.aic = TRUE, show.loglik = TRUE, auto.label = FALSE,
                  dv.labels = c("final", "final.rowdel", "final.case2", "final.reduced"))
```


### Dealing with non-normality


Various transformations of the response variable `wear_time.visit` were tested, including `log()`, `log10()`, and `sqrt()` transformations (not shown). 
Additionally, cubed (`wear_time.visit`$^3$) and Box-Cox transformations were tested.

All tests were done with the model for the reduced data set (i.e., after removal of the outlier and the influential case; `final.reduced`) 
and its data set `theramon_long.reduced`.


```{r}
#| label: mm-cubed-transformation
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

set.seed(2)
options(width = 100)

final.reduced_cubed <- lmerTest::lmer(I(wear_time.visit^3) ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long.reduced, REML = TRUE)
summary(final.reduced_cubed, ddf = "Kenward-Roger")

theramon.reduced_cubed_eff <- Effect(c("visit", "sex", "study_arm"),
                                     final.reduced_cubed,
                                     xlevels = list(visit=seq(1, 6, by=1)),
                                     transformation = list(link=function(x) x^3,
                                                           inverse=function(x) x^(1/3)))
```


Finally, Box-Cox transformation was tested as implemented in `powerTransform()` of the "car" R package yielding $\lambda$ = 2.7, $\gamma$ = 22 [@car].
A Box-Cox transformation was applied to `wear_time.visit` in the `theramon_long.reduced` data frame accordingly.


```{r}
#| label: mm-box-cox-transformation
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

set.seed(2)
options(width = 100)

pwr.theramon <- car::powerTransform(final.reduced, family = "bcnPower")

cat("*** powerTransform*** \n")

summary(pwr.theramon)
# bcn - Box-Cox Power transformation to Normality
# allowing for negative values, lmer fit
# 
# Estimated power, lambda
#      Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
# [1,]       2.7     0.84                1              4.3
# 
# Estimated location, gamma
#      Est.gamma Std.Err. Wald Lower Bound Wald Upper Bound
# [1,]        22       14                0               50
# 
# Likelihood ratio tests about transformation parameters
#                       LRT df    pval
# LR test, lambda = (0)  43  1 5.0e-11
# LR test, lambda = (1)  17  1 3.9e-05

cat("\n")

theramon_long.reduced$wear_time.visit_tran <- car::bcnPower(theramon_long.reduced$wear_time.visit, lambda = 2.7, gamma = 22)

final.reduced_pwr <- lmerTest::lmer(wear_time.visit_tran ~ visit + study_arm + sex + study_arm:sex + (1 | Id), data = theramon_long.reduced, REML = TRUE)

cat("*** Box-Cox-transformed model*** \n")

summary(final.reduced_pwr, ddf = "Kenward-Roger")

theramon.reduced_pwr_eff <- Effect(c("visit", "sex", "study_arm"),
                                   final.reduced_pwr,
                                   xlevels = list(visit=seq(1, 6, by=1)),
                                   transformation = list(link=function(x) car::bcnPower(x, lambda = 2.7, gamma = 22),
                                                         inverse=function(x) car::bcnPowerInverse(x, lambda = 2.7, gamma = 22)))
```


```{r}
#| label: fig-mm-transformations-normality
#| fig-cap: "Normality plots of the models after (a) cubed, and (b) Box-Cox transformation"
#| fig-subcap:
#|   - "Cubed transformation"
#|   - "Box-Cox transformation"
#| error: false
#| warning: false
#| message: false
#| include: true
#| layout-ncol: 2
#| fig-width: 5
#| fig-height: 3
#| fig-align: center

plot(check_normality(final.reduced_cubed))

plot(check_normality(final.reduced_pwr))
```


```{r}
#| label: fig-mm-transformations-normality_data
#| fig-cap: "Effect plot for the interaction of `sex` and `study_arm` in the mixed-effect models fit to the reduced data set. Effect of transformations applied to response variable `wear_time.visit`."
#| fig-subcap:
#|   - "Not transformed"
#|   - "Cubed transformation"
#|   - "Box-Cox transformation"
#| error: false
#| warning: false
#| message: false
#| include: true
#| layout-nrow: 3
#| fig-width: 6
#| fig-height: 4
#| fig-align: center

plot(theramon.reduced_eff, axes = list(y=list(type="response"), x = list(rug = FALSE)),
     ylim = c(8, 20),
     lines = list(multiline=TRUE, lwd = 2),
     confit=list(style = "bands"),
     lattice=list(key.args=list(x=0.9, y = 0.2, corner = c(1,1), padding.text=1.25)),
     xlab = "Visit", ylab = "Wear time (h)", main = "")

plot(theramon.reduced_cubed_eff, axes = list(y=list(type="response"), x = list(rug = FALSE)),
     ylim = c(8, 20), 
     lines = list(multiline=TRUE, lwd = 2),
     confit=list(style = "bands"),
     lattice=list(key.args=list(x=0.9, y = 0.2, corner = c(1,1), padding.text=1.25)),
     xlab = "Visit", ylab = "Wear time (h)", main = "")

plot(theramon.reduced_pwr_eff, axes = list(y=list(type="response"), x = list(rug = FALSE)),
     ylim = c(8, 20), 
     lines = list(multiline=TRUE, lwd = 2),
     confit=list(style = "bands"),
     lattice=list(key.args=list(x=0.9, y = 0.2, corner = c(1,1), padding.text=1.25)),
     xlab = "Visit", ylab = "Wear time (h)", main = "")
```


## `final.reduced` model description


```{r}
#| label: mm_final.reduced-model-report
#| echo: true
#| include: true
#| error: false
#| warning: false
#| message: false
#| classes: styled-output

set.seed(2)
options(width = 110)

report(final.reduced, ci_method = "kenward")
```


@tbl-mm_final-reduced-models_sjplot compares `final` and `final.reduced` models.


```{r}
#| label: tbl-mm_final-reduced-models_sjplot
#| tbl-cap: "Comparision of model summarys of the full `final` and the reduced `final.reduced` model using `sjPlot::tab_model() (Kenward-Roger)`"
#| echo: true
#| include: true
#| error: false
#| warning: false
#| message: false
#| results: asis

set.seed(2)

sjPlot::tab_model(final, final.reduced, show.reflvl = TRUE, p.val = "kr", df.method = "kr",
                  # show.se = TRUE, show.ngroups = TRUE, show.aic = TRUE, show.loglik = TRUE, auto.label = FALSE,
                  dv.labels = c("Accepted model", "Reduced model"))
```


@fig-final_reduced-random_effect depicts the plots from the `final.reduced` model showing the random estimates and marginal effects next to each other.


```{r}
#| label: fig-final_reduced-random_effect
#| fig-cap: 
#|   - "Plots from the `final.reduced` model. (a) Forest plot of random effects by subject's `Id`. 
#|   (b) Plot of the predicted values (marginal effects) for `wear_time.visit` by `visit`, `study_arm` and `sex`. "
#| fig-width: 10
#| fig-height: 6
#| fig-align: center
#| error: false
#| warning: false
#| message: false
#| include: true

set.seed(2)
options(ggplot2.discrete.colour = "npg") 

fig.final.reduced_random <- sjPlot::plot_model(final.reduced, type = "re", colors = "bw", title = "", ci_method = "kenward") + 
  theme_bw() + theme(legend.position = "bottom")
fig.final.reduced_interaction2 <- sjPlot::plot_model(final.reduced, type = "pred", terms = c("visit", "sex", "study_arm"), title = "", ci_method = "kenward") + 
  theme_bw() + theme(legend.position = "bottom")

grid5 <- grid::grid.draw(fig.final.reduced_random + fig.final.reduced_interaction2  + plot_layout(nrow = 1, width = c(1,2))  + plot_annotation(tag_levels = 'a'))

ggsave("Fig_6 panel_final-reduced_intersection.pdf", plot = grid5, device = cairo_pdf, units = "in", width = 10, height = 7, dpi = 300)
ggsave("Fig_6 panel_final-reduced_intersection.png", plot = grid5, device = png, units = "in", width = 10, height = 7, dpi = 300)

```

### Multiple comparisons


```{r}
#| label: LMM-LS-means_finalreduced
#| error: false
#| warning: false
#| message: false
#| include: false
#| echo: false
#| eval: false

# "Least Squares Means table of the reduced model `final.reduced` (Kenward-Roger)."

finalreduced.ls_means <- ls_means(final.reduced, ddf = "Kenward-Roger")

finalreduced.ls_means |> 
  as_tibble() |>
  rename("p" = "Pr(>|t|)") |>
  # [1] "term"       "levels"     "Estimate"   "Std. Error" "df"         "t value"    "lower"      "upper"      "Pr(>|t|)"  
  flextable() |> 
  set_formatter(values = list("p" = p_val_format) ) |>
  colformat_double(j = c("Estimate", "df", "lower", "upper"), digits = 1) |>
  colformat_double(j = c("Std. Error", "t value"), digits = 2) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()

# "Pairwise comparison table of the reduced model `final.reduced` (Kenward-Roger)."

finalreduced.difflsmeans <- difflsmeans(final.reduced, ddf = "Kenward-Roger")

finalreduced.difflsmeans |> 
  as_tibble() |>
  rename("p" = "Pr(>|t|)") |>
  # [1] "term"       "levels"     "Estimate"   "Std. Error" "df"         "t value"    "lower"      "upper"      "Pr(>|t|)"  
  flextable() |> 
  set_formatter(values = list("p" = p_val_format) ) |>
  colformat_double(j = c("df"), digits = 1) |>
  colformat_double(j = c("Estimate", "Std. Error", "t value", "lower", "upper"), digits = 2) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()
```


Multiple comparisons of the `final.reduced` model factors was done as described above (see @sec-lmm-comparisons). 
The *p*-values are based on the t-distribution using degrees of freedom based on Kenward-Roger's method. 
Adjustments for multiple comparisons were not applied, instead unadjusted *p*-values were reported.


```{r}
#| label: test-emmeams_reduced1
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output

options(width = 100)
set.seed(2)

reduced.allsets <- emmeans::lsmeans(final.reduced, ".", lmer.df = "k")

cat("\n")
cat("*** contrast(reduced.allsets, revpairwise, lmer.df = k, adjust = none) ***\n")
cat("\n")
emmeans::contrast(reduced.allsets, "revpairwise", lmer.df = "k", adjust = "none")
cat("\n ... and confidence intervals ... ")
cat("\n")
confint(emmeans::contrast(reduced.allsets, "revpairwise", lmer.df = "k", adjust = "none"))
cat("******************** \n")
cat("\n")
```


The following summarized the contrasts extracted from the `final.reduced' model:


```{r}
#| label: test-emmeams_reduced2
#| error: false
#| warning: false
#| message: false
#| include: true
#| classes: styled-output # <- only for summary() or console output

set.seed(2)
options(width = 150)

# emmeans(final, ~ visit|study_arm*sex) |> contrast(interaction = c("consec", "consec", "consec"), simple = "each", combine = TRUE, adjust = "none")
# emmeans(final, ~ visit|study_arm*sex) |> contrast(interaction = c("consec", "consec", "consec"), simple = "each", combine = TRUE) |> confint(pairs, reverse = TRUE)

cat("\n")
cat("*** lsmeans(final.reduced, ~ visit|study_arm*sex) |> contrast(consec, simple = each, combine = TRUE) NOT adjusted! ***\n")
cat("\n")
emmeans::lsmeans(final.reduced, ~ visit|study_arm*sex, lmer.df = "k") |> 
  emmeans::contrast("consec", simple = "each", combine = TRUE, adjust = "none")
cat("\n")
cat("\n ... and confidence intervals ... ")
cat("\n")
emmeans::lsmeans(final.reduced, ~ visit|study_arm*sex, lmer.df = "k") |> 
  emmeans::contrast("consec", simple = "each", combine = TRUE, adjust = "none") |>
  confint(pairs, reverse = TRUE)
cat("\n")
cat("******************** \n")
cat("\n")
```


### Figures


```{r}
#| label: LMM-effect_reduced_model
#| error: false
#| warning: false
#| message: false
#| include: true
#| eval: false

# tbl-cap: "Predictions of wear time related to visit, sex, and study_arm for the reduced model (Kenward-Roger)."

set.seed(2)

get_model_data(final.reduced, type = "pred", terms = c("visit", "sex", "study_arm"), ci_method = "kenward") |>
  as_data_frame() |>
  select(facet, group, x, predicted, std.error, conf.low, conf.high) |>
  rename(Sex = group,
         Group = facet,
         Visit = x,
         Predicted = predicted,
         "Std. error" = std.error,
         Conf.l = conf.low,
         Conf.h = conf.high) |>
  mutate(Group = recode(Group, 'Group 1 (unaware)' = 'Group 1', 'Group 2 (aware)' = 'Group 2')) |>
  flextable() |> 
  colformat_double(j = c("Predicted", "Std. error", "Conf.l", "Conf.h"), digits = 1) |>
  bold(bold = TRUE, part = "header") |> 
  autofit()
```


Fitting the `final.reduced` LMM to @fig-indiv-mean-WT-across-visits resulted in the following @fig-indiv-mean-WT-across-visits-final_reduced:


```{r}
#| label: fig-indiv-mean-WT-across-visits-final_reduced
#| fig-cap: "Wear time of each individual case across the six follow-up appointments using the 'reduced' model: prediction using fixed effects only and including random effects."
#| echo: true
#| error: false
#| warning: false
#| message: false
#| include: true
#| fig-width: 12
#| fig-height: 14
#| fig-align: center

set.seed(2)
options(ggplot2.discrete.colour = "npg") 

theramon_long.reduced_lmer <- theramon_long.reduced |>
  mutate(final_reduced.fixef = predict(final.reduced, re.form = NA),   # exclude all random effects (cf. lme4::predict.merMod)
         final_reduced.ranef = predict(final.reduced, re.form =~(1|Id))) # include all random effects

theramon_long.joined_lmer <- left_join(theramon_long.lmer, theramon_long.reduced_lmer, by = c("Id", "study_arm", "visit", "wear_time.visit")) |>
  select(Id, study_arm, age.x, sex.x, visit, visit_n.x, wear_time.visit, row.x, cum_days.x, final.fixef, final.ranef, final_reduced.fixef, final_reduced.ranef) |>
  rename(age = age.x,
         sex = sex.x,
         visit_n = visit_n.x,
         row = row.x,
         cum_days = cum_days.x)

fig.wear_time.joined_Id_lmm <- theramon_long.joined_lmer |>
  ggplot(aes(x = visit, y = wear_time.visit, group = 1)) +
  geom_point(shape = 21, size = 2, fill = "white") +
  # geom_smooth(method = "lm", se = FALSE, aes(alpha = "OLS"), colour = "black", linewidth = 0.5) +
  geom_line(aes(y = final_reduced.fixef, alpha = "fixed effects only"), colour = "blue", linewidth = 0.7, linetype = 1) + 
  geom_line(aes(y = final_reduced.ranef, alpha = "fixed and random effects"), colour = "green", linewidth = 0.7, linetype = 1) +
  guides(shape = "none") +
  labs(x = "Wear period", y = "Wear time (h)") + 
  scale_y_continuous(limit = c(0, 30), breaks = seq(0, 30, by = 10), expand = c(0,0)) +
  scale_alpha_manual(name = NULL, values = c(1,1), breaks = c("fixed effects only", "fixed and random effects")) +
  facet_wrap(facet = vars(Id), ncol = 5, labeller = case_names) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  theme(panel.grid.major.y = element_line(colour = "lightgrey", linetype = "dashed"),
        panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank()) +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 30, vjust = 0.8, hjust = 0.8))

fig.wear_time.joined_Id_lmm
```



```{r}
#| label: fig_sX_export
#| error: false
#| warning: false
#| message: false
#| include: false
#| eval: false

ggsave("Fig_SX individuelle_Tragezeit_reduced_Id-wise_LMM.pdf", plot = fig.wear_time.joined_Id_lmm, device = cairo_pdf, units = "in", width = 12, height = 14, dpi = 300)
ggsave("Fig_SX individuelle_Tragezeit_reduced_Id-wise_LMM.png", plot = fig.wear_time.joined_Id_lmm, device = png, units = "in", width = 12, height = 14, dpi = 300)
```


# Citations


```{r}
#| label: citations
#| echo: true
#| include: true
#| error: false
#| warning: false
#| message: false

pkgs <- grateful::cite_packages(output = "table", out.dir = ".", cite.tidyverse = TRUE)
knitr::kable(pkgs)
```


```{r}
#| label: rdata
#| echo: false
#| error: false
#| warning: false
#| message: false

save.image("./Theramon_data-supplement_accepted.RData")
```


# Sessioninfo


```{r}
#| label: sessioninfo
#| echo: false
#| error: false
#| warning: false
#| message: false
#| include: true

sessionInfo()
```
