--- title: "Combination of inflammatory and vascular markers in the febrile phase of dengue is associated with more severe outcomes" author: "Nguyen Lam Vuong" date: "29-Jul-2021" output: github_document --- ```{r load packages, echo=TRUE, message=FALSE, warning=FALSE} library(tidyverse) library(gtsummary) # need to install package 'flextable' library(rms) library(MuMIn) # for best subset selection library(facetscales) # need to install package 'facetscales' from devtools::install_github("zeehio/facetscales") source("Elife ERA functions.R") # to include my functions options(gtsummary.tbl_summary.percent_fun = function(x) sprintf(x * 100, fmt='%1.0f')) # to report percentages without decimal options(knitr.kable.NA = '') # to set NA to '' in kable results theme_set(theme_bw()) # I love black & white theme options(na.action = "na.fail") # for 'dredge' function [MuMIn] ``` ```{r load data, echo=TRUE, message=FALSE, warning=FALSE} # Full data dat0 <- read_csv("Dengue_Biomarkers_data_27Jul2021.csv") %>% mutate(group2 = factor(sev.or.inte, levels=c(0,1), labels=c("Uncomplicated dengue","Severe/moderate dengue")), Country = as.factor(Country), Serotype = as.factor(Serotype), Serology = factor(Serology, levels = c("Probable primary", "Probable secondary", "Inconclusive")), WHO2009 = factor(WHO2009, levels = c("Mild dengue", "Dengue with warning signs", "Severe dengue", "Unknown"))) # Data at enrollment (for Table 1 & Appendix 4-table 1) dat <- dat0 %>% filter(Time == "Enrolment") # Data at enrollment for models ## with inverse probability weights (IPW) for inclusion probability by countries (for the analysis of secondary outcomes) ## transform the biomarker's levels to log-2 and viremia to log-10 ## set all biomarker's levels under the limit of detection (u...=1) to the limit of detection dat1 <- dat %>% mutate(age15 = factor(age15, levels=c("No","Yes"), labels=c("Under 15","15 and above")), Serotype = ifelse(Serotype=="Unknown", NA, as.character(Serotype)), Serotype = ifelse(Serotype=="DENV-1", 1, 2), Serotype2 = Serotype, # for getting results for Appendix5-tables 1, 2 Serotype = factor(Serotype, levels=c(1,2), labels=c("DENV-1","Others")), # set Serotype to DENV-1 and others ipw = ifelse(sev.or.inte==1, 1, ifelse(Country=="Vietnam", (1505-204)/436, ifelse(Country=="Malaysia", (259-29)/58, ifelse(Country=="El Salvador", (306-18)/23, ifelse(Country=="Cambodia", (302-30)/39, NA))))), ipwsd = ipw * 837/2372, VCAM = ifelse(uVCAM==1, log2(0.028), log2(VCAM1)), SDC = log2(SDC1), Ang = ifelse(uAng==1, log2(17.1), log2(Ang2)), IL8 = ifelse(uIL8==1, log2(1.8), log2(IL8)), IP10 = ifelse(uIP10==1, log2(1.18), log2(IP10)), IL1RA = ifelse(uIL1RA==1, log2(18), log2(IL1RA)), CD163 = log2(CD163), TREM = ifelse(uTREM==1, log2(10.65), log2(TREM1)), Fer = log2(Fer), CRP = log2(CRP), Vir = log10(Viremia)) %>% select(Code, Age, age15, Serotype, Serotype2, sev.or.inte, sev.only, sev.or.ws, hospital, VCAM, SDC, Ang, IL8, IP10, IL1RA, CD163, TREM, Fer, CRP, Vir, uVCAM, uAng, uIL8, uIP10, uCD163, uTREM, uVir, ipw, ipwsd) dat1c <- dat1 %>% filter(age15 == "Under 15") ## Data for children dat1a <- dat1 %>% filter(age15 == "15 and above") ## Data for adults # Set references for calculating the ORs and 95% CIs ref0 <- c(median(dat1$VCAM), median(dat1$SDC), median(dat1$Ang), median(dat1$IL8), median(dat1$IP10), median(dat1$IL1RA), median(dat1$CD163), median(dat1$TREM), median(dat1$Fer), median(dat1$CRP)) # Create data for plotting biomarkers (for Figure 2 & Appendix 4-figure 1) ## Enrollment tmp1 <- dat0 %>% filter(Time == "Enrolment") %>% select(Code, group2, Day, Daygr, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP) %>% gather(., "Biomarker", "Result", 5:14) ## Follow up tmp2 <- dat0 %>% filter(Time == "Follow up") %>% select(Code, group2, Day, Daygr, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP) %>% gather(., "Biomarker", "Result", 5:14) ## Merge long data for biomarkers dat_plot <- bind_rows(tmp1, tmp2) %>% mutate(Daygr = factor(Daygr, levels = c("Day 1", "Day 2", "Day 3", "Day 10-20", "Day >20"), labels = c("1", "2", "3", "10-20", ">20")), Biomarker = factor(Biomarker, levels=c("VCAM1", "SDC1", "Ang2", "IL8", "IP10", "IL1RA", "CD163", "TREM1", "Fer", "CRP"), labels=c("VCAM-1 (ng/ml)", "SDC-1 (pg/ml)", "Ang-2 (pg/ml)", "IL-8 (pg/ml)", "IP-10 (pg/ml)", "IL-1RA (pg/ml)", "sCD163 (ng/ml)", "sTREM-1 (pg/ml)", "Ferritin (ng/ml)", "CRP (mg/l)"))) ``` ########################################################################################### #### Table 1. Summary of clinical data by primary outcome. ```{r table 1, echo=TRUE, message=FALSE, warning=FALSE} # All patients t1 <- dat %>% select(group2, Country, Age, Sex, Day, Serotype, Serology, Obesity, Diabetes, WHO2009, hospital) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})", all_categorical() ~ "{n} ({p})"), value = list(Sex ~ "Male"), digits = list(all_continuous() ~ c(0,0)), label = list(Age ~ "Age (years)", Sex ~ "Gender male", Day ~ "Illness day at enrolment", Serology ~ "Immune status", WHO2009 ~ "WHO 2009 classification", hospital ~ "Hospitalization")) %>% add_stat_label(label = c(all_categorical() ~ "n (%)", Age ~ "median (1st, 3rd quartiles)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") # Children (<15 years of age) t2 <- dat %>% filter(age15=="No") %>% select(group2, Country, Age, Sex, Day, Serotype, Serology, Obesity, Diabetes, WHO2009, hospital) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})", all_categorical() ~ "{n} ({p})"), value = list(Sex ~ "Male"), digits = list(all_continuous() ~ c(0,0)), label = list(Age ~ "Age (years)", Sex ~ "Gender male", Day ~ "Illness day at enrolment", Serology ~ "Immune status", WHO2009 ~ "WHO 2009 classification", hospital ~ "Hospitalization")) %>% add_stat_label(label = c(all_categorical() ~ "n (%)", Age ~ "median (1st, 3rd quartiles)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") # Adults (>=15 years of age) t3 <- dat %>% filter(age15=="Yes") %>% select(group2, Country, Age, Sex, Day, Serotype, Serology, Obesity, Diabetes, WHO2009, hospital) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})", all_categorical() ~ "{n} ({p})"), value = list(Sex ~ "Male"), digits = list(all_continuous() ~ c(0,0)), label = list(Age ~ "Age (years)", Sex ~ "Gender male", Day ~ "Illness day at enrolment", Serology ~ "Immune status", WHO2009 ~ "WHO 2009 classification", hospital ~ "Hospitalization")) %>% add_stat_label(label = c(all_categorical() ~ "n (%)", Age ~ "median (1st, 3rd quartiles)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") # Merge tables tbl_merge(tbls = list(t1, t2, t3), tab_spanner = c("**All patients**", "**Children**", "**Adults**")) ``` ########################################################################################### #### Figure 2. Biomarker levels by groups. ```{r fig 2, echo=TRUE, message=FALSE, warning=FALSE, fig.height=8, fig.width=11} tick <- c(0,1,4,10,40,100,200,400,1000,4000,10000,20000,40000,70000) # for y-axis tick labels p2 <- dat_plot %>% ggplot(., aes(Daygr, Result^(1/4), fill=group2, color=group2)) + geom_boxplot(alpha=.5, outlier.size=.9, lwd=.4, fatten=1) + geom_boxplot(alpha=0, outlier.color=NA, color="black", lwd=.4, fatten=1) + facet_wrap(~ Biomarker, scales="free", ncol=5) + scale_y_continuous(breaks=tick^(1/4), labels=tick) + xlab("Illness day (day 1 [n=140]; day 2 [n=390]; day 3 [n=307]; day 10-20 [n=625]; day >20 [n=43])") + theme(axis.title.y=element_blank(), legend.position="top", legend.title=element_blank(), axis.text.y=element_text(size=rel(.8))) p2 #ggsave(filename="Fig 2.pdf", dpi=300, plot=p2, width=11, height=8) ``` ########################################################################################### #### Table 2. Results from models for the primary endpoint (severe or moderate dengue). ```{r table 2, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Use my function 'get_est1' to get ORs and CIs from models tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # Children - single models or1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Children - global model or2c = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo2c = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up2c = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Adults - single models or1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1), # Adults - global model or2a = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo2a = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up2a = get_est2(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Use my function 'get_est1' to get p-values from models tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s1 = get_est1(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), # P overall from single models p.s1_int = get_est1(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1), # P interaction from single models p.g1 = get_est2(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), # P overall from global models p.g1_int = get_est2(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1))) %>% # P interaction from global models ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results into a table res1 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(!is.na(ref1), paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), as.character(bio)), or.sc1 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), # s: single model; c: children or.gc1 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), # g: global model or.sa1 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), # a: adults or.ga1 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc1, or.sa1, p.s1, p.s1_int, or.gc1, or.ga1, p.g1, p.g1_int) names(res1) <- c("", "OR (children - single)", "OR (adults - single)", "P overall (single)", "P interaction (single)", "OR (children - global)", "OR (adults - global)", "P overall (global)", "P interaction (global)") # Report the results knitr::kable(res1) ``` ########################################################################################### #### Figure 3. Results from models for the primary endpoint (severe or moderate dengue). ```{r fig 3, echo=TRUE, message=FALSE, warning=FALSE, fig.height=8, fig.width=11} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get results from models for children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1(out="sev.or.inte", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m2 <- get_pred1(out="sev.or.inte", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m3 <- get_pred1(out="sev.or.inte", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m4 <- get_pred1(out="sev.or.inte", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m5 <- get_pred1(out="sev.or.inte", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m6 <- get_pred1(out="sev.or.inte", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m7 <- get_pred1(out="sev.or.inte", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m8 <- get_pred1(out="sev.or.inte", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m9 <- get_pred1(out="sev.or.inte", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m10 <- get_pred1(out="sev.or.inte", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg1 <- get_pred2(out="sev.or.inte", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg2 <- get_pred2(out="sev.or.inte", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg3 <- get_pred2(out="sev.or.inte", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg4 <- get_pred2(out="sev.or.inte", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg5 <- get_pred2(out="sev.or.inte", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg6 <- get_pred2(out="sev.or.inte", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg7 <- get_pred2(out="sev.or.inte", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg8 <- get_pred2(out="sev.or.inte", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg9 <- get_pred2(out="sev.or.inte", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg10 <- get_pred2(out="sev.or.inte", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "10 years") # Get results from models for adults dd$limits["Adjust to","Age"] <- 25 dat_m1 <- get_pred1(out="sev.or.inte", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m2 <- get_pred1(out="sev.or.inte", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m3 <- get_pred1(out="sev.or.inte", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m4 <- get_pred1(out="sev.or.inte", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m5 <- get_pred1(out="sev.or.inte", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m6 <- get_pred1(out="sev.or.inte", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m7 <- get_pred1(out="sev.or.inte", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m8 <- get_pred1(out="sev.or.inte", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m9 <- get_pred1(out="sev.or.inte", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m10 <- get_pred1(out="sev.or.inte", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg1 <- get_pred2(out="sev.or.inte", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg2 <- get_pred2(out="sev.or.inte", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg3 <- get_pred2(out="sev.or.inte", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg4 <- get_pred2(out="sev.or.inte", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg5 <- get_pred2(out="sev.or.inte", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg6 <- get_pred2(out="sev.or.inte", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg7 <- get_pred2(out="sev.or.inte", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg8 <- get_pred2(out="sev.or.inte", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg9 <- get_pred2(out="sev.or.inte", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg10 <- get_pred2(out="sev.or.inte", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "25 years") # Merge results for children and adults vis1 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 1") # Merge data for plots tmp0 <- dat1 %>% arrange(age15) %>% select(sev.or.inte) tmp <- data.frame(sev.or.inte = rep(tmp0$sev.or.inte, 20)) tmp_p1 <- vis1 %>% arrange(biomarker, model) %>% mutate(value1 = 2^value, age = factor(age, levels=c("10 years", "25 years")), gr = ifelse(model=="Single model" & age=="10 years", 1, ifelse(model=="Single model" & age=="25 years", 2, ifelse(model=="Global model" & age=="10 years", 3, 4))), gr = factor(gr, levels=c(1:4), labels=c("Single children", "Single adults", "Global children", "Global adults")), model = factor(model, levels=c("Single model", "Global model"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(model=="Single model") %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #require(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.125, .25, .5, 1, 2, 4, 8) ybreak2 <- c(.125, .25, .5, 1, 2, 4) scales_y1 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1), `Global model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2), `Global model` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names models <- c(`Single model` = "Single model", `Global model` = "Global model") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.inte==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.inte==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.inte==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.inte==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(models, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.or.inte==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.or.inte==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.or.inte==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.or.inte==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(models, biomarkers))) # Merge plots p3 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.86)) #ggsave(filename="Fig 3.pdf", p3, dpi=300, width=11, height=8) ``` ########################################################################################### #### Table 3. Best combinations of biomarkers associated with severe or moderate dengue for children. ```{r table 3, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + IP10 + IL1RA + CD163 + TREM + Fer + CRP, family=binomial, data=dat1c, x=T, y=T) # Selected model --------------------------------------------------- sel_var <- matrix(0, ncol=length(pred)+1, nrow=5, dimnames=list(NULL, c(pred, "aic"))) for (i in 1:5) { if (i==1) { bs <- dredge(full_mod, rank="AIC") } else { bs <- dredge(full_mod, rank="AIC", m.lim=c(i,i)) } bs_var <- attr(get.models(bs, 1)[[1]]$terms, "term.labels") for (j in 1:(ncol(sel_var)-1)) {sel_var[i,j] <- ifelse(names(sel_var[i,j]) %in% bs_var, 1, 0)} formula <- paste("sev.or.inte~", paste(names(sel_var[i,][sel_var[i,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1c, family = binomial, x = T, y = T) sel_var[i, ncol(sel_var)] <- AIC(sel_mod) } # Report results -------------------------------------------------- out1 <- as.data.frame(sel_var) %>% mutate(AIC = round(aic,1)) %>% select(-aic) for (i in 1:(ncol(out1)-1)) {out1[,i] <- ifelse(out1[,i]==0, NA, "+")} out2 <- as.data.frame(t(out1)) colnames(out2) <- c("Best of all combinations", "Best combination of 2 variables", "Best combination of 3 variables", "Best combination of 4 variables", "Best combination of 5 variables") rownames(out2) <- c("- VCAM-1", "- SDC-1", "- Ang-2", "- IL-8", "- IP-10", "- IL-1RA", "- sCD163", "- sTREM-1", "- Ferritin", "- CRP", "AIC of the selected model") knitr::kable(out2) # For bootstrap results please look at Appendix 7-tables 1, 2 (for the best of all combinations) # and the bootstrap codes for the best combinations of 2, 3, 4, and 5 variables ``` ########################################################################################### #### Table 4. Best combinations of biomarkers associated with severe or moderate dengue for adults. ```{r table 4, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "ns1(IP10)", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + ns1(IP10) + IL1RA + CD163 + TREM + Fer + CRP, family=binomial, data=dat1a, x=T, y=T) # Selected model --------------------------------------------------- sel_var <- matrix(0, ncol=length(pred)+1, nrow=5, dimnames=list(NULL, c(pred, "aic"))) for (i in 1:5) { if (i==1) { bs <- dredge(full_mod, rank="AIC") } else { bs <- dredge(full_mod, rank="AIC", m.lim=c(i,i)) } bs_var <- attr(get.models(bs, 1)[[1]]$terms, "term.labels") for (j in 1:(ncol(sel_var)-1)) {sel_var[i,j] <- ifelse(names(sel_var[i,j]) %in% bs_var, 1, 0)} formula <- paste("sev.or.inte~", paste(names(sel_var[i,][sel_var[i,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1a, family = binomial, x = T, y = T) sel_var[i, ncol(sel_var)] <- AIC(sel_mod) } # Report results -------------------------------------------------- out1 <- as.data.frame(sel_var) %>% mutate(AIC = round(aic,1)) %>% select(-aic) for (i in 1:(ncol(out1)-1)) {out1[,i] <- ifelse(out1[,i]==0, NA, "+")} out2 <- as.data.frame(t(out1)) colnames(out2) <- c("Best of all combinations", "Best combination of 2 variables", "Best combination of 3 variables", "Best combination of 4 variables", "Best combination of 5 variables") rownames(out2) <- c("- VCAM-1", "- SDC-1", "- Ang-2", "- IL-8", "- IP-10", "- IL-1RA", "- sCD163", "- sTREM-1", "- Ferritin", "- CRP", "AIC of the selected model") knitr::kable(out2) # For bootstrap results please look at Appendix 7-tables 3, 4 (for the best of all combinations) # and the bootstrap codes for the best combinations of 2, 3, 4, and 5 variables ``` ########################################################################################### #### Appendix 3. Statistical analysis. Analysis to find the best combination of biomarkers to predict the primary endpoint (aim #2) - Step #1 - AICs of models ```{r Appendix 3, echo=TRUE, message=FALSE, warning=FALSE} # For children screen_child <- data.frame(bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% group_by(bio) %>% do(cbind(., model1_child = screen(mod=1, bio=.$bio, dat=dat1c), model2_child = screen(mod=2, bio=.$bio, dat=dat1c), model3_child = screen(mod=3, bio=.$bio, dat=dat1c), model4_child = screen(mod=4, bio=.$bio, dat=dat1c))) %>% ungroup() %>% mutate(bio = factor(bio, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"))) %>% arrange(bio) # For adults screen_adult <- data.frame(bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% group_by(bio) %>% do(cbind(., model1_adult = screen(mod=1, bio=.$bio, dat=dat1a), model2_adult = screen(mod=2, bio=.$bio, dat=dat1a), model3_adult = screen(mod=3, bio=.$bio, dat=dat1a), model4_adult = screen(mod=4, bio=.$bio, dat=dat1a))) %>% ungroup() %>% mutate(bio = factor(bio, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"))) %>% arrange(bio) # Combine results out <- full_join(screen_child, screen_adult, by="bio") knitr::kable(out) ``` ########################################################################################### #### Appendix 4—table 1. Summary of clinical phenotype of the primary endpoint. ```{r table A4.1, echo=TRUE, message=FALSE, warning=FALSE} dat %>% filter(sev.or.inte == 1) %>% mutate(age15 = factor(age15, levels = c("No", "Yes"), labels = c("Children", "Adults")), mod.only = sev.or.inte - sev.only) %>% # to create 'moderate dengue' select(age15, sev.only, sev.leak, shock, res.dis, sev.neu, sev.bleed, sev.other, sev.liver, mod.only, mod.leak, mod.liver, mod.bleed, mod.other, mod.neu) %>% tbl_summary(by = age15, statistic = list(all_categorical() ~ "{n} ({p})"), label = list(sev.only ~ "Severe dengue, n (%)", sev.leak ~ "- Severe plasma leakage", shock ~ " + Dengue shock syndrome", res.dis ~ " + Respiratory distress", sev.neu ~ "- Severe neurologic involvement", sev.bleed ~ "- Severe bleeding", sev.other ~ "- Severe other major organ failure", sev.liver ~ "- Severe hepatic involvement", mod.only ~ "Moderate dengue, n (%)", mod.leak ~ "- Moderate plasma leakage", mod.liver ~ "- Moderate hepatic involvement", mod.bleed ~ "- Moderate bleeding", mod.other ~ "- Moderate other major organ involvement", mod.neu ~ "- Moderate neurologic involvement")) %>% add_overall() %>% modify_header(label = "", stat_0 = "**All patients (N={N})**", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) ``` ########################################################################################### #### Appendix 4—table 2. Summary of biomarkers’ data (at enrollment) ```{r table A4.2a, echo=TRUE, message=FALSE, warning=FALSE} # All patients t1 <- dat0 %>% filter(Time == "Enrolment") %>% mutate(viremia = Viremia/(10^6)) %>% # to summarize Viremia as 10^6 copies/ml select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP, viremia) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), viremia ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)", viremia ~ "Viremia (10^6 copies/ml)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Children (<15 years of age) t2 <- dat0 %>% filter(Time == "Enrolment") %>% filter(age15 == "No") %>% mutate(viremia = Viremia/(10^6)) %>% # to summarize Viremia as 10^6 copies/ml select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP, viremia) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), viremia ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)", viremia ~ "Viremia (10^6 copies/ml)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Adults (>=15 years of age) t3 <- dat0 %>% filter(Time == "Enrolment") %>% filter(age15 == "Yes") %>% mutate(viremia = Viremia/(10^6)) %>% # to summarize Viremia as 10^6 copies/ml select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP, viremia) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), viremia ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)", viremia ~ "Viremia (10^6 copies/ml)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Merge tables tbl_merge(tbls = list(t1, t2, t3), tab_spanner = c("**All patients**", "**Children**", "**Adults**")) ``` #### Appendix 4—table 2. Summary of biomarkers’ data (at follow-up) ```{r table A4.2b, echo=TRUE, message=FALSE, warning=FALSE} # All patients t1 <- dat0 %>% filter(Time == "Follow up") %>% select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), c(IL8, CRP) ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Children (<15 years of age) t2 <- dat0 %>% filter(Time == "Follow up") %>% filter(age15 == "No") %>% select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), c(IL8, CRP) ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Adults (>=15 years of age) t3 <- dat0 %>% filter(Time == "Follow up") %>% filter(age15 == "Yes") %>% select(group2, VCAM1, SDC1, Ang2, IL8, IP10, IL1RA, CD163, TREM1, Fer, CRP) %>% tbl_summary(by = group2, statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), digits = list(all_continuous() ~ c(0,0), c(IL8, CRP) ~ c(1,1)), label = list(VCAM1 ~ "VCAM-1 (ng/ml)", SDC1 ~ "SDC-1 (pg/ml)", Ang2 ~ "Ang-2 (pg/ml)", IL8 ~ "IL-8 (pg/ml)", IP10 ~ "IP-10 (pg/ml)", IL1RA ~ "IL-1RA (pg/ml)", CD163 ~ "sCD163 (ng/ml)", TREM1 ~ "sTREM-1 (pg/ml)", Fer ~ "Ferritin (ng/ml)", CRP ~ "CRP (mg/l)")) %>% modify_header(label = "", stat_by = "**{level} (N={n})**") %>% modify_footnote(update = everything() ~ NA) # Merge tables tbl_merge(tbls = list(t1, t2, t3), tab_spanner = c("**All patients**", "**Children**", "**Adults**")) ``` ########################################################################################### #### Appendix 4-figure 1. Biomarker levels by individual. ```{r fig A4.1, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} tick <- c(0,1,10,40,100,200,400,1000,4000,10000,20000,40000,70000) # for y-axis tick labels p41 <- ggplot(dat_plot, aes(Day, Result^(1/4), color=group2)) + geom_point(alpha = 0.5, size = .5) + geom_line(aes(group = Code), alpha = 0.07) + facet_wrap(~ Biomarker, scales="free", ncol=5) + scale_y_continuous(breaks=tick^(1/4), labels=tick) + scale_x_continuous(breaks = c(1,3,10,15,20,25,30), limit = c(1,31), name = "Illness day") + theme_bw() + theme(axis.title.y=element_blank(), legend.position="top", legend.title=element_blank(), axis.text.y=element_text(size=rel(.8))) p41 #ggsave(filename="A4_fig1.pdf", plot=p41, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 4—figure 2. Pairwise correlation of biomarker levels at enrollment and age. ```{r fig A4.2, echo=TRUE, message=FALSE, warning=FALSE, fig.height=9.5, fig.width=9.5} collabel <- c("Age", "VCAM-1", "SDC-1", "Ang-2", "IL-8", "IP-10", "IL-1RA", "sCD163", "sTREM-1", "Ferritin", "CRP", "Viremia") # Change biomarker's values to log-2 and viremia to log dat_fig42 <- dat %>% mutate(log2_VCAM1 = as.numeric(log2(VCAM1)), log2_SDC1 = as.numeric(log2(SDC1)), log2_Ang2 = as.numeric(log2(Ang2)), log2_IL8 = as.numeric(log2(IL8)), log2_IP10 = as.numeric(log2(IP10)), log2_IL1RA = as.numeric(log2(IL1RA)), log2_CD163 = as.numeric(log2(CD163)), log2_TREM1 = as.numeric(log2(TREM1)), log2_Fer = as.numeric(log2(Fer)), log2_CRP = as.numeric(log2(CRP)), log_vir = as.numeric(log10(Viremia))) %>% select(Age, log2_VCAM1, log2_SDC1, log2_Ang2, log2_IL8, log2_IP10, log2_IL1RA, log2_CD163, log2_TREM1, log2_Fer, log2_CRP, log_vir) %>% as.data.frame(.) p42 <- GGally::ggpairs(dat_fig42, lower = list(continuous = GGscatterPlot), # GGscatterPlot is a customized function upper = "blank", columnLabels = collabel) + theme(panel.grid.minor = element_blank()) p42 #ggsave(filename="A4_fig2.pdf", plot=p42, width=9.5, height=9.5) ``` ########################################################################################### #### Appendix 5—figure 1. Results from single models for severe/moderate dengue with the interaction with serotype. ```{r fig A5.1, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For DENV-1 dd$limits["Adjust to","Serotype"] <- "DENV-1" ## For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1s(out="sev.or.inte", bio="VCAM", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m2 <- get_pred1s(out="sev.or.inte", bio="SDC", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m3 <- get_pred1s(out="sev.or.inte", bio="Ang", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m4 <- get_pred1s(out="sev.or.inte", bio="IL8", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m5 <- get_pred1s(out="sev.or.inte", bio="IP10", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m6 <- get_pred1s(out="sev.or.inte", bio="IL1RA", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m7 <- get_pred1s(out="sev.or.inte", bio="CD163", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m8 <- get_pred1s(out="sev.or.inte", bio="TREM", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m9 <- get_pred1s(out="sev.or.inte", bio="Fer", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m10 <- get_pred1s(out="sev.or.inte", bio="CRP", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") ## For adults dd$limits["Adjust to","Age"] <- 25 dat_mg1 <- get_pred1s(out="sev.or.inte", bio="VCAM", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg2 <- get_pred1s(out="sev.or.inte", bio="SDC", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg3 <- get_pred1s(out="sev.or.inte", bio="Ang", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg4 <- get_pred1s(out="sev.or.inte", bio="IL8", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg5 <- get_pred1s(out="sev.or.inte", bio="IP10", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg6 <- get_pred1s(out="sev.or.inte", bio="IL1RA", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg7 <- get_pred1s(out="sev.or.inte", bio="CD163", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg8 <- get_pred1s(out="sev.or.inte", bio="TREM", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg9 <- get_pred1s(out="sev.or.inte", bio="Fer", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg10 <- get_pred1s(out="sev.or.inte", bio="CRP", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(serotype = "DENV-1") # For others dd$limits["Adjust to","Serotype"] <- "Others" ## For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1s(out="sev.or.inte", bio="VCAM", age=10, serotype="Others", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m2 <- get_pred1s(out="sev.or.inte", bio="SDC", age=10, serotype="Others", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m3 <- get_pred1s(out="sev.or.inte", bio="Ang", age=10, serotype="Others", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m4 <- get_pred1s(out="sev.or.inte", bio="IL8", age=10, serotype="Others", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m5 <- get_pred1s(out="sev.or.inte", bio="IP10", age=10, serotype="Others", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m6 <- get_pred1s(out="sev.or.inte", bio="IL1RA", age=10, serotype="Others", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m7 <- get_pred1s(out="sev.or.inte", bio="CD163", age=10, serotype="Others", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m8 <- get_pred1s(out="sev.or.inte", bio="TREM", age=10, serotype="Others", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m9 <- get_pred1s(out="sev.or.inte", bio="Fer", age=10, serotype="Others", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m10 <- get_pred1s(out="sev.or.inte", bio="CRP", age=10, serotype="Others", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") ## For adults dd$limits["Adjust to","Age"] <- 25 dat_mg1 <- get_pred1s(out="sev.or.inte", bio="VCAM", age=25, serotype="Others", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg2 <- get_pred1s(out="sev.or.inte", bio="SDC", age=25, serotype="Others", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg3 <- get_pred1s(out="sev.or.inte", bio="Ang", age=25, serotype="Others", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg4 <- get_pred1s(out="sev.or.inte", bio="IL8", age=25, serotype="Others", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg5 <- get_pred1s(out="sev.or.inte", bio="IP10", age=25, serotype="Others", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg6 <- get_pred1s(out="sev.or.inte", bio="IL1RA", age=25, serotype="Others", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg7 <- get_pred1s(out="sev.or.inte", bio="CD163", age=25, serotype="Others", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg8 <- get_pred1s(out="sev.or.inte", bio="TREM", age=25, serotype="Others", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg9 <- get_pred1s(out="sev.or.inte", bio="Fer", age=25, serotype="Others", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg10 <- get_pred1s(out="sev.or.inte", bio="CRP", age=25, serotype="Others", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(serotype = "Others") # Merge data for children and adults vis1 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 1") # Merge data for plots tmp0 <- dat1 %>% filter(!is.na(Serotype)) %>% arrange(Serotype) %>% select(sev.or.inte) tmp <- data.frame(sev.or.inte = rep(tmp0$sev.or.inte, 20)) tmp_p1 <- vis1 %>% arrange(biomarker, age) %>% mutate(value1 = 2^value, serotype = factor(serotype, levels=c("DENV-1", "Others")), age = factor(age, levels=c("Children", "Adults"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #require(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.125, .25, .5, 1, 2, 4, 8) ybreak2 <- c(.125, .25, .5, 1, 2, 4) scales_y1 <- list( `Children` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1), `Adults` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Children` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2), `Adults` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names age <- c(`Children` = "Children", `Adults` = "Adults") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=serotype), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=serotype), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, age=="Children" & sev.or.inte==0 & serotype=="DENV-1"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, age=="Children" & sev.or.inte==1 & serotype=="DENV-1"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, age=="Adults" & sev.or.inte==0 & serotype=="Others"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, age=="Adults" & sev.or.inte==1 & serotype=="Others"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + scale_fill_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(age), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(age, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=serotype), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=serotype), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, age=="Children" & sev.or.inte==0 & serotype=="DENV-1"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, age=="Children" & sev.or.inte==1 & serotype=="DENV-1"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, age=="Adults" & sev.or.inte==0 & serotype=="Others"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, age=="Adults" & sev.or.inte==1 & serotype=="Others"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + scale_fill_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(age), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(age, biomarkers))) # Merge plots p51 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.86)) #ggsave(filename="A5_fig1.pdf", p3, dpi=300, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 5—figure 2. Results from global model for severe/moderate dengue with the interaction with serotype. ```{r fig A5.2, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For DENV-1 dd$limits["Adjust to","Serotype"] <- "DENV-1" ## For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred2s(out="sev.or.inte", bio="VCAM", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m2 <- get_pred2s(out="sev.or.inte", bio="SDC", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m3 <- get_pred2s(out="sev.or.inte", bio="Ang", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m4 <- get_pred2s(out="sev.or.inte", bio="IL8", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m5 <- get_pred2s(out="sev.or.inte", bio="IP10", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m6 <- get_pred2s(out="sev.or.inte", bio="IL1RA", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m7 <- get_pred2s(out="sev.or.inte", bio="CD163", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m8 <- get_pred2s(out="sev.or.inte", bio="TREM", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m9 <- get_pred2s(out="sev.or.inte", bio="Fer", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") dat_m10 <- get_pred2s(out="sev.or.inte", bio="CRP", age=10, serotype="DENV-1", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Children") ## For adults dd$limits["Adjust to","Age"] <- 25 dat_mg1 <- get_pred2s(out="sev.or.inte", bio="VCAM", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg2 <- get_pred2s(out="sev.or.inte", bio="SDC", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg3 <- get_pred2s(out="sev.or.inte", bio="Ang", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg4 <- get_pred2s(out="sev.or.inte", bio="IL8", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg5 <- get_pred2s(out="sev.or.inte", bio="IP10", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg6 <- get_pred2s(out="sev.or.inte", bio="IL1RA", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg7 <- get_pred2s(out="sev.or.inte", bio="CD163", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg8 <- get_pred2s(out="sev.or.inte", bio="TREM", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg9 <- get_pred2s(out="sev.or.inte", bio="Fer", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_mg10 <- get_pred2s(out="sev.or.inte", bio="CRP", age=25, serotype="DENV-1", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='DENV-1')) %>% mutate(age="Adults") dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(serotype = "DENV-1") # For others dd$limits["Adjust to","Serotype"] <- "Others" ## For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred2s(out="sev.or.inte", bio="VCAM", age=10, serotype="Others", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m2 <- get_pred2s(out="sev.or.inte", bio="SDC", age=10, serotype="Others", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m3 <- get_pred2s(out="sev.or.inte", bio="Ang", age=10, serotype="Others", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m4 <- get_pred2s(out="sev.or.inte", bio="IL8", age=10, serotype="Others", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m5 <- get_pred2s(out="sev.or.inte", bio="IP10", age=10, serotype="Others", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m6 <- get_pred2s(out="sev.or.inte", bio="IL1RA", age=10, serotype="Others", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m7 <- get_pred2s(out="sev.or.inte", bio="CD163", age=10, serotype="Others", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m8 <- get_pred2s(out="sev.or.inte", bio="TREM", age=10, serotype="Others", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m9 <- get_pred2s(out="sev.or.inte", bio="Fer", age=10, serotype="Others", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") dat_m10 <- get_pred2s(out="sev.or.inte", bio="CRP", age=10, serotype="Others", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Children") ## For adults dd$limits["Adjust to","Age"] <- 25 dat_mg1 <- get_pred2s(out="sev.or.inte", bio="VCAM", age=25, serotype="Others", dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg2 <- get_pred2s(out="sev.or.inte", bio="SDC", age=25, serotype="Others", dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg3 <- get_pred2s(out="sev.or.inte", bio="Ang", age=25, serotype="Others", dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg4 <- get_pred2s(out="sev.or.inte", bio="IL8", age=25, serotype="Others", dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg5 <- get_pred2s(out="sev.or.inte", bio="IP10", age=25, serotype="Others", dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg6 <- get_pred2s(out="sev.or.inte", bio="IL1RA", age=25, serotype="Others", dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg7 <- get_pred2s(out="sev.or.inte", bio="CD163", age=25, serotype="Others", dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg8 <- get_pred2s(out="sev.or.inte", bio="TREM", age=25, serotype="Others", dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg9 <- get_pred2s(out="sev.or.inte", bio="Fer", age=25, serotype="Others", dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_mg10 <- get_pred2s(out="sev.or.inte", bio="CRP", age=25, serotype="Others", dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$Serotype=='Others')) %>% mutate(age="Adults") dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(serotype = "Others") # Merge data for children and adults vis1 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 1") # Merge data for plots tmp0 <- dat1 %>% filter(!is.na(Serotype)) %>% arrange(Serotype) %>% select(sev.or.inte) tmp <- data.frame(sev.or.inte = rep(tmp0$sev.or.inte, 20)) tmp_p1 <- vis1 %>% arrange(biomarker, age) %>% mutate(value1 = 2^value, serotype = factor(serotype, levels=c("DENV-1", "Others")), age = factor(age, levels=c("Children", "Adults"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #devtools::install_github("zeehio/facetscales") library(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.125, .25, .5, 1, 2, 4, 8) ybreak2 <- c(.125, .25, .5, 1, 2, 4) scales_y1 <- list( `Children` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1), `Adults` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Children` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2), `Adults` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names age <- c(`Children` = "Children", `Adults` = "Adults") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=serotype), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=serotype), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, age=="Children" & sev.or.inte==0 & serotype=="DENV-1"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, age=="Children" & sev.or.inte==1 & serotype=="DENV-1"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, age=="Adults" & sev.or.inte==0 & serotype=="Others"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, age=="Adults" & sev.or.inte==1 & serotype=="Others"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + scale_fill_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(age), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(age, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=serotype), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=serotype), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, age=="Children" & sev.or.inte==0 & serotype=="DENV-1"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, age=="Children" & sev.or.inte==1 & serotype=="DENV-1"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, age=="Adults" & sev.or.inte==0 & serotype=="Others"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, age=="Adults" & sev.or.inte==1 & serotype=="Others"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + scale_fill_manual(values=c("red","blue"), labels=c("DENV-1","Others")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(age), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(age, biomarkers))) # Merge plots p52 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.86)) #ggsave(filename="A5_fig2.pdf", p52, dpi=300, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 5—table 1. Results from single models for severe/moderate dengue with the interaction with serotype. ```{r table A5.1, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # DENV-1 - children or1c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, serotype=1, dat=dat1), lo1c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, serotype=1, dat=dat1), up1c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, serotype=1, dat=dat1), # Others - children or1a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, serotype=2, dat=dat1), lo1a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, serotype=2, dat=dat1), up1a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, serotype=2, dat=dat1), # DENV-1 - adults or2c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, serotype=1, dat=dat1), lo2c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, serotype=1, dat=dat1), up2c = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, serotype=1, dat=dat1), # Others - adults or2a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, serotype=2, dat=dat1), lo2a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, serotype=2, dat=dat1), up2a = get_est1s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, serotype=2, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s1 = get_est1s(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), p.s1_int = get_est1s(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1), p.s1_int_age = get_est1s(out="sev.or.inte", bio=.$bio, est="p int age", dat=dat1), p.s1_int_serotype = get_est1s(out="sev.or.inte", bio=.$bio, est="p int serotype", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res1 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(is.na(ref1), as.character(bio), ifelse(bio!="Vir", paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), paste(" -", round(10^ref2,0), "vs", round(10^ref1,0), sep=" "))), or.sc1 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), or.sa1 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), or.sc2 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), or.sa2 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc1, or.sa1, or.sc2, or.sa2, p.s1, p.s1_int, p.s1_int_age, p.s1_int_serotype) names(res1) <- c("", "OR (children - DENV-1)", "OR (children - others)", "OR (adults - DENV-1)", "OR (adults - others)", "P overall", "P overall interaction", "P interaction with age", "P interaction with serotype") # Report results knitr::kable(res1) ``` ########################################################################################### #### Appendix 5—table 2. Results from global model for severe/moderate dengue with the interaction with serotype. ```{r table A5.2, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # DENV-1 - children or1c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, serotype=1, dat=dat1), lo1c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, serotype=1, dat=dat1), up1c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, serotype=1, dat=dat1), # Others - children or1a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, serotype=2, dat=dat1), lo1a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, serotype=2, dat=dat1), up1a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, serotype=2, dat=dat1), # DENV-1 - adults or2c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, serotype=1, dat=dat1), lo2c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, serotype=1, dat=dat1), up2c = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, serotype=1, dat=dat1), # Others - adults or2a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, serotype=2, dat=dat1), lo2a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, serotype=2, dat=dat1), up2a = get_est2s(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, serotype=2, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s2 = get_est2s(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), p.s2_int = get_est2s(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1), p.s2_int_age = get_est2s(out="sev.or.inte", bio=.$bio, est="p int age", dat=dat1), p.s2_int_serotype = get_est2s(out="sev.or.inte", bio=.$bio, est="p int serotype", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res2 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(is.na(ref1), as.character(bio), ifelse(bio!="Vir", paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), paste(" -", round(10^ref2,0), "vs", round(10^ref1,0), sep=" "))), or.sc1 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), or.sa1 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), or.sc2 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), or.sa2 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc1, or.sa1, or.sc2, or.sa2, p.s2, p.s2_int, p.s2_int_age, p.s2_int_serotype) names(res2) <- c("", "OR (children - DENV-1)", "OR (children - others)", "OR (adults - DENV-1)", "OR (adults - others)", "P overall", "P overall interaction", "P interaction with age", "P interaction with serotype") # Report results knitr::kable(res2) ``` ########################################################################################### #### Appendix 6—figure 1. Results from models for severe dengue endpoint. ```{r fig A6.1, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1(out="sev.only", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m2 <- get_pred1(out="sev.only", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m3 <- get_pred1(out="sev.only", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m4 <- get_pred1(out="sev.only", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m5 <- get_pred1(out="sev.only", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m6 <- get_pred1(out="sev.only", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m7 <- get_pred1(out="sev.only", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m8 <- get_pred1(out="sev.only", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m9 <- get_pred1(out="sev.only", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m10 <- get_pred1(out="sev.only", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg1 <- get_pred2(out="sev.only", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg2 <- get_pred2(out="sev.only", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg3 <- get_pred2(out="sev.only", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg4 <- get_pred2(out="sev.only", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg5 <- get_pred2(out="sev.only", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg6 <- get_pred2(out="sev.only", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg7 <- get_pred2(out="sev.only", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg8 <- get_pred2(out="sev.only", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg9 <- get_pred2(out="sev.only", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg10 <- get_pred2(out="sev.only", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "10 years") # For adults dd$limits["Adjust to","Age"] <- 25 dat_m1 <- get_pred1(out="sev.only", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m2 <- get_pred1(out="sev.only", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m3 <- get_pred1(out="sev.only", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m4 <- get_pred1(out="sev.only", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m5 <- get_pred1(out="sev.only", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m6 <- get_pred1(out="sev.only", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m7 <- get_pred1(out="sev.only", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m8 <- get_pred1(out="sev.only", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m9 <- get_pred1(out="sev.only", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m10 <- get_pred1(out="sev.only", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg1 <- get_pred2(out="sev.only", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg2 <- get_pred2(out="sev.only", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg3 <- get_pred2(out="sev.only", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg4 <- get_pred2(out="sev.only", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg5 <- get_pred2(out="sev.only", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg6 <- get_pred2(out="sev.only", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg7 <- get_pred2(out="sev.only", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg8 <- get_pred2(out="sev.only", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg9 <- get_pred2(out="sev.only", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg10 <- get_pred2(out="sev.only", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "25 years") # Merge results for children and adults vis2 <- rbind(dat_p1, dat_p2) %>% mutate(yhat = ifelse(yhat>log(32), NA, ifelse(yhat<log(1/32), NA, yhat)), lower = ifelse(lower>log(32), log(32), ifelse(lower<log(1/32), log(1/32), lower)), upper = ifelse(upper>log(32), log(32), ifelse(upper<log(1/32), log(1/32), upper)), outcome = "outcome 2") # Merge data for plots tmp0 <- dat %>% arrange(age15) %>% select(sev.only) tmp <- data.frame(sev.only = rep(tmp0$sev.only, 20)) tmp_p1 <- vis2 %>% mutate(value1 = 2^value, age = factor(age, levels=c("10 years", "25 years")), gr = ifelse(model=="Single model" & age=="10 years", 1, ifelse(model=="Single model" & age=="25 years", 2, ifelse(model=="Global model" & age=="10 years", 3, 4))), gr = factor(gr, levels=c(1:4), labels=c("Single children", "Single adults", "Global children", "Global adults")), model = factor(model, levels=c("Single model", "Global model"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(model=="Single model") %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #library(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.0625,.125, .25, .5, 1, 2, 4, 8, 16) ybreak2 <- c(.0625,.125, .25, .5, 1, 2, 4, 8, 16) scales_y1 <- list( `Single model` = scale_y_continuous(limits = c(log(.0624), log(16.01)), breaks = log(ybreak1), labels = ybreak1), `Global model` = scale_y_continuous(limits = c(log(.0624), log(16.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Single model` = scale_y_continuous(limits = c(log(.0624), log(16.01)), breaks = log(ybreak2), labels = ybreak2), `Global model` = scale_y_continuous(limits = c(log(.0624), log(16.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names models <- c(`Single model` = "Single model", `Global model` = "Global model") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.only==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.only==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.only==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.only==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(models, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.only==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.only==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.only==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.only==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(models, biomarkers))) # Merge plots p61 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.85)) #ggsave(filename="A6_fig1.pdf", p61, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 6—figure 2. Results from models for severe dengue or dengue with warning signs endpoint. ```{r fig A6.2, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1(out="sev.or.ws", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m2 <- get_pred1(out="sev.or.ws", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m3 <- get_pred1(out="sev.or.ws", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m4 <- get_pred1(out="sev.or.ws", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m5 <- get_pred1(out="sev.or.ws", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m6 <- get_pred1(out="sev.or.ws", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m7 <- get_pred1(out="sev.or.ws", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m8 <- get_pred1(out="sev.or.ws", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m9 <- get_pred1(out="sev.or.ws", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m10 <- get_pred1(out="sev.or.ws", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg1 <- get_pred2(out="sev.or.ws", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg2 <- get_pred2(out="sev.or.ws", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg3 <- get_pred2(out="sev.or.ws", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg4 <- get_pred2(out="sev.or.ws", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg5 <- get_pred2(out="sev.or.ws", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg6 <- get_pred2(out="sev.or.ws", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg7 <- get_pred2(out="sev.or.ws", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg8 <- get_pred2(out="sev.or.ws", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg9 <- get_pred2(out="sev.or.ws", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg10 <- get_pred2(out="sev.or.ws", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "10 years") # For adults dd$limits["Adjust to","Age"] <- 25 dat_m1 <- get_pred1(out="sev.or.ws", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m2 <- get_pred1(out="sev.or.ws", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m3 <- get_pred1(out="sev.or.ws", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m4 <- get_pred1(out="sev.or.ws", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m5 <- get_pred1(out="sev.or.ws", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m6 <- get_pred1(out="sev.or.ws", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m7 <- get_pred1(out="sev.or.ws", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m8 <- get_pred1(out="sev.or.ws", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m9 <- get_pred1(out="sev.or.ws", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m10 <- get_pred1(out="sev.or.ws", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg1 <- get_pred2(out="sev.or.ws", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg2 <- get_pred2(out="sev.or.ws", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg3 <- get_pred2(out="sev.or.ws", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg4 <- get_pred2(out="sev.or.ws", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg5 <- get_pred2(out="sev.or.ws", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg6 <- get_pred2(out="sev.or.ws", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg7 <- get_pred2(out="sev.or.ws", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg8 <- get_pred2(out="sev.or.ws", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg9 <- get_pred2(out="sev.or.ws", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg10 <- get_pred2(out="sev.or.ws", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "25 years") # Merge data for children and adults vis3 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 3") # Merge data for plots tmp0 <- dat %>% arrange(age15) %>% select(sev.or.ws) tmp <- data.frame(sev.or.ws = rep(tmp0$sev.or.ws, 20)) tmp_p1 <- vis3 %>% mutate(value1 = 2^value, age = factor(age, levels=c("10 years", "25 years")), gr = ifelse(model=="Single model" & age=="10 years", 1, ifelse(model=="Single model" & age=="25 years", 2, ifelse(model=="Global model" & age=="10 years", 3, 4))), gr = factor(gr, levels=c(1:4), labels=c("Single children", "Single adults", "Global children", "Global adults")), model = factor(model, levels=c("Single model", "Global model"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(model=="Single model") %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #library(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.25, .5, 1, 2, 4) ybreak2 <- c(.25, .5, 1, 2, 4) scales_y1 <- list( `Single model` = scale_y_continuous(limits = c(log(.24), log(4.01)), breaks = log(ybreak1), labels = ybreak1), `Global model` = scale_y_continuous(limits = c(log(.24), log(4.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Single model` = scale_y_continuous(limits = c(log(.24), log(4.01)), breaks = log(ybreak2), labels = ybreak2), `Global model` = scale_y_continuous(limits = c(log(.24), log(4.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names models <- c(`Single model` = "Single model", `Global model` = "Global model") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.ws==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.ws==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.ws==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.ws==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(models, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.or.ws==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Single model" & sev.or.ws==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.or.ws==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, model=="Global model" & sev.or.ws==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(models, biomarkers))) # Merge plots p62 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.85)) #ggsave(filename="A6_fig2.pdf", p62, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 6—figure 3. Results from models for hospitalization endpoint. ```{r fig A6.3, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1(out="hospital", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m2 <- get_pred1(out="hospital", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m3 <- get_pred1(out="hospital", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m4 <- get_pred1(out="hospital", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m5 <- get_pred1(out="hospital", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m6 <- get_pred1(out="hospital", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m7 <- get_pred1(out="hospital", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m8 <- get_pred1(out="hospital", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m9 <- get_pred1(out="hospital", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m10 <- get_pred1(out="hospital", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg1 <- get_pred2(out="hospital", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg2 <- get_pred2(out="hospital", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg3 <- get_pred2(out="hospital", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg4 <- get_pred2(out="hospital", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg5 <- get_pred2(out="hospital", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg6 <- get_pred2(out="hospital", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg7 <- get_pred2(out="hospital", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg8 <- get_pred2(out="hospital", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg9 <- get_pred2(out="hospital", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg10 <- get_pred2(out="hospital", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "10 years") # For adults dd$limits["Adjust to","Age"] <- 25 dat_m1 <- get_pred1(out="hospital", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m2 <- get_pred1(out="hospital", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m3 <- get_pred1(out="hospital", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m4 <- get_pred1(out="hospital", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m5 <- get_pred1(out="hospital", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m6 <- get_pred1(out="hospital", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m7 <- get_pred1(out="hospital", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m8 <- get_pred1(out="hospital", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m9 <- get_pred1(out="hospital", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m10 <- get_pred1(out="hospital", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg1 <- get_pred2(out="hospital", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg2 <- get_pred2(out="hospital", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg3 <- get_pred2(out="hospital", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg4 <- get_pred2(out="hospital", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg5 <- get_pred2(out="hospital", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg6 <- get_pred2(out="hospital", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg7 <- get_pred2(out="hospital", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg8 <- get_pred2(out="hospital", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg9 <- get_pred2(out="hospital", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg10 <- get_pred2(out="hospital", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10) %>% mutate(age = "25 years") # Merge data for children and adults vis4 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 4") # Merge data for plots tmp0 <- dat %>% arrange(age15) %>% select(hospital) tmp <- data.frame(hospital = rep(tmp0$hospital, 20)) tmp_p1 <- vis4 %>% mutate(value1 = 2^value, age = factor(age, levels=c("10 years", "25 years")), gr = ifelse(model=="Single model" & age=="10 years", 1, ifelse(model=="Single model" & age=="25 years", 2, ifelse(model=="Global model" & age=="10 years", 3, 4))), gr = factor(gr, levels=c(1:4), labels=c("Single children", "Single adults", "Global children", "Global adults")), model = factor(model, levels=c("Single model", "Global model"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(model=="Single model") %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_10 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP"))) # Modify facets' scales #library(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP) ) ybreak1 <- c(.125, .25, .5, 1, 2, 4, 8) ybreak2 <- c(.125, .25, .5, 1, 2, 4, 8) scales_y1 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1), `Global model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak2), labels = ybreak2), `Global model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names models <- c(`Single model` = "Single model", `Global model` = "Global model") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, model=="Single model" & hospital==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Single model" & hospital==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Global model" & hospital==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, model=="Global model" & hospital==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(models, biomarkers))) # Plot for the last 5 biomarkers p.6.10 <- ggplot(dat_6_10, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_10, model=="Single model" & hospital==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Single model" & hospital==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_10, model=="Global model" & hospital==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_10, model=="Global model" & hospital==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(models, biomarkers))) # Merge plots p63 <- gridExtra::grid.arrange(p.1.5, p.6.10, nrow=2, heights=c(1,.85)) #ggsave(filename="A6_fig3.pdf", p63, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 6—table 1. Results from models for severe dengue endpoint. ```{r table A6.1, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # Children - single models or1c = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo1c = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up1c = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Children - global model or2c = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo2c = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up2c = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Adults - single models or1a = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo1a = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up1a = get_est1(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1), # Adults - global model or2a = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo2a = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up2a = get_est2(out="sev.only", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s2 = get_est1(out="sev.only", bio=.$bio, est="p", dat=dat1), p.s2_int = get_est1(out="sev.only", bio=.$bio, est="p int", dat=dat1), p.g2 = get_est2(out="sev.only", bio=.$bio, est="p", dat=dat1), p.g2_int = get_est2(out="sev.only", bio=.$bio, est="p int", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res2 <- tmp1 %>% group_by(bio) %>% slice(1) %>% left_join(., select(tmp2, -c(sort1, sort2)), by="bio") %>% arrange(sort1) %>% mutate(bio = as.character(bio), or.sc2 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), # s: single model; c: children or.gc2 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), # g: global model or.sa2 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), # a: adults or.ga2 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc2, or.sa2, p.s2, p.s2_int, or.gc2, or.ga2, p.g2, p.g2_int) names(res2) <- c("", "OR (single - children)", "OR (single - adults)", "P overall (single)", "P interaction (single)", "OR (global - children)", "OR (global - adults)", "P overall (global)", "P interaction (global)") # Report results knitr::kable(res2) ``` ########################################################################################### #### Appendix 6—table 2. Results from models for severe dengue or dengue with warning signs endpoint. ```{r table A6.2, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # Children - single models or1c = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo1c = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up1c = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Children - global model or2c = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo2c = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up2c = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Adults - single models or1a = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo1a = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up1a = get_est1(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1), # Adults - global model or2a = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo2a = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up2a = get_est2(out="sev.or.ws", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s3 = get_est1(out="sev.or.ws", bio=.$bio, est="p", dat=dat1), p.s3_int = get_est1(out="sev.or.ws", bio=.$bio, est="p int", dat=dat1), p.g3 = get_est2(out="sev.or.ws", bio=.$bio, est="p", dat=dat1), p.g3_int = get_est2(out="sev.or.ws", bio=.$bio, est="p int", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res3 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(!is.na(ref1), paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), as.character(bio)), or.sc3 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), # s: single model; c: children or.gc3 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), # g: global model or.sa3 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), # a: adults or.ga3 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc3, or.sa3, p.s3, p.s3_int, or.gc3, or.ga3, p.g3, p.g3_int) names(res3) <- c("", "OR (single - children)", "OR (single - adults)", "P overall (single)", "P interaction (single)", "OR (global - children)", "OR (global - adults)", "P overall (global)", "P interaction (global)") # Report results knitr::kable(res3) ``` ########################################################################################### #### Appendix 6—table 3. Results from models for hospitalization endpoint. ```{r table A6.3, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:10), 2), sort2 = c(rep(2,10), rep(3,10)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP"), 2), ref1 = c(ref0-1, ref0), ref2 = c(ref0, ref0+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # Children - single models or1c = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo1c = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up1c = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Children - global model or2c = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo2c = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up2c = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Adults - single models or1a = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo1a = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up1a = get_est1(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1), # Adults - global model or2a = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo2a = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up2a = get_est2(out="hospital", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:10), sort2 = rep(1,10), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") ) %>% group_by(sort1) %>% do(cbind(., p.s4 = get_est1(out="hospital", bio=.$bio, est="p", dat=dat1), p.s4_int = get_est1(out="hospital", bio=.$bio, est="p int", dat=dat1), p.g4 = get_est2(out="hospital", bio=.$bio, est="p", dat=dat1), p.g4_int = get_est2(out="hospital", bio=.$bio, est="p int", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res4 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(!is.na(ref1), paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), as.character(bio)), or.sc4 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), # s: single model; c: children or.gc4 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), # g: global model or.sa4 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), # a: adults or.ga4 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc4, or.sa4, p.s4, p.s4_int, or.gc4, or.ga4, p.g4, p.g4_int) names(res4) <- c("", "OR (single - children)", "OR (single - adults)", "P overall (single)", "P interaction (single)", "OR (global - children)", "OR (global - adults)", "P overall (global)", "P interaction (global)") # Report results knitr::kable(res4) ``` ########################################################################################### #### Appendix 7—table 1. Model selection frequencies for children. ```{r table A7.1, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Bootstrap results ------------------------------------------------ ## The following codes were modified from Heinze G. et al (https://doi.org/10.1002/bimj.201700067) bootnum <- 1000 load("bootstrap_results.RData") # Load bootstrap results # Model frequency group_cols <- c(colnames(data.frame(boot_varc))[1:length(pred)]) boot_modfreq <- data.frame(boot_varc) %>% mutate(time = 1) %>% group_by_(.dots = group_cols) %>% summarise(aic_median = round(median(aic),1), aic_025 = round(quantile(aic, 0.025),1), aic_975 = round(quantile(aic, 0.975),1), count = sum(time)) %>% ungroup() %>% mutate(percent = count/bootnum * 100) %>% arrange(desc(count)) %>% slice(1:20) %>% as.data.frame(.) out <- boot_modfreq for (i in 1:(ncol(out)-5)) { out[,i] <- ifelse(out[,i]==0, NA, "+") } knitr::kable(out) ``` ########################################################################################### #### Appendix 7—table 2. Model stability for children. ```{r table A7.2, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + IP10 + IL1RA + CD163 + TREM + Fer + CRP, family=binomial, data=dat1c, x=T, y=T) full_est <- coef(full_mod) full_se <- coef(summary(full_mod))[, "Std. Error"] pred_name <- names(full_est)[-1] # Selected model (for bootstrap results) --------------------------- sel_m <- dredge(full_mod, rank="AIC") sel_var <- matrix(0, ncol = length(pred), nrow = 1, dimnames = list(NULL, pred)) sel_var_tmp <- attr(get.models(sel_m, 1)[[1]]$terms, "term.labels") for (j in 1:ncol(sel_var)) { sel_var[1,j] <- ifelse(names(sel_var[1,j]) %in% sel_var_tmp, 1, 0) } formula <- paste("sev.or.inte~", paste(names(sel_var[1,][sel_var[1,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1, family = binomial, x = T, y = T) sel_est <- coef(sel_mod[1])[c("(Intercept)", pred_name)] sel_est[is.na(sel_est)] <- 0 names(sel_est) <- c("(Intercept)", pred_name) sel_se <- coef(summary(sel_mod))[, "Std. Error"][c("(Intercept)", pred_name)] sel_se[is.na(sel_se)] <- 0 # Bootstrap results ------------------------------------------------ ## The following codes were modified from Heinze G. et al (https://doi.org/10.1002/bimj.201700067) bootnum <- 1000 load("bootstrap_results.RData") # Load bootstrap results boot_01 <- (boot_estc != 0) * 1 boot_inclusion <- apply(boot_01, 2, function(x) sum(x) / length(x) * 100) ## Overview of estimates and measures sqe <- (t(boot_estc) - full_est) ^ 2 rmsd <- apply(sqe, 1, function(x) sqrt(mean(x))) rmsdratio <- rmsd / full_se boot_mean <- apply(boot_estc, 2, mean) boot_meanratio <- boot_mean / full_est boot_relbias <- (boot_meanratio / (boot_inclusion / 100) - 1) * 100 boot_median <- apply(boot_estc, 2, median) boot_025per <- apply(boot_estc, 2, function(x) quantile(x, 0.025)) boot_975per <- apply(boot_estc, 2, function(x) quantile(x, 0.975)) overview <- round(cbind(full_est, full_se, boot_inclusion, sel_est, sel_se, rmsdratio, boot_relbias, boot_median, boot_025per, boot_975per), 4) overview <- overview[order(overview[,"boot_inclusion"], decreasing=T),] knitr::kable(overview) ``` ########################################################################################### #### Appendix 7—table 3. Model selection frequencies for adults. ```{r table A7.3, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "ns1(IP10)", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Bootstrap results ------------------------------------------------ ## The following codes were modified from Heinze G. et al (https://doi.org/10.1002/bimj.201700067) bootnum <- 1000 load("bootstrap_results.RData") # Load bootstrap results ## Model frequency group_cols <- c(colnames(data.frame(boot_vara))[1:length(pred)]) boot_modfreq <- data.frame(boot_vara) %>% mutate(time = 1) %>% group_by_(.dots = group_cols) %>% summarise(aic_median = round(median(aic),1), aic_025 = round(quantile(aic, 0.025),1), aic_975 = round(quantile(aic, 0.975),1), count = sum(time)) %>% ungroup() %>% mutate(percent = count/bootnum * 100) %>% arrange(desc(count)) %>% slice(1:20) %>% as.data.frame(.) out <- boot_modfreq for (i in 1:(ncol(out)-5)) { out[,i] <- ifelse(out[,i]==0, NA, "+") } knitr::kable(out) ``` ########################################################################################### #### Appendix 7—table 4. Model stability for adults. ```{r table A7.4, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "ns1(IP10)", "IL1RA", "CD163", "TREM", "Fer", "CRP") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + ns1(IP10) + IL1RA + CD163 + TREM + Fer + CRP, family=binomial, data=dat1a, x=T, y=T) full_est <- coef(full_mod) full_se <- coef(summary(full_mod))[, "Std. Error"] pred_name <- names(full_est)[-1] # Selected model (for bootstrap results) --------------------------- sel_m <- dredge(full_mod, rank="AIC") sel_var <- matrix(0, ncol = length(pred), nrow = 1, dimnames = list(NULL, pred)) sel_var_tmp <- attr(get.models(sel_m, 1)[[1]]$terms, "term.labels") for (j in 1:ncol(sel_var)) { sel_var[1,j] <- ifelse(names(sel_var[1,j]) %in% sel_var_tmp, 1, 0) } formula <- paste("sev.or.inte~", paste(names(sel_var[1,][sel_var[1,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1, family = binomial, x = T, y = T) sel_est <- coef(sel_mod[1])[c("(Intercept)", pred_name)] sel_est[is.na(sel_est)] <- 0 names(sel_est) <- c("(Intercept)", pred_name) sel_se <- coef(summary(sel_mod))[, "Std. Error"][c("(Intercept)", pred_name)] sel_se[is.na(sel_se)] <- 0 # Bootstrap results ------------------------------------------------ ## The following codes were modified from Heinze G. et al (https://doi.org/10.1002/bimj.201700067) bootnum <- 1000 load("bootstrap_results.RData") # Load bootstrap results boot_01 <- (boot_esta != 0) * 1 boot_inclusion <- apply(boot_01, 2, function(x) sum(x) / length(x) * 100) ## Overview of estimates and measures sqe <- (t(boot_esta) - full_est) ^ 2 rmsd <- apply(sqe, 1, function(x) sqrt(mean(x))) rmsdratio <- rmsd / full_se boot_mean <- apply(boot_esta, 2, mean) boot_meanratio <- boot_mean / full_est boot_relbias <- (boot_meanratio / (boot_inclusion / 100) - 1) * 100 boot_median <- apply(boot_esta, 2, median) boot_025per <- apply(boot_esta, 2, function(x) quantile(x, 0.025)) boot_975per <- apply(boot_esta, 2, function(x) quantile(x, 0.975)) overview <- round(cbind(full_est, full_se, boot_inclusion, sel_est, sel_se, rmsdratio, boot_relbias, boot_median, boot_025per, boot_975per), 4) overview <- overview[order(overview[,"boot_inclusion"], decreasing=T),] knitr::kable(overview) ``` ########################################################################################### #### Appendix 8—figure 1. Results from models for severe/moderate dengue including viremia as a potential biomarker. ```{r fig A8.1, echo=TRUE, message=FALSE, warning=FALSE, fig.height=6.2, fig.width=8.5} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # For children dd$limits["Adjust to","Age"] <- 10 dat_m1 <- get_pred1(out="sev.or.inte", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m2 <- get_pred1(out="sev.or.inte", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m3 <- get_pred1(out="sev.or.inte", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m4 <- get_pred1(out="sev.or.inte", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m5 <- get_pred1(out="sev.or.inte", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m6 <- get_pred1(out="sev.or.inte", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m7 <- get_pred1(out="sev.or.inte", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m8 <- get_pred1(out="sev.or.inte", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m9 <- get_pred1(out="sev.or.inte", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m10 <- get_pred1(out="sev.or.inte", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_m11 <- get_pred1(out="sev.or.inte", bio="Vir", age=10, dat=dat1) %>% rename(value=Vir) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg1 <- get_pred2v(out="sev.or.inte", bio="VCAM", age=10, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg2 <- get_pred2v(out="sev.or.inte", bio="SDC", age=10, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg3 <- get_pred2v(out="sev.or.inte", bio="Ang", age=10, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg4 <- get_pred2v(out="sev.or.inte", bio="IL8", age=10, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg5 <- get_pred2v(out="sev.or.inte", bio="IP10", age=10, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg6 <- get_pred2v(out="sev.or.inte", bio="IL1RA", age=10, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg7 <- get_pred2v(out="sev.or.inte", bio="CD163", age=10, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg8 <- get_pred2v(out="sev.or.inte", bio="TREM", age=10, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg9 <- get_pred2v(out="sev.or.inte", bio="Fer", age=10, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg10 <- get_pred2v(out="sev.or.inte", bio="CRP", age=10, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_mg11 <- get_pred2v(out="sev.or.inte", bio="Vir", age=10, dat=dat1) %>% rename(value=Vir) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="Under 15")) dat_p1 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10, dat_m11,dat_mg11) %>% mutate(age = "10 years") # For adults dd$limits["Adjust to","Age"] <- 25 dat_m1 <- get_pred1(out="sev.or.inte", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m2 <- get_pred1(out="sev.or.inte", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m3 <- get_pred1(out="sev.or.inte", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m4 <- get_pred1(out="sev.or.inte", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m5 <- get_pred1(out="sev.or.inte", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m6 <- get_pred1(out="sev.or.inte", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m7 <- get_pred1(out="sev.or.inte", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m8 <- get_pred1(out="sev.or.inte", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m9 <- get_pred1(out="sev.or.inte", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m10 <- get_pred1(out="sev.or.inte", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_m11 <- get_pred1(out="sev.or.inte", bio="Vir", age=25, dat=dat1) %>% rename(value=Vir) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg1 <- get_pred2v(out="sev.or.inte", bio="VCAM", age=25, dat=dat1) %>% rename(value=VCAM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg2 <- get_pred2v(out="sev.or.inte", bio="SDC", age=25, dat=dat1) %>% rename(value=SDC) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg3 <- get_pred2v(out="sev.or.inte", bio="Ang", age=25, dat=dat1) %>% rename(value=Ang) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg4 <- get_pred2v(out="sev.or.inte", bio="IL8", age=25, dat=dat1) %>% rename(value=IL8) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg5 <- get_pred2v(out="sev.or.inte", bio="IP10", age=25, dat=dat1) %>% rename(value=IP10) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg6 <- get_pred2v(out="sev.or.inte", bio="IL1RA", age=25, dat=dat1) %>% rename(value=IL1RA) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg7 <- get_pred2v(out="sev.or.inte", bio="CD163", age=25, dat=dat1) %>% rename(value=CD163) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg8 <- get_pred2v(out="sev.or.inte", bio="TREM", age=25, dat=dat1) %>% rename(value=TREM) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg9 <- get_pred2v(out="sev.or.inte", bio="Fer", age=25, dat=dat1) %>% rename(value=Fer) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg10 <- get_pred2v(out="sev.or.inte", bio="CRP", age=25, dat=dat1) %>% rename(value=CRP) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_mg11 <- get_pred2v(out="sev.or.inte", bio="Vir", age=25, dat=dat1) %>% rename(value=Vir) %>% select(value, yhat, lower, upper, biomarker, model) %>% slice(which(dat1$age15=="15 and above")) dat_p2 <- rbind(dat_m1,dat_mg1, dat_m2,dat_mg2, dat_m3,dat_mg3, dat_m4,dat_mg4, dat_m5,dat_mg5, dat_m6,dat_mg6, dat_m7,dat_mg7, dat_m8,dat_mg8, dat_m9,dat_mg9, dat_m10,dat_mg10, dat_m11,dat_mg11) %>% mutate(age = "25 years") # Merge data for children and adults vis1 <- rbind(dat_p1, dat_p2) %>% mutate(outcome = "outcome 1") # Merge data for plots tmp0 <- dat %>% arrange(age15) %>% select(sev.or.inte) tmp <- data.frame(sev.or.inte = rep(tmp0$sev.or.inte, 22)) tmp_p1 <- vis1 %>% arrange(biomarker, model) %>% mutate(value1 = 2^value, age = factor(age, levels=c("10 years", "25 years")), gr = ifelse(model=="Single model" & age=="10 years", 1, ifelse(model=="Single model" & age=="25 years", 2, ifelse(model=="Global model" & age=="10 years", 3, 4))), gr = factor(gr, levels=c(1:4), labels=c("Single children", "Single adults", "Global children", "Global adults")), model = factor(model, levels=c("Single model", "Global model"))) %>% bind_cols(., tmp) vline <- tmp_p1 %>% filter(model=="Single model") %>% filter(yhat==0) %>% select(biomarker, value, value1) %>% rename(vline=value, vline1=value1) dat_p <- left_join(tmp_p1, vline, by="biomarker") # Alternative data to limit to 5th - 95th of the values dat_alt <- dat_p %>% filter(value != 0) %>% group_by(biomarker, model) %>% mutate(up = quantile(value, .95), lo = quantile(value, .05)) %>% ungroup() %>% mutate(is.outlier = value<lo | value>up | value<0, value = ifelse(is.outlier, NA, value), value1 = ifelse(is.outlier, NA, value1), yhat = ifelse(is.outlier, NA, yhat), lower = ifelse(is.outlier, NA, lower), upper = ifelse(is.outlier, NA, upper)) dat_1_5 <- dat_alt %>% filter(biomarker %in% c("VCAM", "SDC", "Ang", "IL8", "IP10")) %>% mutate(biomarker = factor(biomarker, levels=c("VCAM", "SDC", "Ang", "IL8", "IP10"))) dat_6_11 <- dat_alt %>% filter(biomarker %in% c("IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir")) %>% mutate(biomarker = factor(biomarker, levels=c("IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir"))) # Modify facets' scales #library(facetscales) xVCAM <- c(1,4,15,60,250,1000,4000) xSDC <- c(1400,2000,2800,4000,5600) xAng <- c(50,100,250,500,1000,2000) xIL8 <- c(5,7,10,14,20,28,40) xIP10 <- c(25,100,400,1600,6400) xIL1RA <- c(1000,2000,4000,8000,16000) xCD163 <- c(75,150,300,600) xTREM <- c(35,50,70,100,140,200) xFer <- c(50,100,200,400,800) xCRP <- c(2.5,5,10,20,40,80) xVir <- c(0:10) scales_x <- list( `VCAM` = scale_x_continuous(breaks = log2(xVCAM), labels = xVCAM), `SDC` = scale_x_continuous(breaks = log2(xSDC), labels = xSDC), `Ang` = scale_x_continuous(breaks = log2(xAng), labels = xAng), `IL8` = scale_x_continuous(breaks = log2(xIL8), labels = xIL8), `IP10` = scale_x_continuous(breaks = log2(xIP10), labels = xIP10), `IL1RA` = scale_x_continuous(breaks = log2(xIL1RA), labels = xIL1RA), `CD163` = scale_x_continuous(breaks = log2(xCD163), labels = xCD163), `TREM` = scale_x_continuous(breaks = log2(xTREM), labels = xTREM), `Fer` = scale_x_continuous(breaks = log2(xFer), labels = xFer), `CRP` = scale_x_continuous(breaks = log2(xCRP), labels = xCRP), `Vir` = scale_x_continuous(breaks = xVir, labels = xVir) ) ybreak1 <- c(.125, .25, .5, 1, 2, 4, 8) ybreak2 <- c(.125, .25, .5, 1, 2, 4) scales_y1 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1), `Global model` = scale_y_continuous(limits = c(log(.124), log(8.01)), breaks = log(ybreak1), labels = ybreak1) ) scales_y2 <- list( `Single model` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2), `Global model` = scale_y_continuous(limits = c(log(.124), log(4.01)), breaks = log(ybreak2), labels = ybreak2) ) # Set facets' names models <- c(`Single model` = "Single model", `Global model` = "Global model") biomarkers <- c(`VCAM` = "VCAM-1 (ng/ml)", `SDC`= "SDC-1 (pg/ml)", `Ang` = "Ang-2 (pg/ml)", `IL8` = "IL-8 (pg/ml)", `IP10` = "IP-10 (pg/ml)", `IL1RA` = "IL-1RA (pg/ml)", `CD163` = "sCD163 (ng/ml)", `TREM` = "sTREM-1 (pg/ml)", `Fer` = "Ferritin (ng/ml)", `CRP` = "CRP (mg/l)", `Vir` = "Log10 viremia") # Plot for the first 5 biomarkers p.1.5 <- ggplot(dat_1_5, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.inte==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Single model" & sev.or.inte==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.inte==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_1_5, model=="Global model" & sev.or.inte==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="top", legend.title=element_blank(), axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y1), labeller = as_labeller(c(models, biomarkers))) # Plot for the last 5 biomarkers p.6.11 <- ggplot(dat_6_11, aes(x=value)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=age), alpha=.15) + geom_vline(aes(xintercept=vline), color="black", linetype="dashed", alpha=.4) + geom_hline(yintercept=0, color="black", linetype="dashed", alpha=.4) + geom_line(aes(y=yhat, color=age), size=.5, alpha=.7) + geom_rug(data=filter(dat_6_11, model=="Single model" & sev.or.inte==0 & age=="10 years"), sides="b", alpha=.2, color="red") + geom_rug(data=filter(dat_6_11, model=="Single model" & sev.or.inte==1 & age=="10 years"), sides="t", alpha=.2, color="red") + geom_rug(data=filter(dat_6_11, model=="Global model" & sev.or.inte==0 & age=="25 years"), sides="b", alpha=.2, color="blue") + geom_rug(data=filter(dat_6_11, model=="Global model" & sev.or.inte==1 & age=="25 years"), sides="t", alpha=.2, color="blue") + scale_color_manual(values=c("red","blue"), labels=c("Children","Adults")) + scale_fill_manual(values=c("red","blue"), labels=c("Children","Adults")) + ylab("Odds ratio") + theme(legend.position="none", axis.title.x=element_blank(), axis.text.x=element_text(angle=45)) + facet_grid_sc(cols=vars(biomarker), rows=vars(model), scales=list(x=scales_x, y=scales_y2), labeller = as_labeller(c(models, biomarkers))) # Merge plots p81 <- gridExtra::grid.arrange(p.1.5, p.6.11, nrow=2, heights=c(1,.86)) #ggsave(filename="A8_fig1.pdf", p81, dpi=300, width=8.5, height=6.2) ``` ########################################################################################### #### Appendix 8—table 1. Results from models for severe/moderate dengue including viremia as a potential biomarker. ```{r table A8.1, echo=TRUE, message=FALSE, warning=FALSE} # Set datadist for 'lrm' function [rms] dd <- datadist(dat1); options(datadist="dd") # Set references to calculate results ref0v <- c(ref0, median(dat1$Vir)) # Get ORs and 95% CIs tmp1 <- data.frame( sort1 = rep(c(1:11), 2), sort2 = c(rep(2,11), rep(3,11)), bio = rep(c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir"), 2), ref1 = c(ref0v-1, ref0v), ref2 = c(ref0v, ref0v+1) ) %>% arrange(sort1, sort2) %>% group_by(sort1, sort2) %>% do(cbind(., # Children - single models or1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up1c = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Children - global model or2c = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=10, dat=dat1), lo2c = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=10, dat=dat1), up2c = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=10, dat=dat1), # Adults - single models or1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up1a = get_est1(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1), # Adults - global model or2a = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="OR", age=25, dat=dat1), lo2a = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="loCI", age=25, dat=dat1), up2a = get_est2v(out="sev.or.inte", bio=.$bio, ref1=.$ref1, ref2=.$ref2, est="upCI", age=25, dat=dat1))) %>% ungroup() for (i in 6:17) {tmp1[[i]] <- sprintf("%.2f", round(tmp1[[i]],2))} # Get p-values tmp2 <- data.frame( sort1 = c(1:11), sort2 = rep(1,11), bio = c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir") ) %>% group_by(sort1) %>% do(cbind(., p.s1 = get_est1(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), p.s1_int = get_est1(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1), p.g1 = get_est2v(out="sev.or.inte", bio=.$bio, est="p", dat=dat1), p.g1_int = get_est2v(out="sev.or.inte", bio=.$bio, est="p int", dat=dat1))) %>% ungroup() for (i in 4:7) {tmp2[[i]] <- ifelse(tmp2[[i]]<0.001, "<0.001", sprintf("%.3f", round(tmp2[[i]],3)))} # Combine results res1 <- bind_rows(tmp1, tmp2) %>% arrange(sort1, sort2) %>% mutate(bio = ifelse(is.na(ref1), as.character(bio), ifelse(bio!="Vir", paste(" -", round(2^ref2,0), "vs", round(2^ref1,0), sep=" "), paste(" -", round(ref2,1), "vs", round(ref1,1), sep=" "))), or.sc1 = ifelse(is.na(lo1c), NA, paste(or1c, " (", lo1c, "-", up1c, ")", sep="")), # s: single model; c: children or.gc1 = ifelse(is.na(lo2c), NA, paste(or2c, " (", lo2c, "-", up2c, ")", sep="")), # g: global model or.sa1 = ifelse(is.na(lo1a), NA, paste(or1a, " (", lo1a, "-", up1a, ")", sep="")), # a: adults or.ga1 = ifelse(is.na(lo2a), NA, paste(or2a, " (", lo2a, "-", up2a, ")", sep=""))) %>% select(bio, or.sc1, or.sa1, p.s1, p.s1_int, or.gc1, or.ga1, p.g1, p.g1_int) names(res1) <- c("", "OR (single - children)", "OR (single - adults)", "P overall (single)", "P interaction (single)", "OR (global - children)", "OR (global - adults)", "P overall (global)", "P interaction (global)") # Report results knitr::kable(res1) ``` ########################################################################################### #### Appendix 8—table 2. Best combinations of biomarkers associated with severe or moderate dengue for children. ```{r table A8.2, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + IP10 + IL1RA + CD163 + TREM + Fer + CRP + Vir, family=binomial, data=dat1c, x=T, y=T) # Selected model --------------------------------------------------- sel_var <- matrix(0, ncol=length(pred)+1, nrow=5, dimnames=list(NULL, c(pred, "aic"))) for (i in 1:5) { if (i==1) { bs <- dredge(full_mod, rank="AIC") } else { bs <- dredge(full_mod, rank="AIC", m.lim=c(i,i)) } bs_var <- attr(get.models(bs, 1)[[1]]$terms, "term.labels") for (j in 1:(ncol(sel_var)-1)) {sel_var[i,j] <- ifelse(names(sel_var[i,j]) %in% bs_var, 1, 0)} formula <- paste("sev.or.inte~", paste(names(sel_var[i,][sel_var[i,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1c, family = binomial, x = T, y = T) sel_var[i, ncol(sel_var)] <- AIC(sel_mod) } # Report results -------------------------------------------------- out1 <- as.data.frame(sel_var) %>% mutate(AIC = round(aic,1)) %>% select(-aic) for (i in 1:(ncol(out1)-1)) {out1[,i] <- ifelse(out1[,i]==0, NA, "+")} out2 <- as.data.frame(t(out1)) colnames(out2) <- c("Best of all combinations", "Best combination of 2 variables", "Best combination of 3 variables", "Best combination of 4 variables", "Best combination of 5 variables") rownames(out2) <- c("- VCAM-1", "- SDC-1", "- Ang-2", "- IL-8", "- IP-10", "- IL-1RA", "- sCD163", "- sTREM-1", "- Ferritin", "- CRP", "- Viremia", "AIC of the selected model") knitr::kable(out2) ``` ########################################################################################### #### Appendix 8—table 3. Best combinations of biomarkers associated with severe or moderate dengue for adults. ```{r table A8.3, echo=TRUE, message=FALSE, warning=FALSE} # EPV -------------------------------------------------------------- pred <- c("VCAM", "SDC", "Ang", "IL8", "ns1(IP10)", "IL1RA", "CD163", "TREM", "Fer", "CRP", "Vir") # Estimate full model ---------------------------------------------- full_mod <- glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + ns1(IP10) + IL1RA + CD163 + TREM + Fer + CRP + Vir, family=binomial, data=dat1a, x=T, y=T) # Selected model --------------------------------------------------- sel_var <- matrix(0, ncol=length(pred)+1, nrow=5, dimnames=list(NULL, c(pred, "aic"))) for (i in 1:5) { if (i==1) { bs <- dredge(full_mod, rank="AIC") } else { bs <- dredge(full_mod, rank="AIC", m.lim=c(i,i)) } bs_var <- attr(get.models(bs, 1)[[1]]$terms, "term.labels") for (j in 1:(ncol(sel_var)-1)) {sel_var[i,j] <- ifelse(names(sel_var[i,j]) %in% bs_var, 1, 0)} formula <- paste("sev.or.inte~", paste(names(sel_var[i,][sel_var[i,]==1]), collapse = "+")) sel_mod <- glm(formula, data = dat1a, family = binomial, x = T, y = T) sel_var[i, ncol(sel_var)] <- AIC(sel_mod) } # Report results -------------------------------------------------- out1 <- as.data.frame(sel_var) %>% mutate(AIC = round(aic,1)) %>% select(-aic) for (i in 1:(ncol(out1)-1)) {out1[,i] <- ifelse(out1[,i]==0, NA, "+")} out2 <- as.data.frame(t(out1)) colnames(out2) <- c("Best of all combinations", "Best combination of 2 variables", "Best combination of 3 variables", "Best combination of 4 variables", "Best combination of 5 variables") rownames(out2) <- c("- VCAM-1", "- SDC-1", "- Ang-2", "- IL-8", "- IP-10", "- IL-1RA", "- sCD163", "- sTREM-1", "- Ferritin", "- CRP", "- Viremia", "AIC of the selected model") knitr::kable(out2) ``` ########################################################################################### #### Appendix 9—table 1. Results of variable selection for children. ```{r table A9.1, eval=FALSE, include=FALSE} # Full model pred <- c("VCAM", "SDC", "Ang", "IL8", "IP10", "IL1RA", "CD163", "TREM", "Fer", "CRP") formula <- paste("sev.or.inte ~", paste(pred, collapse = "+")) full_mod <- glm(formula, data = dat1c, family = binomial, x = T, y = T) # Null model mg <- glm(sev.or.inte ~ 1, data=dat1c, x=T, y=T, family="binomial") # Backward elimination s1 <- step(glm(formula, data = dat1c, family = binomial, x=T,y=T), direction = "backward", scope = list(upper = formula), trace = 0) summary(s1) # Forward selection s2 <- step(glm(mg, data = dat1c, family = binomial, x=T,y=T), direction = "forward", scope = list(upper = formula, lower = mg), trace = 0) summary(s2) # Stepwise forward s3 <- step(glm(mg, data = dat1c, family = binomial, x=T,y=T), direction = "both", scope = list(upper = formula, lower = mg), trace = 0) summary(s3) # Stepwise backward s4 <- step(glm(formula, data = dat1c, family = binomial, x=T,y=T), direction = "both", scope = list(upper = formula, lower = mg), trace = 0) summary(s4) # Augmented backward elimination require(abe) s5 <- abe(glm(formula, data=dat1c, family=binomial, x=T, y=T), criterion="AIC", tau=0.1, data=dat1c) summary(s5) # Bayesian projection ## Thanks Vehtari A. et al. for the codes (https://avehtari.github.io/modelselection/diabetes.html) ## The codes took around 10 minutes to complete require(projpred) require(rstanarm) p <- 11 n <- dim(dat)[1] p0 <- 2 # prior guess for the number of relevant variables tau0 <- p0/(p-p0) * 1/sqrt(n) hs_prior <- hs(df=1, global_df=1, global_scale=tau0) # alternative horseshoe prior on weights t_prior <- student_t(df = 9, location = 0, scale = 2.5) post2 <- stan_glm(sev.or.inte ~ VCAM + SDC + Ang + IL8 + IP10 + IL1RA + CD163 + TREM + Fer + CRP, data = dat1c, family = binomial(link = "logit"), prior = hs_prior, prior_intercept = t_prior, seed = 51086, adapt_delta = 0.9) varsel2 <- cv_varsel(post2, method='forward', nloo = n) plot(varsel2, stats = c('elpd', 'pctcorr')) (nv <- suggest_size(varsel2, alpha=0.2)) proj2 <- project(varsel2, nv = nv, ns = 4000) round(colMeans(as.matrix(proj2)),1) round(posterior_interval(as.matrix(proj2)),1) ``` ########################################################################################### #### Appendix 9—table 2. Results of variable selection for adults. ```{r table A9.2, eval=FALSE, include=FALSE} # Full model pred <- c("VCAM", "SDC", "Ang", "IL8", "ns1(IP10)", "IL1RA", "CD163", "TREM", "Fer", "CRP") formula <- paste("sev.or.inte ~", paste(pred, collapse = "+")) full_mod <- glm(formula, data = dat1a, family = binomial, x = T, y = T) # Null model mg <- glm(sev.or.inte ~ 1, data=dat1a, x=T, y=T, family="binomial") # Backward elimination s1 <- step(glm(formula, data = dat1a, family = binomial, x=T,y=T), direction = "backward", scope = list(upper = formula), trace = 0) summary(s1) # Forward selection s2 <- step(glm(mg, data = dat1a, family = binomial, x=T,y=T), direction = "forward", scope = list(upper = formula, lower = mg), trace = 0) summary(s2) # Stepwise forward s3 <- step(glm(mg, data = dat1a, family = binomial, x=T,y=T), direction = "both", scope = list(upper = formula, lower = mg), trace = 0) summary(s3) # Stepwise backward s4 <- step(glm(formula, data = dat1a, family = binomial, x=T,y=T), direction = "both", scope = list(upper = formula, lower = mg), trace = 0) summary(s4) # Augmented backward elimination require(abe) s5 <- abe(glm(formula, data=dat1a, family=binomial, x=T, y=T), criterion="AIC", tau=0.1, data=dat1a) summary(s5) # Bayesian projection ## Thanks Vehtari A. et al. for the codes (https://avehtari.github.io/modelselection/diabetes.html) ## The codes took around 15 minutes to complete require(projpred) require(rstanarm) p <- 11 n <- dim(dat)[1] p0 <- 2 # prior guess for the number of relevant variables tau0 <- p0/(p-p0) * 1/sqrt(n) hs_prior <- hs(df=1, global_df=1, global_scale=tau0) # alternative horseshoe prior on weights t_prior <- student_t(df = 10, location = 0, scale = 2.5) require(MASS) post2 <- stan_gamm4(sev.or.inte ~ VCAM + SDC + Ang + IL8 + s(IP10) + IL1RA + CD163 + TREM + Fer + CRP, data = dat1a, family = binomial(link = "logit"), prior = hs_prior, prior_intercept = t_prior, seed = 51086, adapt_delta = 0.9) # Projection predictive variable selection varsel2 <- cv_varsel(post2, method='forward', nloo = n) plot(varsel2, stats = c('elpd', 'pctcorr')) (nv <- suggest_size(varsel2, alpha=0.2)) proj2b <- project(varsel2, nv = nv, ns = 4000) round(colMeans(as.matrix(proj2b)),1) round(posterior_interval(as.matrix(proj2b)),1) ```