Shiny applications are a great way to show your complicated models in an interactive way. In recent years, many different examples have been published online showcasing a wide range of functionality. Shiny applications are also being increasingly used in pharmacometrics, to quickly see the effect of a different dosing regimen or even develop a Bayesian TDM application.
Getting started
If you are completely new to Shiny, I can recommend visiting a few pages to get the basics:
As a basis, a shiny application consists of two parts/scripts. The user interface (UI) and the server script.
UI script: Consists of all the graphical components needed to interact with the model. This can be buttons, checkboxes, textboxes, labels, graphical output, tables, etc.
Server script: Main script where all the calculations are being performed, being updated by the input from the user. This scripts consists of input components (retrieved from the UI script), reactive expressions (see later), and outputs (such as graphs and tables). We will use the script from Part 1 as the basis for our Shiny application.
Reactive expressions
In order to change our previous script for PK simulations with mrgsolve to work with Shiny, we need to make use of reactive expressions. Reactive expressions are expressions, ´objects´, that can change based on the user input. For example, the dose, the parameter estimates, the simulation time are all examples of objects that need be included in a reactive expression.
R scripts – ui.R
If you use RStudio, you can easily create a new Shiny Web Application which would automatically generate the ui.R and the server.R script. When you click the Run App button, this default Shiny app is being generated:
Even though the Old Faithful Geyser Data is interesting, we have different problems to solve… However, this app provides a basis which we can further adapt. Like in part 1 of this blog, we need to define which options we need to have to fulfill all of our criteria. You can check with the R script of part 1 to see how I decided what I need.
Select from a list:
- Oral or intravenous administration
Button:
- Download the simulated dataset
Numerical input:
- Population parameters
- ETA’s and Sigma’s
- Dose
- Number of samples for simulation
- Simulation time
- Min and max probabilities for figure
Graph
- The PK simulation
We can now start updating the sidebarlayout in our ui.R script to incorporate all the components. Previously, it contained a sliderInput to modify the number of bins used in the histogram.
# Sidebar with a slider input for number of bins sidebarLayout( sidebarPanel( sliderInput("bins", "Number of bins:", min = 1, max = 50, value = 30) ),
Now, we want to change this to include headers, buttons, labels, checkboxes, and numerical input. Let’s start with the dosing details. We need to select the route of administration and be able to set the dose:
# Sidebar sidebarLayout( sidebarPanel( h4("Dosing details"), selectInput('admin', label='Route of administration', c("Depot","I.V. bolus")), numericInput("dose", label='Dose', value = 10,min=0) ),
Then, we add the simulation details and a number of numeric inputs for all the parameters. Disclaimer: I fully agree that the design of the interface can be much better but it does the job.
# Sidebar sidebarLayout( sidebarPanel( h4("Dosing details"), hr(), selectInput('admin', label='Route of administration', c("Depot","I.V. bolus")), numericInput("dose", label='Dose', value = 100,min=0), br(), h4("Simulation details"), hr(), numericInput("nsim", label='Number of samples for simulation', value = 1000,min=0), numericInput("simtime", label='Simulation time', value = 12,min=0), numericInput("minprob", label='Minimal probability', value = 0.25,min=0,max=1), numericInput("maxprob", label='Maximal probability', value = 0.75,min=0,max=1), br(), h4("Population parameters"), hr(), numericInput("ka", label='Absorption rate constant', value = 10,min=0), numericInput("v1", label='Central volume', value = 10,min=0), numericInput("v2", label='Peripheral volume', value = 10,min=0), numericInput("q1", label='Inter-compartmental clearance', value = 10,min=0), numericInput("cl", label='Clearance', value = 10,min=0), br(), h4("IIV"), hr(), numericInput("etaka", label='ETA - Absorption rate constant', value = 10,min=0), numericInput("etav1", label='ETA - Central volume', value = 10,min=0), numericInput("etav2", label='ETA - Peripheral volume', value = 10,min=0), numericInput("etaq1", label='ETA - Inter-compartmental clearance', value = 10,min=0), numericInput("etacl", label='ETA - Clearance', value = 10,min=0), br(), h4("Residual variability"), hr(), numericInput("sigmaprop", label='Proportional', value = 0.01,min=0), numericInput("sigmaadd", label='Additive', value = 0.1,min=0) ),
Now, in the main panel, we want to have the graph as an output and add the button here where we can download our simulated dataset. Luckily Shiny knows that this is a common feature which is why we can use the downloadButton function.
mainPanel( plotOutput("simPlot"), downloadButton("downloadData", "Download") )
This code means that we need to have two output objects in our server.R script, the simPlot, and the downloadData.
Remember that you need to provide a unique name to all the components that you are adding to your application. We can then call these components from the server script to see the user input.
R scripts – server.R
The RStudio template contains just a simple graphical output to start with. No reactive expressions. Nothing complicated.
library(shiny) # Define server logic required to draw a histogram shinyServer(function(input, output) { output$distPlot <- renderPlot({ # generate bins based on input$bins from ui.R x <- faithful[, 2] bins <- seq(min(x), max(x), length.out = input$bins + 1) # draw the histogram with the specified number of bins hist(x, breaks = bins, col = 'darkgray', border = 'white') }) })
Now we want to start implementing our previous script into the server side. The way I approach this is by identifying which parts of the code need to be updated when a user changes the input. Unfortunately for us, almost the complete script is dependent on the user input. From the dose that is given to the model that is used.
As a guideline, every reactive expression results in a single dataframe or object that will be updated. In our script, we however need to have two outputs: The dataframe with the simulated values, which we need to be able to download, and the dataframe that will be used for the generation of the graph.
For starters, we create a reactive expression in our script in which we copy the R script until we have a dataframe with the simulations:
sim_dataframe <- reactive({ ############################################### ### Dosing oral <- T dose <- 100 # mg ############################################### ## Parameter estimates ka <- 1 # Absorption rate constant cl <- 1 # Clearance vd <- 1 # Central distribution volume vd2 <- 1 # Peripheral distribution volume q1 <- 1 # Inter-compartmental clearance (set to 0 for 1-cmt model) ############################################### ## Insert etas and sigmas etaka <- 0.09 etacl <- 0.09 etavd <- 0.09 etavd2 <- 0.09 etaq1 <- 0.09 ################################################# ## Insert residual error sigmaprop <- 0.01 # Proportional error sigmaadd <- 0 # Additive error ###############################3 ## Simulation info nsamples <- 1000 ### Number of simulated individuals sim_time <- 12 ## Time of simulation # set probabilities of ribbon in figure minprobs=0.1 maxprobs=0.9 ################################################### ############# Set Dosing objects if(oral){ ## Oral dose Administration <- as.data.frame(ev(ID=1:nsamples,ii=24, cmt=1, addl=9999, amt=dose, rate = 0,time=0)) }else{ ## IV BOLUS Administration <- as.data.frame(ev(ID=1:nsamples,ii=24, cmt=2, addl=9999, amt=dose, rate = 0,time=0)) } ## Sort by ID and time data <- Administration[order(Administration$ID,Administration$time),] ###################################################################### ### Load in model for mrgsolve mod <- mread_cache("popPK") ## Specify the omegas and sigmas matrices omega <- cmat(etaka, 0, etacl, 0,0,etavd, 0,0,0,etavd2, 0,0,0,0,etaq1) sigma <- cmat(sigmaprop, 0,sigmaadd) ## Set parameters in dataset data$TVKA <- ka data$TVCL <- cl data$TVVC <- vd data$TVVP1 <- vd2 data$TVQ1 <- q1 ################################################################# ###### Perform simulation out <- mod %>% data_set(data) %>% omat(omega) %>% smat(sigma) %>% mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE) # Simulate 100 observations, independent of the total simulated time ### Output of simulation to dataframe df <- as.data.frame(out) })
If there is anything unclear about this script, go back to part 1.
We can then change the objects that are currently fixed based on the user input. We can do this by retrieving the information with input$'component name'. For example, input$ka will give the value of the numericinput of the absorption rate constant.
The same piece of code with this modifications can be found below:
sim_dataframe <- reactive({ ############################################### ### Dosing if(input$admin =='Depot'){ oral <- T }else{ oral <- F } dose <- input$dose ############################################### ## Parameter estimates ka <- input$ka # Absorption rate constant cl <- input$cl # Clearance vd <- input$v1 # Central distribution volume vd2 <- input$v2 # Peripheral distribution volume q1 <- input$q1 # Inter-compartmental clearance (set to 0 for 1-cmt model) ############################################### ## Insert etas and sigmas etaka <- input$etaka etacl <- input$etacl etavd <- input$etav1 etavd2 <- input$etav2 etaq1 <- input$etaq1 ################################################# ## Insert residual error sigmaprop <- input$sigmaprop # Proportional error sigmaadd <- input$sigmaadd # Additive error ###############################3 ## Simulation info nsamples <- input$nsim ### Number of simulated individuals sim_time <- input$simtime ## Time of simulation # set probabilities of ribbon in figure minprobs=input$minprob maxprobs=input$maxprob ################################################### ############# Set Dosing objects if(oral){ ## Oral dose Administration <- as.data.frame(ev(ID=1:nsamples,ii=24, cmt=1, addl=9999, amt=dose, rate = 0,time=0)) }else{ ## IV BOLUS Administration <- as.data.frame(ev(ID=1:nsamples,ii=24, cmt=2, addl=9999, amt=dose, rate = 0,time=0)) } ## Sort by ID and time data <- Administration[order(Administration$ID,Administration$time),] ###################################################################### ### Load in model for mrgsolve mod <- mread_cache("popPK") ## Specify the omegas and sigmas matrices omega <- cmat(etaka, 0, etacl, 0,0,etavd, 0,0,0,etavd2, 0,0,0,0,etaq1) sigma <- cmat(sigmaprop, 0,sigmaadd) ## Set parameters in dataset data$TVKA <- ka data$TVCL <- cl data$TVVC <- vd data$TVVP1<- vd2 data$TVQ1 <- q1 ################################################################# ###### Perform simulation out <- mod %>% data_set(data) %>% omat(omega) %>% smat(sigma) %>% mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE) # Simulate 100 observations, independent of the total simulated time ### Output of simulation to dataframe df <- as.data.frame(out) return(df) })
This piece of code will result in the reactive dataframe called sim_dataframe which will be updated every time a user changes the input. If we want to use this dataframe later in our script, we need to use the following syntax: sim_dataframe()
Since objects cannot be transferred between reactive expressions, we need to define the min and maxprobs again here. We use the simulated dataset for the summary calculations as follows:
sum_stat <- reactive({ ################################################################# ## Calculate summary statistics # set probabilities of ribbon in figure minprobs=input$minprob maxprobs=input$maxprob sum_stat <- sim_dataframe() %>% group_by(time) %>% summarise(Median_C=median(DV), Low_percentile=quantile(DV,probs=minprobs), High_percentile=quantile(DV,probs=maxprobs) ) return(sum_stat) })
Notice how I use the sim_dataframe() as the first chunk of the piping function.
Finally, our output, the graph that we want is similar to the code as we developed in part 1, with the only difference being the brackets behind sum_stat(). In this piece of code, we just say that the output$simPlot in the UI, needs to contain the generated ggplot object.
output$simPlot <- renderPlot({ ggplot(sum_stat(), aes(x=time,y=Median_C)) + ## Add ribbon for variability geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+ ## Add median line geom_line(size=2) + # scale_y_log10()+ # Set axis and theme ylab(paste("Concentration",sep=""))+ xlab("Time after dose (h)")+ theme_bw()+ # Remove legend theme(legend.position="none") })
So what about our second output, the download button? As I mentioned previously, Shiny is capable of handling these kind of things with limited lines of code. We can specify a default filename and we specify the contents (in our case the sim_dataframe()).
# Downloadable csv of selected dataset ---- output$downloadData <- downloadHandler( filename = function() { "Simulated_dataset.csv" }, content = function(file) { write.csv(sim_dataframe(), file, row.names = FALSE) } )
Congratulations! You just created your pharmacometric Shiny application! As you may have noticed, starting with the development of a basis script, without any Shiny components, saves a lot of time in the end. Development of a Shiny application can result in some strange error messages and you do not want to waste a lot of time debugging your scripts.
Fancy user interface
In order to get away from the very standard Shiny template, a very interesting package is shinythemes. This allows you to change the look of your application with just the click of a button.
Deploying your Shiny application for free
Information on how to make your app publicly available can be found on the following website:
https://shiny.rstudio.com/deploy/
The free package on shinyapps.io allows 25 hours of use per month.
Conclusion
Lets take a look of what the functionality of the Shiny app should be, which we specified in Part 1:
- Pharmacokinetic simulations of a simple 2-compartment model with oral or I.V. administration
- Manually specify the parameter estimates (theta’s, omega’s and sigma’s)
- Generate a figure of the simulation
- Have the option to save the simulated data to a .csv file
As you now can see from your own application, we fulfilled all the points! Of course, there are many remaining items that can be improved. For example, now, when a single number is being changed, the whole simulation is being repeated. You can add an action button to your application which holds all actions until you press it (especially for large simulations this makes sense). Also, I load the model code in the reactive expression, which means that it is being loaded on every change (whereas the model remains the same...). Furthermore, you can write a code to generate pdf reports from your app, by using RMarkdown. The options are limitless.
If you know the basics, it is just a small step from the app that we have now created to the one that is located on my website. Good luck!
You can download the created ui.R and server.R scripts below
COMMENT
Any suggestions or typo’s? Leave a comment or contact me at info@pmxsolutions.com!
Hi,
Would you be able to upload the R scripts for “simulation application for pharmacokinetic models”?
Thank you,
Nova
Hi Nova,
Could you clarify which scripts you would like to see? The R scripts of this post are located at the bottom and of course in the text.
The simulation application on the top menu is located at my GitHub:
https://www.pmxsolutions.com/2019/12/08/pmx-sim-is-now-open-source/
Michiel