marginalizing-gam.R

``` r
# from https://github.com/ASKurz/non-linear-RCT/blob/master/Simulate%20the%20data.Rmd

library(tidyverse)
library(faux)
#> 
#> ************
#> Welcome to faux. For support and examples visit:
#> https://debruine.github.io/faux/
#> - Get and set global package options with: faux_options()
#> ************
#> 
#> Attaching package: 'faux'
#> The following object is masked from 'package:purrr':
#> 
#>     %||%

# initial sample size per group
n <- 54
# simulate MSEL and VABS values for the A/M group
set.seed(1)
pe <- rnorm_multi(
  n = n,
  mu = c(78.28, 65.03, 44.28, 36.63, 38.33),
  sd = c(16.35, 29.19, 33.73, 31.80, 28.49), 
  r = .5, 
  varnames = list("caps.0", "caps.1", "caps.2", "caps.3", "caps.4")
)

# initial sample size per group
n <- 54
# simulate MSEL and VABS values for the A/M group
set.seed(2)
vret <- rnorm_multi(
  n = n,
  mu = c(80.44, 71.19, 57.07, 56.64, 53.50),
  sd = c(16.23, 23.27, 32.32, 31.50, 28.07), 
  r = .5, 
  varnames = list("caps.0", "caps.1", "caps.2", "caps.3", "caps.4")
)

d <- bind_rows(pe, vret) %>% 
  mutate(id = 1:n(),
         tx = rep(c("pe", "vret"), each = n() / 2)) %>% 
  pivot_longer(contains("caps"), values_to = "caps") %>% 
  mutate(wave = str_remove(name, "caps.") %>% as.double()) %>% 
  mutate(caps = ifelse(caps < 0, 0,
                       ifelse(caps > 123, 123, caps))) %>% 
  mutate(weeks = case_when(
    wave == 0 ~ 0,
    wave == 1 ~ 2.5,
    wave == 2 ~ 5,
    wave == 3 ~ 12,
    wave == 4 ~ 26
  )) %>% 
  mutate(group = ifelse(tx == "pe", "a", "b"))


d %>% 
  ggplot(aes(x = weeks)) +
  geom_line(aes(y = caps, group = id),
            size = 1/6, alpha = 1/2) +
  geom_line(data = . %>% 
              group_by(weeks, group) %>% 
              summarise(caps = mean(caps)),
            aes(y = caps),
            size = 1.5) +
  geom_linerange(data = . %>% 
                   group_by(weeks, group) %>% 
                   summarise(m = mean(caps),
                             s = sd(caps)),
                 aes(y = m, ymin = m - s, ymax = m + s)
  ) +
  scale_x_continuous(breaks = c(0, 2.5, 5, 12, 26),
                     labels = c("0", "2.5", "5", "12", "26")) +
  labs(subtitle = "Dark lines are the group means and the vertical lines are plus/minus one standard deviation.\nThe thinner lines are the participant-level trajectories.",
       y = "dv") +
  facet_wrap(~ group, labeller = label_both) +
  theme(panel.grid = element_blank())
#> `summarise()` has grouped output by 'weeks'. You can override using the
#> `.groups` argument.
#> `summarise()` has grouped output by 'weeks'. You can override using the
#> `.groups` argument.
```

![](https://i.imgur.com/yjukGoC.png)

``` r
d %>% 
  mutate(wave = factor(wave,))
#> # A tibble: 540 × 7
#>       id tx    name    caps wave  weeks group
#>    <int> <chr> <chr>  <dbl> <fct> <dbl> <chr>
#>  1     1 pe    caps.0  81.2 0       0   a    
#>  2     1 pe    caps.1  66.1 1       2.5 a    
#>  3     1 pe    caps.2  84.8 2       5   a    
#>  4     1 pe    caps.3  26.5 3      12   a    
#>  5     1 pe    caps.4  68.1 4      26   a    
#>  6     2 pe    caps.0  61.7 0       0   a    
#>  7     2 pe    caps.1  66.8 1       2.5 a    
#>  8     2 pe    caps.2  66.3 2       5   a    
#>  9     2 pe    caps.3   0   3      12   a    
#> 10     2 pe    caps.4  58.7 4      26   a    
#> # … with 530 more rows
ggplot(aes(x = wave, y = caps)) +
  geom_boxplot() +
  facet_wrap(~ tx)
#> Error in `fortify()`:
#> ! `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class uneval.
#> Did you accidentally pass `aes()` to the `data` argument?

##### new code starts here -----------------------------------------------------

library(mgcv)
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> 
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> 
#> This is mgcv 1.8-40. For overview type 'help("mgcv-package")'.
library(gratia)

data <- d |> 
  mutate_at(vars(id, group, wave, tx, name), as.factor)

model <- gam(
  caps ~ group + s(weeks, by = group, k = 4) + s(weeks, id, bs = 'fs', k = 4),
  data = data,
  method = "REML",
  family = "gaussian"
)

grid <- data |> 
  distinct(group) |> 
  expand_grid(
    weeks = 1:30,
    id = factor(1) # doesn't matter since we exclude participant later, but errors without this
  )

grid
#> # A tibble: 60 × 3
#>    group weeks id   
#>    <fct> <int> <fct>
#>  1 a         1 1    
#>  2 a         2 1    
#>  3 a         3 1    
#>  4 a         4 1    
#>  5 a         5 1    
#>  6 a         6 1    
#>  7 a         7 1    
#>  8 a         8 1    
#>  9 a         9 1    
#> 10 a        10 1    
#> # … with 50 more rows

conditional_differences <- fitted_samples(
  model,
  newdata = grid,
  exclude = "s(id)",
  n = 10,
  seed = 27
) |> 
  left_join(
    mutate(grid, row = row_number()),
    by = "row"
  ) |> 
  pivot_wider(
    id_cols = c("draw", "weeks"),
    names_from = "group",
    values_from = "fitted"
  ) |> 
  mutate(
    diff = a - b
  ) |> 
  group_by(weeks) |> 
  summarize(
    lower = quantile(diff, 0.025),
    avg = mean(diff),
    upper = quantile(diff, 0.975), 
    .groups = "drop"
  ) |> 
  select(weeks, avg, lower, upper)

conditional_differences |> 
  ggplot(aes(x = weeks, y = avg, ymin = lower, ymax = upper)) +
  geom_line() +
  geom_ribbon(alpha = 0.3) +
  theme_classic()
```

![](https://i.imgur.com/vNLSnmn.png)

<sup>Created on 2022-05-10 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)</sup>