Hi!
I'm running a differential equation containing model for a drug and something seems to be amiss. The UI contains all the parameters that could be altered and used in the model simulation which will show up as a plot.
The error is very non-specific. Any help would be appreciated.
attaching code for your reference.
#ui.r
fixedPage(
#Application Title and Logo
fixedRow(
column(10,
h2("Drug Simulation tool"),
h6("Drug simulation tool"),
offset = 1, align = "center")
), #Brackets closing "fixedRow"
hr(), #Add a break with a horizontal line
####################################################################################
#Sidebar panel with widgets
sidebarLayout(
sidebarPanel(
#Slider input for BSA
sliderInput("BSA",
"BSA (m2):",
min = 1.08,
max = 2.16,
value = 1.08,
step = 0.01),
textOutput("DOSE"),
#Slider input for Dose
sliderInput("CDOSE",
"Dose administered as short infusion (mg):",
min = 500,
max = 1000,
value = 1,
step = 1),
#Numeric output for infusion rate
textOutput("RATE"),
radioButtons("ISM",
"Sex:",
choices = list("Female" = 0,
"Male" = 1),
selected = 1),
#Slider input for weight
sliderInput("Weight",
"Weight (kg):",
min = 0,
max = 150,
value = 12,
step = 1),
#Slider input for CRCL
sliderInput("CRCL",
"CRCL (ml/min):",
min = 0,
max = 150,
value = 12,
step = 5),
#Slider input for ALB
sliderInput("ALB",
"Albumin (g/dL):",
min = 2.4,
max = 4.8,
value = 2.5,
step = 0.1),
#Slider input for number of individuals
sliderInput("n",
"Number of Individuals:",
min = 10,
max = 1000,
value = 10,
step = 10),
br(),
#Selection box for prediction intervals
selectInput("CI",
"Prediction Intervals %:",
choices = list("80% - 10th and 90th percentiles" = 1,
"90% - 5th and 95th percentiles" = 2,
"95% - 2.5th and 97.5th percentiles" = 3),
selected = 1),
#Button to initiate simulation
submitButton("Simulate"),
align = "left"), #Brackets closing "sidebarPanel"
mainPanel(
#Plot output for concentration-time profile
plotOutput("plotCONC", height = 650, width = 750),
align = "center") #Brackets closing "mainPanel"
) #Brackets closing "sidebarLayout"
) #Brackets closing "fixedPage"
#server.r
#Load package libraries
library(shiny)
library(ggplot2)
library(deSolve)
library(plyr)
library(reshape2)
library(compiler)
library(doParallel)
library(grid) #Plotting
library(MASS) #mvrnorm function
library(MBESS) #cor2cov function
#Code for functions and variables which are not reactive (not dependent on "input$X")
#Setting up cores to run parallel processes, thus increasing speed
#Set up a cluster of cores to run the application over
cl <- makePSOCKcluster(detectCores() - 1)
#detectCores() searches for the number of cores that the local machine has
#Contents with makePSOCKcluster brackets can be changed to a whole number if you
#want to assign an exact number
#List packages that are required to be sent to each core for the parallel process
#The foreach package always needs to be included
#This example uses the .parallel argument in ddply which calls a function that uses
#lsoda from the deSolve package
clusterEvalQ(cl, list(library(deSolve), library(foreach)))
#Registers the parallel backend with the foreach package (automatically loaded when
#doParallel is loaded)
registerDoParallel(cl)
#ggplot2 theme for plotting
theme_custom <- theme_set(theme_grey(18))
#TIME range - times where a concentration will be calculated
time <- seq(from = 0, to = 24, by = 0.25)
#----------------------------------------------------------------------------------------
#Define user-input dependent functions for output
shinyServer(function(input, output){
#Reactive expression to generate a reactive data frame
#This is called whenever the input changes
all.data <- reactive({
#Calculate x% Confidence Intervals
CI <- input$CI
#Assign probabilities for 80% confidence intervals
if (CI == 1) {
CIlow <- 0.1
CIhi <- 0.9
}
#Assign probabilities for 90% confidence intervals
if (CI == 2) {
CIlow <- 0.05
CIhi <- 0.95
}
#Assign probabilities for 95% confidence intervals
if (CI == 3) {
CIlow <- 0.025
CIhi <- 0.975
}
#Function for calculating median, upper and lower confidence intervals for x
sumfuncx <- function(x) {
stat1 <- median(x)
stat2 <- quantile(x, probs=CIlow, names=F)
stat3 <- quantile(x, probs=CIhi, names=F)
stat4 <- length(x)
result <- c("median"=stat1, "low"=stat2, "hi"=stat3, "n"=stat4)
result
}
#----------------------------------------------------------------------------------------
#Step 1. Create a parameter dataframe with ID and parameter values for each individual
#Define a population
n <- input$n
par.data <- seq(from = 1, to = n, by = 1)
par.data <- data.frame(par.data)
names(par.data)[1] <- "ID"
WT <- input$Weight
SEX <- as.factor(input$ISM)
CRCL<- input$CRCL
ALB<- input$ALB
BSA<- input$BSA
#Define parameter values
#Thetas
CLPOP <- 3.14687 #Clearance L/hour
V1POP <- 5.46589 #Central volume L
V2POP <- 6.75549 #Peripheral volume L
QPOP <- 5.7451 #Inter-compartmental clearance L/h
COV1 <- 0.520245 #Covariate effect of CRCL on clearance (power model)
COV2 <- 1.22484 #Covariate effect of gender on CL
COV3 <- -1.51075 #Covariate effect of ALB on V1
#Omegas (as SD)
ETA1SD <- 0.520860826
ETA2SD <- 0.507793265
ETA3SD <- 1.222178383
ETA4SD <- 0.881369956
#Specify a correlation matrix for ETA's
R12 <- 0.387845089 #Correlation coefficient for CL-V1
R13 <- 0 #Correlation coefficient for CL-V2
R23 <- 0.462686988 #Correlation coefficient for V1-V2
R14 <- 0 #Correlation coefficient for CL-Q
R24 <- 0 #Correlation coefficient for V1-Q
R34 <- 0 #Correlation coefficient for V2-Q
#Epsilons (as SD)
EPS1SD <- 0.287751281 #Proportional residual error
EPS2SD <- 0.643975931 #Additive residual error
#Calculate ETA values for each subject
CORR <- matrix(c(1,R12,R13, R14,R12,1,R23,R24,R13,R23,1,R34,R14,R24,R34,1),4,4)
#Specify the between subject variability for CL, V1, V2
SDVAL <- c(ETA1SD,ETA2SD,ETA3SD,ETA4SD)
#Use this function to turn CORR and SDVAL into a covariance matrix
OMEGA <- cor2cov(CORR,SDVAL)
#Now use multivariate rnorm to turn the covariance matrix into ETA values
ETAmat <- mvrnorm(n = n, mu = c(0,0,0,0), OMEGA)
ETA1 <- ETAmat[,1]
ETA2 <- ETAmat[,2]
ETA3 <- ETAmat[,3]
ETA4 <- ETAmat[,4]
#Define individual parameter values
if(SEX==1){
par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1*COV2
}
if(SEX==0){
par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1
}
par.data$V1 <- V1POP*exp(ETA2)*(WT/58)*(ALB/3.6)^COV3
par.data$V2 <- V2POP*exp(ETA3)*(WT/58)
par.data$Q <- QPOP*exp(ETA4)*(WT/58)^0.75
#Calculate rate-constants for differential equation solver
K12 <- Q/V1
K21 <- Q/V2
K10 <- CL/V1
#------------------------------------------------------------------------------------------
#Step 2. Make a time sequence and specify the dose information for a system of differential equations
#There are a number of ways that doses can be coded for the deSolve package - see help for deSolve
#Creating continuous infusion (for 24 hours) - this use the approxfun function to make a "forcing function" for infusion rate in the differential equations
CDOSE <- input$CDOSE #mg
CTinf <- 0.16
CRATE <- CDOSE/CTinf #mg/h
#hours
CTIMEinf <- c(0,CTinf,24) #100 - make sure the function works long after the 24 hours
CRATEinf <- c(CRATE,0,0)
#Define an interpolation function that returns rate when given time - "const"
#i.e. at TIME = 0, RATE = CRATE. At TIME = 0.16, RATE = 0 (infusion stopped)
Cstep.doseinf <- approxfun(CTIMEinf, CRATEinf, method = "const")
#The time sequence must include all "event" times for deSolve so add them here and sort
#Do not repeat a time so use "unique" as well
#Function containing differential equations for amount in each compartment
DES <- function(T, A, THETA) {
RateC <- Cstep.doseinf(T) #Infusion rate
K12 <- THETA[1]
K21 <- THETA[2]
K10 <- THETA[3]
dA <- vector(length = 2)
dA[1] = RateC -K12*A[1] +K21*A[2] -K10*A[1] #Central compartment
dA[2] = K12*A[1] -K21*A[2] #Peripheral compartment
list(dA)
}
#Compile DES function - it's called by lsoda for each individual in the dataset
DES.cmpf <- cmpfun(DES)
#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 3. Run the differential equation solver for each patient in par.data
#Function for simulating concentrations for the ith patient
simulate.conc <- function(par.data) {
#List of parameters from input for the differential equation solver
THETAlist <- c("K12" = par.data$K12,"K21 "= par.data$K21,"K10" = par.data$K10)
#Set initial compartment conditions
A_0 <- c(A1 = 0, A2 = 0)
#Run differential equation solver for simulated variability data
sim.data <- lsoda(A_0, time, DES.cmpf, THETAlist)
sim.data <- as.data.frame(sim.data)
}
#Compile simulate.conc function - it's called by ddply for each individual in the dataset
simulate.conc.cmpf <- cmpfun(simulate.conc)
#Apply simulate.conc.cmpf to each individual in par.data
#Whilst maintaining their individual values for V1, SEX and WT for later calculations
sim.data <- ddply(par.data, .(ID,V1), simulate.conc.cmpf, .parallel=TRUE)
#Calculate individual concentrations in the central compartment
sim.data$IPRE <- sim.data$A1/sim.data$V1
#---------------------------------------------------------------------------------------------------------------------------------------------
# Add residual error
#Use random number generator to simulate residuals from a normal distribution
nobs <- n*length(time) #number of observations = number of subjects * number of time points per subject
EPS1 <- rnorm(nobs, mean = 0, sd = EPS1SD) #Proportional residual error
EPS2 <- rnorm(nobs, mean = 0, sd = EPS2SD) #Additive residual error
sim.data$DV <- sim.data$IPRE*(1 + EPS1) + EPS2
statsCONC <- ddply(sim.data, .(time), function(sim.data) sumfuncx(sim.data$DV))
names(statsCONCS)[c(2,3,4)] <- c("median","low","hi")
sim.data<-statsCONCS
})
#---------------------------------------------------------------------------------------------------------------------------------------------
output$plotCONC <- renderPlot({
#
# theme_bw2 <- theme_set(theme_bw(base_size = 16))
# theme_bw2 <- theme_update(plot.margin = unit(c(1.1,1.1,3,1.1), "lines"),
# axis.title.x=element_text(size = 16, vjust = 0),
# axis.title.y=element_text(size = 16, vjust = 0, angle = 90),
# strip.text.x=element_text(size = 14),
# strip.text.y=element_text(size = 14, angle = 90))
#
plotobj <- ggplot(all.data())
plotobj <- plotobj + stat_summary(aes(x = time, y = DV), fun.ymin = CIlow, fun.ymax = CIhi, geom = "ribbon", fill = "blue", alpha = 0.2)
plotobj <- plotobj + stat_summary(aes(x = time, y = IPRE), fun.ymin = low, fun.ymax = hi, geom = "ribbon", fill = "red", alpha = 0.3)
plotobj <- plotobj + stat_summary(aes(x = time, y = IPRE), fun.y = median, geom = "line", size = 1, colour = "red")
plotobj <- plotobj + theme_bw2
plotobj <- plotobj + scale_y_continuous("Concentration (mg/L) \n")
plotobj <- plotobj + scale_x_continuous("\nTime (hours)")
print(plotobj)
#Use custom ggplot2 theme
})
output$DOSE <- renderText({
paste("Dose to be administered =", signif(input$BSA*500, digits = 3) ,"mg")
}) #Brackets closing "renderText" function
output$RATE <- renderText({
paste("Infusion rate =", signif(input$CDOSE*6, digits = 3) ,"mg/hr")
}) #Brackets closing "renderText" function
})