swh:1:snp:531311172177ecad060ca11b9c3752edb33ce261
Tip revision: 847d8e0f564eeb3f075b443205fb3384598bc2b4 authored by Nguyen Lam Vuong on 30 July 2021, 23:46:27 UTC
Delete .DS_Store
Delete .DS_Store
Tip revision: 847d8e0
Elife ERA codes.Rmd
---
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)
```
