https://medium.com/@ericpgreen/my-wife-couldnt-sleep-550e206cd80a

ICP.R

# Plot code for https://medium.com/@ericpgreen/my-wife-couldnt-sleep-550e206cd80a
# Eric Green, @ericpgreen

library(tidyverse)
library(ggrepel)
library(lubridate)

# Puljic et al. (2015). The risk of infant and fetal death by each additional week of expectant management in intrahepatic cholestasis of pregnancy by gestational age. American Journal of Obstetrics and Gynecology, 212(5), 667-e1.
# create dataframe from published table
puljic <- data.frame(ga=c(rep(seq(from=34, to=40, by=1), 2)),
                     group=c(rep("ICP", 7),
                             rep("Control", 7)),
                     sb=c(2.3, 4.4, 6.8, 8, 4.7, 11.1, 26.5,
                          1.7, 1.9, 2.1, 2.3, 3.2, 4.2, 5.8),
                     sbl=c(0, 0, 0, 0, 0, 0, 0, 
                           1.5, 1.7, 1.9, 2.1, 2.9, 3.8, 5.2),
                     sbu=c(6.2, 9.9, 13.8, 16, 11.9, 25.1, 56.5,
                           1.9, 2.1, 2.3, 2.5, 3.5, 4.5, 6.4),
                     id=c(22.2, 26.9, 4.7, 12.3, 13.7, 18.3, 22.5,
                          42.1, 27.1, 22.9, 18, 11.8, 9.8, 10.4),
                     idl=c(10.2, 13.5, 0, 2.4, 1.5, .5, 0, 
                           41.1, 26.4, 22.2, 17.3, 11.3, 9.3, 9.8),
                     idu=c(34.2, 40.3, 10.5, 22.3, 26.0, 36.2, 50.2,
                           43, 27.9, 23.6, 18.6, 12.3, 10.3, 11)
                     )

# add a label for plotting with ggrepel
puljic <- puljic %>%
  mutate(label = if_else(ga == 40, as.character(group), NA_character_))

# create another dataframe from published table
puljic2 <- data.frame(ga=c(rep(seq(from=34, to=40, by=1), 2)),
                      group=c(rep("Expectant Management", 7),
                              rep("Delivery", 7)),
                      id=c(29.2, 9.1, 19.2, 21.7, 23.1, 33.6, 25.18,
                           22.2, 26.9, 4.7, 12.3, 13.7, 18.3, 22.5),
                      idl=c(15.5, 1.4, 7.6, 8.5, 7.2, 9.5, 0,
                            10.2, 13.5, 0, 2.4, 1.5, .5, 0),
                      idu=c(43, 16.9, 30.8, 35, 38.9, 57.8, 54.51,
                            34.2, 40.3, 10.5, 22.3, 26, 36.2, 50.2)
)

# add a label for plotting with ggrepel
puljic2 <- puljic2 %>%
  mutate(label = if_else(ga == 40, as.character(group), NA_character_))

# plot stillbirth
ggplot(puljic, aes(x=factor(ga), y=sb, group=group, color=group)) +
  geom_line(aes(linetype=group), size=2) +
  geom_ribbon(aes(ymin=puljic$sbl, 
                  ymax=puljic$sbu,
                  fill=puljic$group), linetype="dotted", alpha=0.1) +
  geom_label_repel(aes(label=label, fill=label),
                   nudge_y = -2.5, 
                   nudge_x = 1,
                   size=5, color='white',
                   force=1.5, segment.color='#bbbbbb',
                   na.rm = TRUE) +
  scale_color_manual(values=c("#53357a", "#edbe40")) +
  scale_fill_manual(values=c("#53357a", "#edbe40")) +
  geom_point(size=4, color="#7f7e7a", fill="white", shape=21, stroke = 2) +
  theme_light() + 
  ggtitle("Risk of stillbirth in women with and without ICP") +
  xlab("Weeks gestation") +
  ylab("Stillbirths (per 10,000 ongoing pregnancies)") +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold"))

# plot infant death
ggplot(puljic, aes(x=factor(ga), y=id, group=group, color=group)) +
  geom_line(aes(linetype=group), size=2) +
  geom_ribbon(aes(ymin=puljic$idl, 
                  ymax=puljic$idu,
                  fill=puljic$group), linetype="dotted", alpha=0.1) +
  geom_label_repel(aes(label=label, fill=label),
                   nudge_y = -2.5, 
                   nudge_x = 1,
                   size=5, color='white',
                   force=1.5, segment.color='#bbbbbb',
                   na.rm = TRUE) +
  scale_color_manual(values=c("#53357a", "#edbe40")) +
  scale_fill_manual(values=c("#53357a", "#edbe40")) +
  geom_point(size=4, color="#7f7e7a", fill="white", shape=21, stroke = 2) +
  theme_light() + 
  ggtitle("Risk of infant death in women with and without ICP") +
  xlab("Weeks gestation") +
  ylab("Infant deaths (per 10,000 live births)") +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold"))

# plot delivery vs expectant management
ggplot(puljic2, aes(x=factor(ga), y=id, group=group, color=group)) +
  geom_line(aes(linetype=group), size=2) +
  geom_ribbon(aes(ymin=puljic2$idl, 
                  ymax=puljic2$idu,
                  fill=puljic2$group), linetype="dotted", alpha=0.1) +
  geom_label_repel(aes(label=label, fill=label),
                   nudge_y = -2.5, 
                   nudge_x = 1,
                   size=5, color='white',
                   force=1.5, segment.color='#bbbbbb',
                   na.rm = TRUE) +
  scale_color_manual(values=c("#53357a", "#edbe40")) +
  scale_fill_manual(values=c("#53357a", "#edbe40")) +
  geom_point(size=4, color="#7f7e7a", fill="white", shape=21, stroke = 2) +
  theme_light() + 
  ggtitle("Risk of delivery vs expectant management among women with ICP") +
  xlab("Weeks gestation") +
  ylab("Deaths (per 10,000)") +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold"))


(1) delivery 48h after steroid administration; 
(2) amniocentesis with steroid administration if fetal lung maturity is negative with delivery 48h after steroids; 
(3) amniocentesis with re-testing in one week if fetal lung maturity is negative or 
(4) immediate delivery

# Lo, J. O., Shaffer, B. L., Allen, A. J., Little, S. E., Cheng, Y. W., & Caughey, A. B. (2015). Intrahepatic cholestasis of pregnancy and timing of delivery. The Journal of Maternal-Fetal & Neonatal Medicine, 28(18), 2254-2258.
# create dataframe from published table
lo <- data.frame(ga=c(rep(35, 4), rep(36, 4), rep(37, 4), rep(38, 4)),
                 amino=c(rep(c(0, 1, 1, 0), 4)),
                 steriods=c(rep(c(1, 1, 0, 0), 4)),
                 strategy=c("Delivery 48h after steriods",
                            "Amniocentesis and delivery 48h after steriods",
                            "Amniocentesis with re-testing, no steriods",
                            "Immediate delivery"),
                 qaly=c(565381,
                        565403,
                        565590,
                        565640,
                        565616,
                        565330,
                        565833,
                        565879,
                        564937,
                        564607,
                        565159,
                        565202,
                        564820,
                        564192,
                        564597,
                        564862),
                 iufd=c(37,
                        44,
                        31,
                        29,
                        66,
                        77,
                        60,
                        58,
                        95,
                        109,
                        88,
                        87,
                        87,
                        139,
                        124,
                        116),
                 rds=c(430,
                       111,
                       491,
                       638,
                       219,
                       53,
                       244,
                       328,
                       26,
                       8,
                       29,
                       40,
                       28,
                       11,
                       26,
                       40),
                 cp=c(50,
                      46,
                      50,
                      50,
                      31,
                      29,
                      31,
                      31,
                      25,
                      23,
                      25,
                      25,
                      13,
                      12,
                      13,
                      13))

# construct difference in qaly, add label for plotting with ggrepel
loTidy <- lo %>%
  mutate(qaly2 = qaly-min(qaly)) %>%
  gather(key = "outcome", value = "value", qaly:qaly2) %>%
  mutate(label = if_else(ga == 40, as.character(strategy), NA_character_))

# plot outcomes
ggplot(loTidy %>% 
         filter(outcome!="qaly" & outcome!="qaly2") %>%
         mutate(outcome = recode(outcome,
                                 "cp" = "Cerebral palsy",
                                 "iufd" = "Stillbirth",
                                 "rds" = "Respiratory Distress Syndrome")), 
       aes(x=factor(ga), y=value, fill=strategy)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  #scale_fill_brewer(palette = "Dark2") +
  scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9", "#53357a")) +
  theme_light() + 
  ggtitle("Simulated outcomes for 10,000 pregnant women with ICP") +
  xlab("Weeks gestation") +
  ylab("Events") +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold")) +
  facet_wrap(~outcome, scales="free_y", ncol=3) +
  theme(legend.position="bottom") +
  guides(fill = guide_legend(title="",
                             nrow=2))

# plot qaly
ggplot(loTidy %>% 
         filter(outcome=="qaly2") %>%
         mutate(outcome = recode(outcome,
                                 "qaly2" = "Additional QALYs (relative to min)")), 
       aes(x=factor(ga), y=value, fill=strategy)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  #scale_fill_brewer(palette = "Dark2") +
  scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9", "#53357a")) +
  theme_light() + 
  ggtitle("Simulated additional QALYs for 10,000 pregnant women with ICP") +
  xlab("Weeks gestation") +
  ylab("QALYs") +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold")) +
  theme(legend.position="bottom") +
  guides(fill = guide_legend(title="",
                             nrow=2))

# Eve's bile acid results
eve <- data.frame(date = ymd(c("2018-05-24", 
                               "2018-06-07",
                               "2018-06-14",
                               "2018-06-18", 
                               "2018-06-25",
                               "2018-06-28",
                               "2018-07-02",
                               "2018-07-05"
                              )),
                  ba = c(11.5, # "normal", 3.7-14.5 ref for preg female 15-45
                         25.6,
                         85.6,
                         57.1, 
                         76.4,
                         55.2,
                         67.1,
                         38.4))

ggplot(eve, aes(x = date, y = ba)) + 
  geom_hline(yintercept=10, linetype="dashed", color = "#53357a", size=1) +
  geom_hline(yintercept=20, linetype="dashed", color = "#53357a", size=1) +
  geom_hline(yintercept=40, linetype="dashed", color = "#53357a", size=1) +
  geom_hline(yintercept=100, linetype="dashed", color = "#53357a", size=1) +
  geom_vline(xintercept=ymd("2018-07-01"), #35
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-06-24"), #34
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-06-17"), #33
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-06-10"), #32
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-06-03"), #31
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-05-27"), #30
             linetype="solid", 
             color = "grey", size=0.5, alpha=0.5) +
  geom_vline(xintercept=ymd("2018-06-08"), 
             linetype="dotted", 
             color = "#db3458", size=0.5, alpha=0.5) +
  geom_line(size=2, color="#edbe40") +
  geom_line(size=1, color="#edbe40") +
  geom_point(size=4, color="#7f7e7a", fill="white", shape=21, stroke = 2) +
  theme_light() + 
  ylim(0, 120) +
  annotate(geom="text", x=ymd("2018-06-15"), y=15, 
           label="ICP",
           color="#53357a",
           fontface = 2) +
  annotate(geom="text", x=ymd("2018-06-15"), y=25, 
           label="Moderate ICP",
           color="#53357a",
           fontface = 2) +
  annotate(geom="text", x=ymd("2018-06-15"), y=45, 
           label="Severe ICP",
           color="#53357a",
           fontface = 2) +
  annotate(geom="text", x=ymd("2018-06-15"), y=105, 
           label="Marker of elevated risk of stillbirth",
           color="#53357a",
           fontface = 2) +
  annotate(geom="text", x=ymd("2018-05-27"), y=120, 
           label="w30",
           color="grey") +
  annotate(geom="text", x=ymd("2018-06-03"), y=120, 
           label="w31",
           color="grey") +
  annotate(geom="text", x=ymd("2018-06-10"), y=120, 
           label="w32",
           color="grey") +
  annotate(geom="text", x=ymd("2018-06-17"), y=120, 
           label="w33",
           color="grey") +
  annotate(geom="text", x=ymd("2018-06-24"), y=120, 
           label="w34",
           color="grey") +
  annotate(geom="text", x=ymd("2018-07-01"), y=120, 
           label="w35",
           color="grey") +
  annotate(geom="text", x=ymd("2018-06-11"), y=95, 
           label="Started UDCA",
           color="#db3458") +
  annotate(geom="text", x=ymd("2018-06-06"), y=29, 
           label="Diagnosed with ICP",
           color="#7f7e7a",
           hjust = 1) +
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        title=element_text(size=16,face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  ggtitle("Eve's bile acid levels over time") +
  xlab("Date") +
  ylab("Bile acid level (µmol/L)")