##Lancet paper : Once-weekly semaglutide versus placebo in patients with alcohol use disorder and comorbid obesity: a randomised, double-blind, placebo-controlled trial ##Creating dataset

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
set.seed(777)

#ASDL Subject level data

n <- 180

adsl <- tibble(
  USUBJID = paste0("SUBJ-", sprintf("%03d", 1:n)),
  
  TRT = sample(c("SEMAGLUTIDE", "PLACEBO"), n, replace = TRUE),
  
  AGE = round(rnorm(n, 45, 10)),
  
  SEX = sample(c("M", "F"), n, replace = TRUE),
  
  BMI = round(rnorm(n, 35, 5), 1),  # obese population
  
  BASE_DRINKS = round(rnorm(n, 25, 8)),  # drinks/week
  
  BASE_WEIGHT = round(rnorm(n, 105, 15), 1)
)

##Longitudinal efficacy dataset

visits <- c("BASE", "W4", "W8", "W12")

adeff <- expand.grid(
  USUBJID = adsl$USUBJID,
  VISIT = visits
) %>%
  as_tibble() %>%
  
  left_join(adsl, by = "USUBJID") %>%
  
  mutate(
    
    DRINKS = case_when(
      
      VISIT == "BASE" ~ BASE_DRINKS,
      
      TRT == "SEMAGLUTIDE" & VISIT == "W4" ~ BASE_DRINKS - rnorm(n(), 6, 3),
      TRT == "SEMAGLUTIDE" & VISIT == "W8" ~ BASE_DRINKS - rnorm(n(), 10, 4),
      TRT == "SEMAGLUTIDE" & VISIT == "W12" ~ BASE_DRINKS - rnorm(n(), 14, 5),
      
      TRT == "PLACEBO" & VISIT == "W4" ~ BASE_DRINKS - rnorm(n(), 2, 3),
      TRT == "PLACEBO" & VISIT == "W8" ~ BASE_DRINKS - rnorm(n(), 4, 4),
      TRT == "PLACEBO" & VISIT == "W12" ~ BASE_DRINKS - rnorm(n(), 5, 5)
    ),
    
    WEIGHT = case_when(
      VISIT == "BASE" ~ BASE_WEIGHT,
      
      TRT == "SEMAGLUTIDE" ~ BASE_WEIGHT - rnorm(n(), 5, 2),
      TRUE ~ BASE_WEIGHT - rnorm(n(), 1, 2)
    )
  )

##Column drinks above is the number of drinks at that visit. ##Deriving change from baseline

adeff <- adeff %>%
  
  group_by(USUBJID) %>%
  
  mutate(
    CHG_DRINKS = DRINKS - DRINKS[VISIT == "BASE"][1],
    CHG_WEIGHT = WEIGHT - WEIGHT[VISIT == "BASE"][1]
  ) %>%
  
  ungroup()

##In real world scenarios, baseline values are confused. So just reinforcing the baseline.

adeff <- adeff %>%
  group_by(USUBJID) %>%
  mutate(
    BASE_DRINKS2 = first(DRINKS)
  ) %>%
  ungroup()

##Primary analysis ##Mean change at week 12

week12 <- adeff %>%
  filter(VISIT == "W12")

week12 %>%
  group_by(TRT) %>%
  summarise(
    Mean_Change = mean(CHG_DRINKS),
    SD = sd(CHG_DRINKS)
  )
## # A tibble: 2 × 3
##   TRT         Mean_Change    SD
##   <chr>             <dbl> <dbl>
## 1 PLACEBO           -4.56  4.82
## 2 SEMAGLUTIDE      -13.1   5.61

##Similarly mean change in weight

wk12 <- 
  adeff |> 
  filter(VISIT == "W12")
wk12 |> 
  group_by(TRT) |> 
  summarise(
    Mean_change = mean(CHG_WEIGHT),
      SD = sd(CHG_WEIGHT)
  )
## # A tibble: 2 × 3
##   TRT         Mean_change    SD
##   <chr>             <dbl> <dbl>
## 1 PLACEBO           -1.19  1.89
## 2 SEMAGLUTIDE       -5.40  2.06

##Comparison of means

t.test(CHG_DRINKS ~ TRT, data = week12)
## 
##  Welch Two Sample t-test
## 
## data:  CHG_DRINKS by TRT
## t = 10.852, df = 158.65, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group PLACEBO and group SEMAGLUTIDE is not equal to 0
## 95 percent confidence interval:
##   7.006352 10.124194
## sample estimates:
##     mean in group PLACEBO mean in group SEMAGLUTIDE 
##                  -4.55774                 -13.12301

##Visualization

adeff %>%
  group_by(VISIT, TRT) %>%
  summarise(mean_drinks = mean(DRINKS)) %>%
  ggplot(aes(x = VISIT, y = mean_drinks, color = TRT, group = TRT)) +
  geom_line() +
  geom_point()
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by VISIT and TRT.
## ℹ Output is grouped by VISIT.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(VISIT, TRT))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

##Safety dataset

ae <- tibble(
  USUBJID = sample(adsl$USUBJID, 250, replace = TRUE),
  
  AEDECOD = sample(c("NAUSEA", "VOMITING", "HEADACHE"), 250, replace = TRUE),
  
  AESEV = sample(c("MILD", "MODERATE", "SEVERE"), 250, replace = TRUE)
)

##Now to summarise the safety dataset

ae %>%
  count(AEDECOD, AESEV)
## # A tibble: 9 × 3
##   AEDECOD  AESEV        n
##   <chr>    <chr>    <int>
## 1 HEADACHE MILD        31
## 2 HEADACHE MODERATE    28
## 3 HEADACHE SEVERE      21
## 4 NAUSEA   MILD        26
## 5 NAUSEA   MODERATE    29
## 6 NAUSEA   SEVERE      23
## 7 VOMITING MILD        23
## 8 VOMITING MODERATE    32
## 9 VOMITING SEVERE      37

##Simulating dropouts

set.seed(999)

adeff <- adeff %>%
  group_by(USUBJID) %>%
  mutate(
    dropout = sample(c(0,1), 1, prob = c(0.8, 0.2)),
    VISITNUM = case_when(
      VISIT == "BASE" ~ 0,
      VISIT == "W4" ~ 4,
      VISIT == "W8" ~ 8,
      VISIT == "W12" ~ 12
    )
  ) %>%
  mutate(
    DRINKS = ifelse(dropout == 1 & VISITNUM > 4, NA, DRINKS)
  ) %>%
  ungroup()

##Creating a PARAM structure

adparam <- adeff %>%
  mutate(
    PARAM = "Weekly Alcohol Consumption",
    PARAMCD = "DRINKS",
    AVAL = DRINKS
  ) %>%
  select(USUBJID, TRT, VISIT, VISITNUM, PARAM, PARAMCD, AVAL)

##Adding baseline and change

adparam <- adparam %>%
  group_by(USUBJID, PARAMCD) %>%
  arrange(VISITNUM) %>%
  mutate(
    BASE = first(AVAL),
    CHG = AVAL - BASE
  ) %>%
  ungroup()

##Create w12 data

w12 <- adparam %>%
  filter(VISIT == "W12")

##Run ANCOVA

ancova <- lm(CHG ~ TRT + BASE, data = w12)
summary(ancova)
## 
## Call:
## lm(formula = CHG ~ TRT + BASE, data = w12)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.2019  -3.4337  -0.6076   3.7704  11.2994 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -5.97631    1.51884  -3.935 0.000131 ***
## TRTSEMAGLUTIDE -8.81816    0.87551 -10.072  < 2e-16 ***
## BASE            0.04646    0.05614   0.828 0.409333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.146 on 139 degrees of freedom
##   (38 observations deleted due to missingness)
## Multiple R-squared:  0.4223, Adjusted R-squared:  0.414 
## F-statistic:  50.8 on 2 and 139 DF,  p-value: < 2.2e-16

##Repeated measures model

mmrm <- nlme::lme(
  CHG ~ TRT * VISIT,
  random = ~1 | USUBJID,
  data = adparam,
  na.action = na.exclude
)

summary(mmrm)
## Linear mixed-effects model fit by REML
##   Data: adparam 
##        AIC      BIC    logLik
##   3407.593 3452.145 -1693.796
## 
## Random effects:
##  Formula: ~1 | USUBJID
##          (Intercept) Residual
## StdDev: 0.0002712695 3.376178
## 
## Fixed effects:  CHG ~ TRT * VISIT 
##                             Value Std.Error  DF    t-value p-value
## (Intercept)              0.000000 0.3393187 458   0.000000       1
## TRTSEMAGLUTIDE           0.000000 0.5058264 178   0.000000       1
## VISITW4                 -2.206709 0.4798691 458  -4.598565       0
## VISITW8                 -4.720130 0.5093358 458  -9.267224       0
## VISITW12                -4.814247 0.5093358 458  -9.452008       0
## TRTSEMAGLUTIDE:VISITW4  -3.759950 0.7153466 458  -5.256124       0
## TRTSEMAGLUTIDE:VISITW8  -4.903423 0.7622834 458  -6.432546       0
## TRTSEMAGLUTIDE:VISITW12 -8.730995 0.7622834 458 -11.453739       0
##  Correlation: 
##                         (Intr) TRTSEMAGLUTIDE VISITW4 VISITW8 VISITW1
## TRTSEMAGLUTIDE          -0.671                                       
## VISITW4                 -0.707  0.474                                
## VISITW8                 -0.666  0.447          0.471                 
## VISITW12                -0.666  0.447          0.471   0.444         
## TRTSEMAGLUTIDE:VISITW4   0.474 -0.707         -0.671  -0.316  -0.316 
## TRTSEMAGLUTIDE:VISITW8   0.445 -0.664         -0.315  -0.668  -0.297 
## TRTSEMAGLUTIDE:VISITW12  0.445 -0.664         -0.315  -0.297  -0.668 
##                         TRTSEMAGLUTIDE:VISITW4 TRTSEMAGLUTIDE:VISITW8
## TRTSEMAGLUTIDE                                                       
## VISITW4                                                              
## VISITW8                                                              
## VISITW12                                                             
## TRTSEMAGLUTIDE:VISITW4                                               
## TRTSEMAGLUTIDE:VISITW8   0.469                                       
## TRTSEMAGLUTIDE:VISITW12  0.469                  0.440                
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -4.073919e+00 -4.571239e-01 -9.297148e-10  4.225142e-01  3.389614e+00 
## 
## Number of Observations: 644
## Number of Groups: 180

##Make efficacy table

w12 %>%
  group_by(TRT) %>%
  summarise(
    N = n(),
    Mean = mean(CHG, na.rm=TRUE),
    SD = sd(CHG, na.rm=TRUE)
  )
## # A tibble: 2 × 4
##   TRT             N   Mean    SD
##   <chr>       <int>  <dbl> <dbl>
## 1 PLACEBO        99  -4.81  4.79
## 2 SEMAGLUTIDE    81 -13.5   5.54