Error in Shiny app containing differential equation solving tool shiny


#1

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
	    
	  })




#2

Thanks for posting and for including your example! Can you include the error message / logs as well?


#3

Logs/ error messages are as follows-

Listening on http://127.0.0.1:5412
Warning: Error in reactive:all.data: object 'Q' not found
Stack trace (innermost first):
112: reactive:all.data [C:\Users\ms1225\Documents\NONMEM articles\RShiny\App codes\desolve\query/server.R#152]
101: all.data
100: ggplot
99: renderPlot [C:\Users\ms1225\Documents\NONMEM articles\RShiny\App codes\desolve\query/server.R#254]
89: reactive:plotObj
78: plotObj
77: origRenderFunc
76: output$plotCONC
1: runApp


#4

That error seems to point very clearly to server.R#152 (that's line 152 of your server.R file) and server.R#254 as well. It also looks like you are somewhere referencing an object Q that does not exist in whatever context it is being referenced (object 'Q' not found). It looks to me like there are a handful of places that could be, so it's probably worth just paring down your example a bit so that you can find / isolate the place where that mistaken reference is happening.

This article might help explore some debugging best practices.


#5

Thank you for your help. I used the debugging tools to rectify most of my code.

Here is my updated code

#ui.r
fixedPage(

	#Application Title 
	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"),   #BSA Based dosing. This will show the output of the dose to be 
			                      #administered on usual regimen 500 mg/ m2
			
			#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"),          

			#Radio button to choose sex of patient
			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 Albumin
			sliderInput("ALB",
			            "Albumin (g/dL):",
			            min = 2.4,
			            max = 4.8,
			            value = 2.5,
			            step = 0.1),
			
			
			#Slider input for number of individuals to carry out simulation on
			sliderInput("n",
						"Number of Individuals:",
						min = 10,
						max = 1000,
						value = 10,
						step = 10),
			
			br(),
			
			#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)    #for application
	library(ggplot2)  #plotting 
	library(deSolve)  #differential equation solver
	library(plyr)     #Split and rearrange data, ddply function
	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
	
# Creates a set of copies of R running in parallel and communicating over sockets.
	#package- parallel
	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- These functions provide several ways to parallelize computations using a cluster.
	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 concentrations 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({

	    #----------------------------------------------------------------------------------------	
	    #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"
	    
	    #Input for weight from UI
	    WT <- input$Weight
	    
	    #Input for sex from UI
	    SEX <- input$ISM
	    
	    #Input for CRCL form UI
	    CRCL<- input$CRCL
	    
	    #Input of Albumin level from UI
	    ALB<- input$ALB
	    
	    #Input for BSA from UI
	    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
	    
	    #Epsilons (as SD)
	    EPS1SD <- 0.287751281  #Proportional residual error
	    EPS2SD <- 0.643975931  #Additive residual error	
	    
	    OMEGA <- matrix(c(0.271,0.102,0,0,0.102,0.257,0.291,0,0,0.291,1.5,0,0,0,0,0.774),4,4)
	    
	    #Now use multivariate rnorm to turn the covariance matrix into ETA values
	    #mvrnorm- package MASS- Produces one or more samples from the specified 
	    #multivariate normal distribution
	    #n-number of samples required
	    #mu- vector giving means of the variables
	    #Sigma-omega-a positive-definite symmetric matrix 
	    #specifying the covariance matrix of the variables.
	    
	    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
	    #SEX=male
	    #Clearance
	     if(SEX==1){
	       par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1*COV2
	     }
	    #SEX=Female
	    #Clearance
	    if(SEX==0){
	      par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1
	    }
	    #Central volume
	    par.data$V1 <- V1POP*exp(ETA2)*(WT/58)*(ALB/3.6)^COV3
	    
	    #Peripheral volume
	    par.data$V2 <- V2POP*exp(ETA3)*(WT/58)
	    
	    #Inter-compartment clearance
	    par.data$Q  <- QPOP*exp(ETA4)*(WT/58)^0.75
	    
	    #Calculate rate-constants for differential equation solver
	    par.data$K12 <- par.data$Q/par.data$V1
	    par.data$K21 <- par.data$Q/par.data$V2
	    par.data$K10 <- par.data$CL/par.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
	   
	    #------------------------------------------------------------------------------------------	
	    #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        #h
	    CRATE <- CDOSE/CTinf   #mg/h
	    
	    
	    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)
	    #approxfun- package stats, Return a list of points which linearly interpolate given data points, or a function 
	    #performing the linear (or constant) interpolation.
	    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
	    #specified above
	    
	    #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	
	    #package- compiler-provide an interface to a byte code compiler for R
	    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
	      #lsoda- package desolve
	      #Solving initial value problems for stiff or non-stiff systems of 
	      #first-order ordinary differential equations (ODEs).
	      #The R function lsoda provides an interface to the FORTRAN ODE solver
	      
	      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,V2), simulate.conc.cmpf, .parallel=TRUE)
	    
	    #Calculate individual concentrations in the central compartment	
	    sim.data$IPRE <- sim.data$A1/sim.data$V1
	    
	    #---------------------------------------------------------------------------------------------------------------------------------------------
	    
	    #simulating observations (Y)
	    sim.data$DV <- sim.data$IPRE*(1 + EPS1) + EPS2
})
	    #---------------------------------------------------------------------------------------------------------------------------------------------
	  #   #Step 5.  Draw a plot of the simulated data
	   CIlow <- function(x) quantile(x, probs = 0.05)
	   CIhi <- function(x) quantile(x, probs = 0.95)

   output$plotCONC <- renderPlot({
   plotobj <- ggplot(data=all.data())
   plotobj <- plotobj + geom_line(aes(x = TIME, y = DV), colour = "red", size = 1)
   #plotobj <- plotobj + geom_ribbon(aes(x = time, ymin = CIlow, ymax = CIhi), fill = "red", alpha = 0.3)
   plotobj <- plotobj + scale_y_continuous("Concentration (microg/mL) \n", lim=c(0,600))
   plotobj <- plotobj + scale_x_continuous("\nTime (hours)", breaks=c(0,8,16,24))
   print(plotobj)

	    })
	    
	    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
	    
	  }) #bracket closing shiny server

Error logs

Listening on http://127.0.0.1:6899
Warning: <anonymous>: ... may be used in an incorrect context: ‘.fun(piece, ...)’

Warning: <anonymous>: ... may be used in an incorrect context: ‘.fun(piece, ...)’

Warning: Error in : ggplot2 doesn't know how to deal with data of class numeric
Stack trace (innermost first):
    105: fortify.default
    104: fortify
    103: structure
    102: ggplot.data.frame
    101: ggplot.default
    100: ggplot
     99: renderPlot [C:\Users\ms1225\Downloads\Shiny_code_2/server.R#222]
     89: <reactive:plotObj>
     78: plotObj
     77: origRenderFunc
     76: output$plotCONC
      1: runApp

Basically, I'm running a model that has been described by differential equations to R for simulating different scenarios.

I don't understand the current error.

Thank you !! :slight_smile:


#6

This is a common R error - it basically means you are passing the wrong type into the function (i.e. maybe ggplot2 is expecting a data.frame, but you are passing it a numeric vector - the ol' square peg vs. round hole). I definitely recommend testing functions outside of the reactive context before pushing things into Shiny. Debugging in a Shiny app is much more tricky than in an interactive session, so it is worth testing things interactively first.

Besides that, google searching the error messages can be helpful!


#7

Thank you so much! Will do! :smile: