--- title: "

Visualization of DLQI Mixed Models Results by Agustin Calatroni" output: html_document: self_containded: TRUE code_download: yes toc: false --- ```{r setup, include=FALSE} options(width = 200) knitr::opts_chunk$set(echo = FALSE, eval = TRUE, warning = FALSE, message = FALSE, comment = NA, cache = FALSE) ``` ```{r packages} pacman::p_load(rio, tidyverse) pacman::p_load(lme4, emmeans) pacman::p_load(ggdist, distributional) pacman::p_load(gtsummary) #devtools::install_github("teunbrand/ggh4x") library(ggh4x) pacman::p_load(ggdendro) pacman::p_load(ggtext) pacman::p_load(patchwork) pacman::p_load(gt) ``` ```{r import data} df1 <- import("https://raw.githubusercontent.com/VIS-SIG/Wonderful-Wednesdays/master/data/2021/2021-01-13/ww2020_dlqi.csv") %>% drop_na(DLQI_SCORE) %>% mutate(VISIT_i = as.integer(as.factor(VISIT))) export(df1, "dlqi_mixed.csv") ``` ```{r eval=FALSE, include=FALSE} # fun code to derive change df1_chg <- import("https://raw.githubusercontent.com/VIS-SIG/Wonderful-Wednesdays/master/data/2021/2021-01-13/ww2020_dlqi.csv") %>% select(USUBJID, VISIT, TRT, starts_with("DLQI")) %>% arrange(USUBJID, VISIT) %>% group_by(USUBJID, TRT) %>% # filter(n() > 1) %>% # keep participants with 2 visits transmute(across(starts_with("DLQI"), ~.x - lag(.x), .names = "{.col}_CHG")) %>% rowwise() %>% filter( sum(is.na(c_across(starts_with("DLQI")))) < 11 ) ``` ```{r variable labels} dl <- c(DLQI101 = 'How Itchy, Sore, Painful, Stinging

DLQI101', DLQI102 = 'How Embarrassed, Self Conscious

DLQI102', DLQI103 = 'Interfered Shopping, Home, Yard

DLQI103', DLQI104 = 'Influenced Clothes You Wear

DLQI104', DLQI105 = 'Affected Social, Leisure Activity

DLQI105', DLQI106 = 'Made It Difficult to Do Any Sports

DLQI106', DLQI107 = 'Prevented Working or Studying

DLQI107', DLQI108 = 'Problem Partner, Friends, Relative

DLQI108', DLQI109 = 'Caused Any Sexual Difficulties

DLQI109', DLQI110 = 'How Much a Problem is Treatment

DLQI110', DLQI_SCORE = '**DLQI Total Score**

DLQI_SCORE') cl <- c(MEAN = "TRT by VISIT", DIFF = "Within TRT between VISIT Differences", DELTA = "Between TRT between VISIT Difference") ``` ```{r eval=FALSE, include=FALSE} # Mini Mixed-Model example m <- lmer(DLQI_SCORE ~ VISIT*TRT + (1|USUBJID), data = df1) e <- emmeans(m, ~ VISIT*TRT) s <- eff_size(e, sigma = sigma(m), edf = 700, method = 'identity') c <- contrast(e, method = list(diff_a = c(-1, 1, 0, 0), diff_b = c( 0, 0, -1, 1), diff_d = c(1, -1, -1, 1))) f <- eff_size(c, sigma = sigma(m), edf = 700, method = 'identity') summary(e, infer = TRUE) summary(c, infer = TRUE) # Mini Repeated Mixed-Model example pacman::p_load(nlme) m <- gls(DLQI_SCORE ~ VISIT*TRT, data = df1, method = "REML", correlation = nlme::corSymm( form = ~VISIT_i|USUBJID), weights = nlme::varIdent(form = ~1|VISIT)) e <- emmeans(m, ~ VISIT*TRT, mode = "df.error") s <- eff_size(e, sigma = sigma(m), edf = 700, method = 'identity') c <- contrast(e, method = list(diff_a = c(-1, 1, 0, 0), diff_b = c( 0, 0, -1, 1), diff_d = c(1, -1, -1, 1))) f <- eff_size(c, sigma = sigma(m), edf = 700, method = 'identity') summary(e, infer = TRUE) summary(c, infer = TRUE) # Mini bayesian example pacman::p_load(brms) m <- brm(DLQI_SCORE ~ VISIT*TRT + (1 | USUBJID), data = df1, family = gaussian()) e <- emmeans(m, ~ VISIT*TRT, point.est = mean) s <- contrast(e, method = list(diff_a = c(-1, 1, 0, 0), diff_b = c( 0, 0, -1, 1), diff_d = c(1, -1, -1, 1))) summary(e, point.est = mean, infer = TRUE) summary(s, point.est = mean, infer = TRUE) ``` ```{r run ALL models} df2 <- df1 %>% pivot_longer(cols = contains("DLQI"), names_to = "VAR", values_to = "VAL") %>% nest_by(VAR) %>% mutate(m = list( lmer(VAL ~ VISIT*TRT + (1|USUBJID), data = data) ), e = list( emmeans(m, ~ VISIT*TRT) ), s = list( eff_size(e, sigma = sigma(m), edf = 700, method = 'identity') ), c = list( contrast(e, method = list(diff_a = c(-1, 1, 0, 0), diff_b = c( 0, 0, -1, 1), diff_d = c(1, -1, -1, 1))) ), f = list( eff_size(c, sigma = sigma(m), edf = 700, method = 'identity') ), r = list( bind_rows( summary(e, infer = TRUE) %>% as.data.frame(), summary(c, infer = TRUE) %>% as.data.frame() %>% rename(emmean = estimate, TRT = contrast) ) %>% mutate(contrast = paste(VISIT, TRT)) ), g = list( bind_rows( summary(s, infer = TRUE) %>% as.data.frame(), summary(f, infer = TRUE) %>% as.data.frame() ) ), ) ``` ```{r unnest raw} df3 <- df2 %>% select(VAR, r) %>% unnest(r) %>% mutate(contrast = str_remove(contrast, "NA ")) %>% mutate(GRP = case_when( contrast %in% c("Baseline A", "Week 16 A", "Baseline B", "Week 16 B") ~ "MEAN", contrast %in% c("diff_a", "diff_b") ~ "DIFF", contrast %in% c("diff_d") ~ "DELTA") ) %>% mutate(TRT = case_when( contrast %in% c("Baseline A", "Week 16 A", "diff_a") ~ "A", contrast %in% c("Baseline B", "Week 16 B", "diff_b") ~ "B") ) %>% mutate(VISIT = case_when( contrast %in% c("Baseline A", "Baseline B") ~ "Baseline", contrast %in% c("Week 16 A", "Week 16 B") ~ "Week 16") ) ``` ```{r unnest Standardized} df3s <- df2 %>% select(VAR, g) %>% unnest(g) %>% mutate(GRP = case_when( contrast %in% c("Baseline A", "Week 16 A", "Baseline B", "Week 16 B") ~ "MEAN", contrast %in% c("diff_a", "diff_b") ~ "DIFF", contrast %in% c("diff_d") ~ "DELTA") ) %>% mutate(TRT = case_when( contrast %in% c("Baseline A", "Week 16 A", "diff_a") ~ "A", contrast %in% c("Baseline B", "Week 16 B", "diff_b") ~ "B") ) %>% mutate(VISIT = case_when( contrast %in% c("Baseline A", "Baseline B") ~ "Baseline", contrast %in% c("Week 16 A", "Week 16 B") ~ "Week 16") ) ``` ```{r dendogram1} y1 <- df1 %>% #filter(VISIT == "Baseline") %>% select(starts_with("DLQI")) c1 <- cor(y1, method="pearson") d1 <- as.dist(1-c1) h1 <- hclust(d1, method = "complete") gg_d1 <- ggdendrogram(h1, rotate = TRUE) + labs(title = " ") + theme_dendro() + theme(plot.margin = margin(0, 0, 0, 0, "pt")) df3s <- df3s %>% mutate(VAR = factor(VAR, levels = rev(h1$labels[h1$order]) ), GRP = factor(GRP, levels = c("MEAN","DIFF","DELTA")) ) df3 <- df3 %>% mutate(VAR = factor(VAR, levels = rev(h1$labels[h1$order]) ), GRP = factor(GRP, levels = c("MEAN","DIFF","DELTA")) ) ``` ```{r table data} df4 <- df3 %>% select(VAR, C = contrast, E = emmean, L = lower.CL, U = upper.CL, P = p.value) %>% mutate(C = factor(C, labels = c("B_A","B_B","D_A","D_B","D_D","W_A","W_B")) %>% fct_relevel("B_A","B_B","W_A","W_B","D_A","D_B","D_D") ) %>% pivot_wider(id_cols = VAR, names_from = C, values_from = c(E, L, U, P)) df4s <- df3s %>% select(VAR, C = contrast, E = effect.size, L = lower.CL, U = upper.CL, P = p.value) %>% mutate(C = factor(C, labels = c("B_A","B_B","D_A","D_B","D_D","W_A","W_B")) %>% fct_relevel("B_A","B_B","W_A","W_B","D_A","D_B","D_D") ) %>% pivot_wider(id_cols = VAR, names_from = C, values_from = c(E, L, U, P)) df5 <- bind_rows(df4 %>% mutate(MODEL = "Raw"), df4s %>% mutate(MODEL = "Standardized")) ``` ```{r gg_f1} # Figure F1 gg_f1 <- ggplot(data = df3s, aes(x = effect.size, y = VAR)) + geom_vline(data = tribble(~GRP, ~VAL, "MEAN", NA, "DIFF", 0, "DELTA",0) %>% mutate(GRP = factor(GRP) %>% fct_relevel("MEAN","DIFF","DELTA")), aes(xintercept = VAL), color = 'gray50' ) + geom_vline(data = df3s %>% filter(GRP == "MEAN", VISIT == "Baseline") %>% select(GRP, VAR, effect.size) %>% group_by(GRP, VAR) %>% summarize(xint = mean(effect.size)), aes(xintercept = xint), color = 'gray50') + stat_dist_halfeye(ggsubset(TRT == "A"), mapping = aes(fill1 = contrast, color1 = contrast, dist = dist_student_t(df = df, mu = effect.size, sigma = SE)), .width = c(0.95), alpha = 0.5, justification = -0.10, size = 1, side = 'top') + stat_dist_halfeye(ggsubset(TRT == "B"), mapping = aes(fill2 = contrast, color2 = contrast, dist = dist_student_t(df = df, mu = effect.size, sigma = SE)), .width = c(0.95), alpha = 0.5, justification = 1.10, size = 1, side = 'bottom') + stat_dist_halfeye(ggsubset(GRP == "DELTA"), mapping = aes(dist = dist_student_t(df = df, mu = effect.size, sigma = SE)), .width = c(0.95), size = 1, side = 'both') + geom_richtext(ggsubset(GRP == "DELTA"), mapping = aes(y = 1, x = 3.7, label = str_glue("{style_ratio(effect.size)} ({style_ratio(lower.CL)}, {style_ratio(upper.CL)})
{style_pvalue(p.value, 2, prepend_p = TRUE)}")), size = 2.5, hjust = "inward", color = 'gray75', text.color = 'black', fill = 'white') + scale_listed(scalelist = list( scale_fill_manual(values = RColorBrewer::brewer.pal(9,"Blues")[c(5,7,8)], aesthetics = "fill1"), scale_fill_manual(values = RColorBrewer::brewer.pal(9,"Greens")[c(5,7,8)], aesthetics = "fill2"), scale_color_manual(values = RColorBrewer::brewer.pal(9,"Blues")[c(5,7,8)], aesthetics = "color1"), scale_color_manual(values = RColorBrewer::brewer.pal(9,"Greens")[c(5,7,8)], aesthetics = "color2")), replaces = c("fill", "fill", "color", "color")) + facet_grid(VAR ~ GRP, scales = 'free_y', switch = 'y', labeller = labeller(VAR = dl, GRP = cl)) + #force_panelsizes(cols = c(3,2,2)) + scale_y_discrete(position = "right") + labs(x = "Standardized effect-size (_Cohen's d_)", y = NULL) + theme_bw() + theme(plot.margin = margin(5.5, 0, 5.5, 5.5, "pt"), legend.position = "none", strip.text.y.left = element_markdown(angle = 0, hjust = 1), strip.placement = "outside", axis.ticks.y.right = element_blank(), axis.text.y.right = element_blank(), axis.title.x = element_markdown(hjust = 0.5) ) ``` #### {.tabset .tabset-pills} ##### Figure ```{r gg_f1 combined, fig.height = 9.5, fig.width = 12} gg_f1 + gg_d1 + plot_layout(widths = c(6,1)) + plot_annotation( title = "Phase 3 Randomized Clinical Trial Mixed-Model Results evaluating quality of life **DLQI Total Scale** & Individual Subscales", subtitle = "Standardized Effects (MEAN) and 95% CI for Placebo (n=150) and Active (n=450) at Baseline|Baseline and Week 16|Week 16 vertical lines are mean results at Baseline
Standardized Difference (DELTA) and 95% CI between VISIT within TRT (i.e. Active (Week 16 - Baseline) overlapping densities to the vertical lines not significant different from 0
Standardized Difference (DELTA-DELTA) and 95% CI between TRT between VISIT
Hierarchical Clustering dendrogram results of Pearson correlation between **DLQI Total** & Individual Subscales", theme = theme(plot.title.position = "plot", plot.title = element_markdown(lineheight = 1.1), plot.subtitle = element_markdown(lineheight = 1.1)) ) ``` ##### Table ```{r} # 1-ESTIMATE Estimate, Lower, Upper, Pvalue # 2-TIME Baseline, Week, Diff # 3-GROUP A, B, Diff dl2 <- data.frame(LAB = dl) %>% rownames_to_column(var = "VAR") %>% mutate(LAB = str_remove(LAB, "
")) df4s %>% select(-P_B_A, -P_B_B, -P_W_A, -P_W_B) %>% arrange(VAR) %>% ungroup() %>% left_join(dl2) %>% gt() %>% fmt_markdown(columns = TRUE) %>% fmt_number(columns = starts_with(c("E_","L_","U_")), decimals = 2) %>% fmt( columns = starts_with("P"), fns = function(x) { gtsummary::style_pvalue(x, digits = 2, prepend_p = TRUE) } ) %>% cols_merge( columns = vars(E_B_A, L_B_A, U_B_A), hide_columns = vars(L_B_A, U_B_A), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_W_A, L_W_A, U_W_A), hide_columns = vars(L_W_A, U_W_A), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_D_A, L_D_A, U_D_A, P_D_A), hide_columns = vars(L_D_A, U_D_A, P_D_A), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_A), rows = P_D_A < 0.05) ) %>% cols_merge( columns = vars(E_B_B, L_B_B, U_B_B), hide_columns = vars(L_B_B, U_B_B), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_W_B, L_W_B, U_W_B), hide_columns = vars(L_W_B, U_W_B), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_D_B, L_D_B, U_D_B, P_D_B), hide_columns = vars(L_D_B, U_D_B, P_D_B), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_B), rows = P_D_B < 0.05) ) %>% cols_merge( columns = vars(E_D_D, L_D_D, U_D_D, P_D_D), hide_columns = vars(L_D_D, U_D_D, P_D_D), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_D), rows = P_D_D < 0.05) ) %>% # colors tab_style( style = list( cell_fill(color = "#6BAED6", alpha = 0.5) ), locations = cells_body( columns = vars(E_B_A) ) ) %>% tab_style( style = list( cell_fill(color = "#2171B5", alpha = 0.5) ), locations = cells_body( columns = vars(E_W_A) ) ) %>% tab_style( style = list( cell_fill(color = "#08519C", alpha = 0.5) ), locations = cells_body( columns = vars(E_D_A) ) ) %>% tab_style( style = list( cell_fill(color = "#74C476", alpha = 0.5) ), locations = cells_body( columns = vars(E_B_B) ) ) %>% tab_style( style = list( cell_fill(color = "#238B45", alpha = 0.5) ), locations = cells_body( columns = vars(E_W_B) ) ) %>% tab_style( style = list( cell_fill(color = "#006D2C", alpha = 0.5) ), locations = cells_body( columns = vars(E_D_B) ) ) %>% tab_style( style = list( cell_fill(color = "#bfbfbf") ), locations = cells_body( columns = vars(LAB, E_D_D) ) ) %>% cols_move_to_start(columns = vars(LAB)) %>% cols_hide(columns = vars(VAR)) %>% cols_label( E_B_A = "Baseline", E_W_A = "Week 16", E_B_B = "Baseline", E_W_B = "Week 16", E_D_A = "Placebo", E_D_B = "Active", E_D_D = "Difference" ) %>% cols_width(starts_with(c("E_")) ~ px(135), vars(LAB) ~ px(200)) %>% cols_align("right", columns = vars(LAB)) %>% tab_spanner(label = "Placebo", vars(E_B_A, E_W_A)) %>% tab_spanner(label = "Active", vars(E_B_B, E_W_B)) %>% tab_spanner(label = "Within TRT between VISIT", vars(E_D_A, E_D_B)) %>% tab_spanner(label = "Overall", vars(E_D_D)) %>% tab_options(table.font.size = px(12), data_row.padding = px(5), table_body.hlines.width = px(5), table_body.hlines.color = "white") %>% opt_table_lines(extent = "default") %>% tab_header(title = html("Phase 3 Randomized Clinical Trial Mixed-Model Results evaluating quality of life DLQI Total Scale & Individual Subscales"), subtitle = "Standardized Effects (95% CI) and Standardized Differences (95% CI) and p-values") %>% tab_stubhead(label = "") %>% opt_align_table_header( align = "left") ``` ##### Table (Raw) ```{r} df4 %>% select(-P_B_A, -P_B_B, -P_W_A, -P_W_B) %>% arrange(df4s$VAR) %>% ungroup() %>% left_join(dl2) %>% gt() %>% fmt_markdown(columns = TRUE) %>% fmt_number(columns = starts_with(c("E_","L_","U_")), decimals = 2) %>% fmt( columns = starts_with("P"), fns = function(x) { gtsummary::style_pvalue(x, digits = 2, prepend_p = TRUE) } ) %>% cols_merge( columns = vars(E_B_A, L_B_A, U_B_A), hide_columns = vars(L_B_A, U_B_A), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_W_A, L_W_A, U_W_A), hide_columns = vars(L_W_A, U_W_A), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_D_A, L_D_A, U_D_A, P_D_A), hide_columns = vars(L_D_A, U_D_A, P_D_A), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_A), rows = P_D_A < 0.05) ) %>% cols_merge( columns = vars(E_B_B, L_B_B, U_B_B), hide_columns = vars(L_B_B, U_B_B), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_W_B, L_W_B, U_W_B), hide_columns = vars(L_W_B, U_W_B), pattern = "{1} ({2}, {3})" ) %>% cols_merge( columns = vars(E_D_B, L_D_B, U_D_B, P_D_B), hide_columns = vars(L_D_B, U_D_B, P_D_B), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_B), rows = P_D_B < 0.05) ) %>% cols_merge( columns = vars(E_D_D, L_D_D, U_D_D, P_D_D), hide_columns = vars(L_D_D, U_D_D, P_D_D), pattern = "{1} ({2}, {3})
{4}" ) %>% tab_style( style = list( cell_text(weight = "bold") ), locations = cells_body( columns = vars(E_D_D), rows = P_D_D < 0.05) ) %>% # colors tab_style( style = list( cell_fill(color = "#6BAED6", alpha = 0.5) ), locations = cells_body( columns = vars(E_B_A) ) ) %>% tab_style( style = list( cell_fill(color = "#2171B5", alpha = 0.5) ), locations = cells_body( columns = vars(E_W_A) ) ) %>% tab_style( style = list( cell_fill(color = "#08519C", alpha = 0.5) ), locations = cells_body( columns = vars(E_D_A) ) ) %>% tab_style( style = list( cell_fill(color = "#74C476", alpha = 0.5) ), locations = cells_body( columns = vars(E_B_B) ) ) %>% tab_style( style = list( cell_fill(color = "#238B45", alpha = 0.5) ), locations = cells_body( columns = vars(E_W_B) ) ) %>% tab_style( style = list( cell_fill(color = "#006D2C", alpha = 0.5) ), locations = cells_body( columns = vars(E_D_B) ) ) %>% tab_style( style = list( cell_fill(color = "#bfbfbf") ), locations = cells_body( columns = vars(LAB, E_D_D) ) ) %>% cols_move_to_start(columns = vars(LAB)) %>% cols_hide(columns = vars(VAR)) %>% cols_label( E_B_A = "Baseline", E_W_A = "Week 16", E_B_B = "Baseline", E_W_B = "Week 16", E_D_A = "Placebo", E_D_B = "Active", E_D_D = "Difference" ) %>% cols_width(starts_with(c("E_")) ~ px(135), vars(LAB) ~ px(200)) %>% cols_align("right", columns = vars(LAB)) %>% tab_spanner(label = "Placebo", vars(E_B_A, E_W_A)) %>% tab_spanner(label = "Active", vars(E_B_B, E_W_B)) %>% tab_spanner(label = "Within TRT between VISIT", vars(E_D_A, E_D_B)) %>% tab_spanner(label = "Overall", vars(E_D_D)) %>% tab_options(table.font.size = px(12), data_row.padding = px(5), table_body.hlines.width = px(5), table_body.hlines.color = "white") %>% tab_header(title = html("Phase 3 Randomized Clinical Trial Mixed-Model Results evaluating quality of life DLQI Total Scale & Individual Subscales"), subtitle = "Effects (95% CI) and Differences (95% CI) and p-values") %>% tab_stubhead(label = "") %>% opt_align_table_header( align = "left") ```