##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