Comparing Exposures with back to back Density/Boxplots

In this post, I will cover visualizing simulated drug exposures head to head compare using some ggplot2 tricks. First we use the NHANES body weight (kg), illustrate the effects of taking into account survey weights in the quantile regression showing the 5^th, 50^th, and 95^th percentiles. I also overlay the simulated demographics (red points).

Code

library(ggplot2)
library(tidyr)
library(dplyr)
library(patchwork)
library(nhanesgamlss)
library(gganimate)

factor.with.units <- function(x, units) {
  l <- sort(unique(x))
  y <- factor(x, levels=l, labels=paste(l, units))
}


# note https://cran.r-project.org/web/packages/RNHANES/vignettes/introduction.html
# RNHANES help in applying survey weights and in downloading NHANES data
# I am using a bmxdata.csv
# bmxdata <- nhanes_load_data("BMX", "2013-2014", cache = "./nhanes_data",
#                              demographics = TRUE,recode = TRUE)

bmxdata <- read.csv("bmxdata.csv")
pedcov4000 <- read.csv("pedcov4000.csv")  %>% 
         filter(AGEY<=18) %>% 
      mutate(RIAGENDR = ifelse(SEX=="boys","Male","Female"))
ggplot(bmxdata %>% 
         filter(RIDAGEYR<=18),aes(RIDAGEYR,BMXWT))+
  geom_point(aes(alpha=WTMEC2YR/1000) )+
  geom_quantile(method = "rqss",quantiles = c(0.05,0.5,0.95),lambda= 10,
  aes(weight = WTMEC2YR,col="rqss-weighted"),linewidth=1.5, alpha = 0.6)+
  geom_quantile(method = "rqss",quantiles = c(0.05,0.5,0.95),lambda= 10,
  aes(col="rqss"),linewidth=1.5, alpha = 0.6)+
  geom_point(data = pedcov4000,aes(x=AGEY,y = WT),col="red",alpha=0.1)+
  facet_grid(~RIAGENDR)+
  labs(x="Years", y ="Weight (kg)",color="quantile regression",
       alpha="NHANES Weights")+
    scale_color_manual(values=c("#093B6D", "#EF761B"))+
    theme_bw()+
    theme(strip.placement = "outside",
          strip.text.y.left = element_text(angle=0),
          strip.background = ggplot2::element_rect(fill = "#475c6b"), 
          strip.text = ggplot2::element_text(face = "bold", 
                                             color = "white"))

I then plug these demographics into an mrgsolve model and simulate flat fixed dose of 150 mg versus a 2 mg/kg dose ( 150 mg for a 75 kg individual). The resulting areas under the curves (AUC’s) are compared using boxplots.

Code

library(mrgsolve)

codepedcov <- '
$PARAM
KA=0.365, CL = 2, Vc=22,  Vp1 =200  ,Qp1= 4
WT=75, AGE=18, SEX = 1
$CMT GUT CENT PER1
$MAIN
double KAi = KA;
double CLi = CL *pow((WT/75), 0.75)*exp(ETA(1));
double Vci = Vc *pow((WT/75),    1)*exp(ETA(2));
double Vp1i = Vp1 *pow((WT/75),    1);
double Qp1i = Qp1 *pow((WT/75), 0.75);
double Keli  = CLi/Vci    ;
double Kpt1i =Vp1i/Vci    ;
double Ktp1i =Qp1i/Vp1i   ;
$OMEGA
0.6
$OMEGA
0.6

$ODE
dxdt_GUT  = -KAi*GUT;
dxdt_CENT =  KAi*GUT-Kpt1i*CENT-Keli*CENT+Ktp1i*PER1;
dxdt_PER1 =          Kpt1i*CENT          -Ktp1i*PER1;

$TABLE
double CP    = CENT/ Vci;
double CPER1 = PER1/Vp1i;
$CAPTURE CP CPER1 CLi
'
modcovpedsim <- mcode("codepedcov", codepedcov)
idata <- data_frame(ID=1:nrow(pedcov4000),
                    WT=pedcov4000$WT,
                    AGE=pedcov4000$AGEY,
                    SEX=as.double(as.factor(pedcov4000$SEX )),
                    TRT =1 
)
idata2 <- data_frame(ID=1:nrow(pedcov4000),
                     WT=pedcov4000$WT,
                     AGE=pedcov4000$AGEY,
                     SEX=as.double(as.factor(pedcov4000$SEX )),
                     TRT =2 
)
idata <- rbind(idata,idata2)
rm(idata2)
ev1 <- ev(time=0, amt=150, cmt=1)
data.dose <- ev(ev1)
data.dose<-as.data.frame(data.dose)
data.all<-merge(idata,data.dose)
data.all$amt<- ifelse( data.all$TRT==1,150,2*round(data.all$WT,0))

outpedsim<-modcovpedsim %>%
  data_set(data.all) %>%
  carry.out(TRT,WT,AGE,SEX,CLi) %>%
  mrgsim(end=24, delta=1)
outpedsim<-as.data.frame(outpedsim)
outpedsim <- outpedsim %>% 
  arrange(ID,time,WT)

out.ped.nca <- outpedsim %>% 
  group_by(ID,TRT,SEX,WT,AGE)%>% 
  summarise (Cmax = max(CP,na.rm = TRUE),
             Clast= CP[n()],
             AUC= sum(diff(time ) *na.omit(lead(CP) + CP)) / 2,
             CLi= median(CLi)) 

out.ped.nca$TRTC <- ifelse(out.ped.nca$TRT==1,"fixed Dose","perkg Dose")
out.ped.nca$WTC <- cut(out.ped.nca$WT,breaks = c(15,30,45,60,125))
out.ped.nca$`Weight Category`  <- factor.with.units(out.ped.nca$WTC, "{kg}")

out.ped.nca$REGIMEN <-  "a.BASE"
ggplot(out.ped.nca, aes(x=`Weight Category`, AUC))+
  annotate(geom="rect",xmin=-Inf,xmax=Inf,ymin=0.25,ymax=1.5,alpha=0.1,fill="blue")+
  geom_boxplot(aes(fill=TRTC,color = TRTC))+
  facet_grid(REGIMEN~`Weight Category`,
             labeller = labeller(`Weight Category`= label_value,
                                 REGIMEN = label_value),scales="free")+
  theme_bw()+
  theme(axis.text.x = element_blank(),legend.position="bottom",
        axis.ticks.x=element_blank(),
        strip.text=element_text(size=12),
        plot.title = element_text(size=14),
        axis.text.y=element_text(size=12),axis.title.y = element_text(size=16))+
  guides(color=guide_legend(reverse=TRUE),fill=guide_legend(reverse=TRUE))+
  theme(strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", 
                                           color = "white"))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.5)) +
  labs(y="Area Under the Curve (mg*h/mL)",x="")+
  coord_trans(y = "log")+
  guides(color=guide_legend(reverse=TRUE),
         fill=guide_legend(reverse=TRUE))

Ideally we want to keep what we compare close together i.e. the fixed Dose exposures next to the perKg Dose as shown in the second plot above. As a pharmacometrician, I like to see the distribution not just the boxplot so I come up with a hybrid density/boxplot visual below which I will call the back to back density boxplot.

Code

statsdata <- out.ped.nca %>%
  group_by(TRTC ,WTC,`Weight Category` ) %>%
  dplyr::summarize(
    low=quantile(AUC, probs = c(0.05) ),
    quart=quantile(AUC, probs = c(0.25) ),
    med=quantile(AUC, probs = c(0.5) ),
    seven=quantile(AUC, probs = c(0.75) ),
    up=quantile(AUC, probs = c(0.95) )
  )

statsdata$REGIMEN <-"a.BASE"

backtobackhorizontal <- ggplot(out.ped.nca, aes(x = AUC,fill=TRTC ,col=TRTC ),alpha=0.1) +
facet_grid(REGIMEN~`Weight Category`,
           labeller = labeller(`Weight Category`= label_value,
                               REGIMEN = label_value),scales="free_y")+
        annotate(geom="rect",
             ymin=-Inf,ymax=Inf,xmin=0.25,xmax=1.5,alpha=0.1,fill="blue")+
  stat_density(data=out.ped.nca[out.ped.nca$TRT ==1,],aes(ymax = 0.5*..scaled..,  ymin = rep(0)),
               colour = "transparent",geom = "ribbon", position = "identity",alpha=0.1) +
  stat_density(data=out.ped.nca[out.ped.nca$TRT ==2,],aes(ymin = -0.5*..scaled..,  ymax = rep(0)),
               colour = "transparent",geom = "ribbon", position = "identity",alpha=0.1) +
  geom_point(data=statsdata[statsdata$TRTC=="perkg Dose",]  , aes(x=med),y=-0.1)+
  geom_point(data=statsdata[statsdata$TRTC!="perkg Dose",]  , aes(x=med),y=0.1)+
  geom_errorbarh(data=statsdata[statsdata$TRTC =="perkg Dose",]  ,
                 aes(x=med,y=-0.1,xmax = up, xmin = low),width = 0.1)+
  geom_errorbarh(data=statsdata[statsdata$TRTC !="perkg Dose",]  ,
                 aes(x=med,y=0.1,xmax = up, xmin = low),width = 0.1)+
  geom_errorbarh(data=statsdata[statsdata$TRTC =="perkg Dose",]  ,
                 aes(x=med,y=-0.1,xmax = seven, xmin = quart),width = 0.15)+
  geom_errorbarh(data=statsdata[statsdata$TRTC !="perkg Dose",]  ,
                 aes(x=med,y=0.1,xmax = seven, xmin =quart),width = 0.15)+
  labs(color="",fill="")

backtobackhorizontal+
  coord_flip()+
  labs(col="Dosing Type",fill="Dosing Type",
       title="Comparing AUC Disributions:\nFixed versus a Weight-Based Dosing Regimens")+
  scale_x_continuous("Area Under the Curve (mg*h/mL)") +
  ylab("")+ 
  theme_bw()+
  theme(axis.text.x = element_blank(),legend.position="bottom",
        axis.ticks.x=element_blank(),
        strip.text=element_text(size=12),
        plot.title = element_text(size=14),
        axis.text.y=element_text(size=12),axis.title.y = element_text(size=16))+
  guides(color=guide_legend(reverse=TRUE),fill=guide_legend(reverse=TRUE))+
  theme(strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", 
                                           color = "white"))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.2))

We see that the fixed dosing has the densities crossing above the therapeutic range [0.25-1.5] mg*h/mL.
Next, we simulate two fixed dosing strategy reducing the fixed doses by a percentage:

Dose Strategy b: dose reduced by 80%, 50% and 40 % for the (45,60], (30,45], and (15,30] kg, respectively.
Dose Strategy c: dose reduced by 75%, 50% and 30 % for the (45,60], (30,45], and (15,30] kg, respectively.

We compare the outputs using boxplots:

Code

# simulate dose strategy 1 80 % 50 % 40 %

ev1 <- ev(time=0,amt=150, cmt=1)
data.dose <- ev(ev1)
data.dose<-as.data.frame(data.dose)
data.all<-merge(idata,data.dose)
data.all$amt<- ifelse( data.all$TRT==1,150,2*round(data.all$WT,0))
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<60 ,150*0.8,data.all$amt)
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<45 ,150*0.5,data.all$amt)
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<30 ,150*0.4,data.all$amt)

outpedsim<-modcovpedsim %>%
  data_set(data.all) %>%
  carry.out(TRT,WT,AGE,SEX,CLi) %>%
  mrgsim(end=24, delta=1)
outpedsim<-as.data.frame(outpedsim)
outpedsim <- outpedsim %>% 
  arrange(ID,time,WT)

out.ped.nca.1 <- outpedsim %>% 
  group_by(ID,TRT,SEX,WT,AGE)%>% 
  summarise (Cmax = max(CP,na.rm = TRUE),
             Clast= CP[n()],
             AUC= sum(diff(time ) *na.omit(lead(CP) + CP)) / 2,
             CLi= median(CLi)) 
out.ped.nca.1$TRTC <- ifelse(out.ped.nca.1$TRT==1,"fixed Dose","perkg Dose")
out.ped.nca.1$WTC <- cut(out.ped.nca.1$WT,breaks = c(15,30,45,60,125))
out.ped.nca.1$`Weight Category`  <- factor.with.units(out.ped.nca.1$WTC, "{kg}")


# simulate dose strategy 2 75 % 50 % 30 %

ev1 <- ev(time=0,amt=150, cmt=1)
data.dose <- ev(ev1)
data.dose<-as.data.frame(data.dose)
data.all<-merge(idata,data.dose)
data.all$amt<- ifelse( data.all$TRT==1,150,2*round(data.all$WT,0))
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<60 ,150*0.75,data.all$amt)
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<45 ,150*0.5,data.all$amt)
data.all$amt<- ifelse( data.all$TRT==1&data.all$WT<30 ,150*0.30,data.all$amt)


outpedsim<-modcovpedsim %>%
  data_set(data.all) %>%
  carry.out(TRT,WT,AGE,SEX,CLi) %>%
  mrgsim(end=24, delta=1)
outpedsim<-as.data.frame(outpedsim)
outpedsim <- outpedsim %>% 
  arrange(ID,time,WT)

out.ped.nca.2 <- outpedsim %>% 
  group_by(ID,TRT,SEX,WT,AGE)%>% 
  summarise (Cmax = max(CP,na.rm = TRUE),
             Clast= CP[n()],
             AUC= sum(diff(time ) *na.omit(lead(CP) + CP)) / 2,
             CLi= median(CLi)) 
out.ped.nca.2$TRTC <- ifelse(out.ped.nca.2$TRT==1,"fixed Dose","perkg Dose")
out.ped.nca.2$WTC <- cut(out.ped.nca.2$WT,breaks = c(15,30,45,60,125))
out.ped.nca.2$`Weight Category`  <- factor.with.units(out.ped.nca.2$WTC, "{kg}")

statsdata1 <- out.ped.nca.1 %>%
  group_by(TRTC ,WTC,`Weight Category` ) %>%
  dplyr::summarize(
    low=quantile(AUC, probs = c(0.05) ),
    quart=quantile(AUC, probs = c(0.25) ),
    med=quantile(AUC, probs = c(0.5) ),
    seven=quantile(AUC, probs = c(0.75) ),
    up=quantile(AUC, probs = c(0.95) )
  )

statsdata2 <- out.ped.nca.2 %>%
  group_by(TRTC ,WTC,`Weight Category` ) %>%
  dplyr::summarize(
    low=quantile(AUC, probs = c(0.05) ),
    quart=quantile(AUC, probs = c(0.25) ),
    med=quantile(AUC, probs = c(0.5) ),
    seven=quantile(AUC, probs = c(0.75) ),
    up=quantile(AUC, probs = c(0.95) )
  )
statsdata1$REGIMEN <-"b.OPTIMIZED1"
statsdata2$REGIMEN <-"c.OPTIMIZED2"
statsdata1all <-rbind(statsdata,statsdata1,statsdata2)

out.ped.nca$REGIMEN <-  "a.BASE"
out.ped.nca.1$REGIMEN <-"b.OPTIMIZED1"
out.ped.nca.2$REGIMEN <-"c.OPTIMIZED2"

out.ped.nca$WTX <- ifelse(out.ped.nca$`Weight Category`=="(15,30] {kg}",(15+30)/2,NA)
out.ped.nca$WTX <- ifelse(out.ped.nca$`Weight Category`=="(30,45] {kg}",(45+30)/2,out.ped.nca$WTX)
out.ped.nca$WTX <- ifelse(out.ped.nca$`Weight Category`=="(45,60] {kg}",(45+60)/2,out.ped.nca$WTX)
out.ped.nca$WTX <- ifelse(out.ped.nca$`Weight Category`=="(60,125] {kg}",(60+125)/2,out.ped.nca$WTX)

out.ped.nca.2$WTX <- ifelse(out.ped.nca.2$`Weight Category`=="(15,30] {kg}",(15+30)/2,NA)
out.ped.nca.2$WTX <- ifelse(out.ped.nca.2$`Weight Category`=="(30,45] {kg}",(45+30)/2,out.ped.nca.2$WTX)
out.ped.nca.2$WTX <- ifelse(out.ped.nca.2$`Weight Category`=="(45,60] {kg}",(45+60)/2,out.ped.nca.2$WTX)
out.ped.nca.2$WTX <- ifelse(out.ped.nca.2$`Weight Category`=="(60,125] {kg}",(60+125)/2,out.ped.nca.2$WTX)
out.ped.nca.1$WTX <- ifelse(out.ped.nca.1$`Weight Category`=="(15,30] {kg}",(15+30)/2,NA)
out.ped.nca.1$WTX <- ifelse(out.ped.nca.1$`Weight Category`=="(30,45] {kg}",(45+30)/2,out.ped.nca.1$WTX)
out.ped.nca.1$WTX <- ifelse(out.ped.nca.1$`Weight Category`=="(45,60] {kg}",(45+60)/2,out.ped.nca.1$WTX)
out.ped.nca.1$WTX <- ifelse(out.ped.nca.1$`Weight Category`=="(60,125] {kg}",(60+125)/2,out.ped.nca.1$WTX)


out.ped.nca.all <-rbind(out.ped.nca,out.ped.nca.1,out.ped.nca.2)

ggplot(out.ped.nca.all, aes(x=`Weight Category`, AUC))+
  annotate(geom="rect",xmin=-Inf,xmax=Inf,ymin=0.25,ymax=1.5,alpha=0.1,fill="blue")+
  geom_boxplot(aes(fill=TRTC,color = TRTC))+
  facet_grid(REGIMEN~`Weight Category`,
             labeller = labeller(`Weight Category`= label_value,
                                 REGIMEN = label_value),scales="free")+
  theme_bw()+
  theme(axis.text.x = element_blank(),legend.position="bottom",
        axis.ticks.x=element_blank(),
        strip.text=element_text(size=12),
        plot.title = element_text(size=14),
        axis.text.y=element_text(size=12),axis.title.y = element_text(size=16))+
  guides(color=guide_legend(reverse=TRUE),fill=guide_legend(reverse=TRUE))+
  theme(strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", 
                                           color = "white"))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.5)) +
  labs(y="Area Under the Curve (mg*h/mL)",x="")+
  coord_trans(y = "log")+
  guides(color=guide_legend(reverse=TRUE),
         fill=guide_legend(reverse=TRUE))

We also re-do the back to back density boxplot comparing all three regimens:

Code

backtobackhorizontalall <- ggplot(out.ped.nca.all, aes(x = AUC,fill=TRTC ,col=TRTC ),alpha=0.1) +
  annotate(geom="rect",ymin=-Inf,ymax=Inf,xmin=0.25,xmax=1.5,alpha=0.1,fill="blue")+
  stat_density(data=out.ped.nca.all[out.ped.nca.all$TRT ==1,],aes(ymax = 0.5*..scaled..,  ymin = rep(0)),
               colour = "transparent",geom = "ribbon", position = "identity",alpha=0.1) +
  stat_density(data=out.ped.nca.all[out.ped.nca.all$TRT ==2,],aes(ymin = -0.5*..scaled..,  ymax = rep(0)),
               colour = "transparent",geom = "ribbon", position = "identity",alpha=0.1) +
  geom_point(data=statsdata1all[statsdata1all$TRTC=="perkg Dose",]  , aes(x=med),y=-0.1)+
  geom_point(data=statsdata1all[statsdata1all$TRTC!="perkg Dose",]  , aes(x=med),y=0.1)+
  geom_errorbarh(data=statsdata1all[statsdata1all$TRTC =="perkg Dose",]  ,
                 aes(x=med,y=-0.1,xmax = up, xmin = low),width = 0.1)+
  geom_errorbarh(data=statsdata1all[statsdata1all$TRTC !="perkg Dose",]  ,
                 aes(x=med,y=0.1,xmax = up, xmin = low),width = 0.1)+
  geom_errorbarh(data=statsdata1all[statsdata1all$TRTC =="perkg Dose",]  ,
                 aes(x=med,y=-0.1,xmax = seven, xmin = quart),width = 0.15)+
  geom_errorbarh(data=statsdata1all[statsdata1all$TRTC !="perkg Dose",]  ,
                 aes(x=med,y=0.1,xmax = seven, xmin =quart),width = 0.15)+
  labs(color="",fill="",y="")+
  coord_flip()+
  scale_x_continuous("Area Under the Curve (mU*h/mL)") +
  labs(col="Dosing Type",fill="Dosing Type",
       title="Comparing AUC Disributions:\nFixed versus a Weight-Based Dosing Regimens by Weight Category")+
  theme(axis.text.x = element_blank(),legend.position="bottom", axis.ticks.x=element_blank(),
        strip.text=element_text(size=12),
        plot.title = element_text(size=14),
        axis.text.y=element_text(size=12),axis.title.y = element_text(size=16))+
  facet_grid(REGIMEN~`Weight Category`,labeller = labeller(`Weight Category`= label_value,
                                                           REGIMEN = label_value
                                                           ),scales="free_y")+
  guides(color=guide_legend(reverse=TRUE),
         fill=guide_legend(reverse=TRUE))
backtobackhorizontalall+ 
  theme_bw()+
  theme(axis.text.x = element_blank(),legend.position="bottom",
        axis.ticks.x=element_blank(),
        strip.text=element_text(size=12),
        plot.title = element_text(size=14),
        axis.text.y=element_text(size=12),axis.title.y = element_text(size=16))+
  guides(color=guide_legend(reverse=TRUE),fill=guide_legend(reverse=TRUE))+
  theme(strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", 
                                           color = "white"))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.2))

As a bonus, I show how gganimate can be used to animate and transition the boxplots across regimens:

Code

labeldata<- rbind(
  data.frame(TRTC = c(rep("fixed Dose",4),rep("perkg Dose",4)),
             x=c((15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2,
                 (15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2), 
             y = rep(Inf,8),
             label= c("//","//","//","   150 mg   ",
                      "//","//","//","   2 mg/kg   "),
             REGIMEN = "a.BASE"),
  data.frame(TRTC = c(rep("fixed Dose",4),rep("perkg Dose",4)),
             x=c((15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2,
                 (15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2), 
             y = rep(Inf,8),
             label= c("x0.4","x0.5","x0.80","   150 mg   ",
                      "//","//","//","   2 mg/kg   "),
             REGIMEN = "b.OPTIMIZED1" ),
  data.frame(TRTC = c(rep("fixed Dose",4),rep("perkg Dose",4)),
             x=c((15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2,
                 (15+30)/2,(45+30)/2,(45+60)/2,(60+125)/2), 
             y = rep(Inf,8),
             label= c("x0.30","x0.5","x0.75","   150 mg   ",
                      "//","//","//","   2 mg/kg   "),
             REGIMEN = "c.OPTIMIZED2" )
)

ggplot(out.ped.nca.all, aes(WT, AUC,col = TRTC, fill = TRTC))+
  annotate(geom="rect",xmin=-Inf,xmax=Inf,ymin=0.25,ymax=1.5,alpha=0.1,fill="blue")+
  geom_vline(data=data.frame(xintercept= c(15,30,45,60) ),
             aes(xintercept=xintercept) )+
  geom_boxplot(aes(x=WTX,group=WTX),alpha=0.5)+
  geom_point(alpha=0.01,size=1)+
  geom_label(data=labeldata,aes(x,y,label=label),vjust="inward",fill="white",
             show.legend = FALSE)+
  facet_grid(REGIMEN~TRTC,labeller=label_value)+
  scale_x_continuous(breaks=c(15,30,45,60,125))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.2))+
  labs(x="Weight (kg)",y="Area Under the Curve (mg*h/mL)")+
  theme_bw()+
  theme(strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", 
                                           color = "white"),
        legend.position = "bottom")

Instead of facetting by regimen now we use transition_states(REGIMEN,...):

Code

ggplot(out.ped.nca.all, aes(WT, AUC,col = TRTC, fill = TRTC))+
  annotate(geom="rect",xmin=-Inf,xmax=Inf,ymin=0.25,ymax=1.5,alpha=0.1,fill="blue")+
  geom_vline(data=data.frame(xintercept= c(15,30,45,60) ),aes(xintercept=xintercept) )+
  geom_boxplot(aes(x=WTX,group=WTX),alpha=0.5)+
  geom_point(alpha=0.01,size=1)+
  geom_label(data=labeldata,aes(x,y,label=label),vjust="inward",fill="white")+
  facet_grid(~TRTC,labeller=label_value)+
  scale_x_continuous(breaks=c(15,30,45,60,125))+
  scale_color_manual(values=c("#093B6D", "#EF761B"))+
  scale_fill_manual(values=alpha(c("#093B6D", "#EF761B"),0.2))+
  labs(x="Weight (kg)",y="Area Under the Curve (mg*h/mL)")+
  theme_bw()+
  theme(strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background = ggplot2::element_rect(fill = "#475c6b"), 
        strip.text = ggplot2::element_text(face = "bold", color = "white"),
        legend.position = "bottom")+
  transition_states(
    REGIMEN,
    transition_length = 2,
    state_length = 1
  ) +
  enter_fade() + 
  exit_shrink() +
  ease_aes('sine-in-out')

Code

#anim_save("./boxplot.gif",type="cairo-png")