Introduction and methodology

d <- read_rds("BrainFlow.RDS") %>% 
      select(id, age, sex, af_statusf, measure, tCBFPC, visitf) %>% 
      set_names(c("id", "age", "sex", "status", "measure", "tcbfp", "visit")) %>% 
      mutate(sex = fct_recode(factor(sex),
                              "Male" = "1",
                              "Female" = "2"),
             visit = fct_recode(visit,
                                "Pre" = "1",
                                "Post" = "2"))

In this experiment participant levels of tCBFP were measured before and after a treatment. Participants were grouped by whether the treatment was deemed successful (SR) or not (AF). A linear mixed model (LMM) was used to estimate the treatment effect as well as any difference in effect between the two groups. The LMM was fitted according to the formula

\[ Y_{ij} = \beta_0 + \beta_1 X_i + \beta_2 t_{ij} + \beta_3 t_{ij} X_i + \epsilon_{ij}, \]

where \(Y_{ij}\) is the j’th measurement for subject i, \(X_i\) is the treatment group of subject i, and \(t_ij\) indicates whether \(Y_{ij}\) is a pre- or post measurement for subject i.

Summary statistics

cap <- "Sample summary statistics."
d %>% 
   select(id, status, visit, tcbfp) %>% 
   spread(visit, tcbfp) %>% 
   mutate(change = Post - Pre) %>% 
   group_by(status) %>% 
   summarize(n = n(),
             Change = paste0(round(mean(change, na.rm = T), 3),
                             " (", round(sd(change, na.rm = T), 3), ")"),
             Pre = paste0(round(mean(Pre, na.rm = T), 3),
                          " (", round(sd(Pre, na.rm = T), 3), ")"),
             Post = paste0(round(mean(Post, na.rm = T), 3),
                           " (", round(sd(Post, na.rm = T), 3), ")")) %>% 
   mutate(status = fct_recode(status,
                              "SR" = "SR Group",
                              "AF" = "AF Group")) %>% 
   rename(Group = status) %>% 
   kable(caption = table_cap(cap),
         digits = 3,
         align = c("l", rep("c", 4))) %>% 
   kable_styling() %>% 
   add_header_above(c("", "", "Mean (SD)" = 3))

Pre post change

m <- d %>% 
   lmer(tcbfp ~ visit + status + visit : status + (1 | id): status, data = .)

cap <- "Significance test for overall change in tCBFP."
emmeans(m, pairwise ~ visit | status)$contrasts %>% 
   tidy %>% 
   unite(Contrast, level1, level2) %>% 
   mutate(Contrast = "Post - Pre",
          estimate = -estimate,
          upper = round(estimate + qt(0.975, df) * std.error, 3),
          lower = round(estimate - qt(0.975, df) * std.error, 3)) %>% 
   select(Contrast, status, estimate, lower, upper, p.value) %>% 
   set_names(c("Contrast", "Group", "Estimate", "Lower", "Upper", "p")) %>% 
   kable(caption = table_cap(cap),
         digits = 3,
         align = c("l", rep("c", 5))) %>% 
   kable_styling() %>% 
   add_header_above(c("", "", "", "95% Confidence Interval" = 2, ""))

cap <- "Significance test for difference in change between groups."
tidy(m, effects = "fixed", conf.int = TRUE) %>% 
   filter(str_detect(term, ":")) %>% 
   select(term, estimate, conf.low, conf.high, p.value) %>% 
   mutate(term = "AF - SR") %>% 
   set_names(c("Contrast", "Estimate", "Lower", "Upper", "p")) %>% 
   kable(caption = table_cap(cap),
         digits = 3,
         align = c("l", rep("c", 4))) %>% 
   kable_styling() %>% 
   add_header_above(c("", "", "95% Confidence Interval" = 2, ""))

Observed and fitted trajectories

Figure 1 below shows observed and fitted trajectories for each participant. The trend for the SR group seems to be positive but there is still uncertainty.

cap <- "Observed and fitted trajectories including mean-valued dashed lines for successful and unsuccessful treatments"

d %>% 
   (function(x) {
      x$pred <- predict(m, x)
      x
   }) %>% 
   rename(Predicted = pred, Observed = tcbfp) %>% 
   gather(type, value, Observed, Predicted) %>% 
   ggplot(aes(visit, value)) +
   geom_line(aes(group = id, col = status), alpha = 0.4) +
   stat_summary(geom = "line", 
                aes(group = status, col = status), 
                lty = 2,
                fun.data = mean_se) +
   facet_grid(type ~ status) +
   theme(legend.position = "top") +
   scale_color_jama(guide = "none") +
   scale_fill_jama(guide = "none") +
   labs(x = "Visit",
        y = "tCBFP")

Figure 1: Observed and fitted trajectories including mean-valued dashed lines for successful and unsuccessful treatments

Residuals

Figures 2 and 3 shows that model residuals follow a normal distribution as expected but there are signs of uneven variance conditioned on fitted values.

cap <- "Distribution of LMM residuals."
augment(m) %>% 
      ggplot(aes(.resid)) +
      geom_density() +
      labs(x = "Residual",
           y = "Density")

Figure 2: Distribution of LMM residuals.

cap <- "Distribution of LMM residuals conditioned on predicted value."
augment(m) %>% 
      ggplot(aes(.fitted, .resid)) +
      geom_hline(yintercept = 0, lty = 2, col = "grey30") +
      geom_smooth() +
      geom_point() +
      labs(x = "Fitted",
           y = "Residuals")

Figure 3: Distribution of LMM residuals conditioned on predicted value.