RTutorial - Calibayes

CaliBayes R tutorial

Preliminaries

Before starting this tutorial, you need

A working installation of the R statistical programming language - version 2.9 or later.
The R libraries calibayesR, basisR, and SBMLModels - see SupportPackages for installation instructions;

   1 library(calibayesR)
   2 library(SBMLModels)
   3 data(LogisticModel)

The code below will run automatically via the commands:

   1 library(calibayesR)
   2 demo(YeastGrowthDemo)

The Model

This demo is based on calibayesR version 0.1.8

To demonstrate the CaliBayes system, we will consider a simple logistic model.

See ?LogisticModel for further details.

Using R to generate your initial XML files

First we create our MCMC settings

   1 #Which simulators are available
   2 wsdl = "http://calibayes1.ncl.ac.uk:81/CaliBayesService_v2.wsdl"
   3 listSimulatorMethods(wsdl)
   4 
   5 #Create settings file
   6 settings = createCaliBayesSettings(wsdl, burn=10, thin=2, simulator="copasi-deterministic")

We can save this file for use with our web-services using the saveCaliBayesSettings function.

Next we construct our experiment data, yeast_data:

   1 #Just use the first experiment and a copy of points for illustration
   2 yeast_data = YeastGrowth[[1]][1:3,]
   3 
   4 #Plot of the experiment data for the simple decay model
   5 plot(yeast_data$time, yeast_data$value, xlab="Time", ylab="Population level")
   6 
   7 #yeast_data$value are observations with error, 
   8 #where the error has a Normal distribution with mean 0, and precision tau

Now we will create the prior distribution:

   1 #The number of particles
   2 n = 1000
   3 #Create prior for the initial species level
   4 species=data.frame(S=rlnorm(n,3,sqrt(0.02)))
   5 
   6 #Create prior distribution for the measurement error structure
   7 distributions = c('Gaussian')
   8 errors = data.frame(S.tau=rlnorm(n,- 11, 1))
   9 
  10 #Create priors for the parameters
  11 parameters = data.frame(K=rlnorm(n,9,0.2), r=rlnorm(n,1,0.6))
  12 
  13 #Create a distribution object that contains our prior
  14 prior = createCaliBayesDistribution(parameters, species, distributions, errors, FALSE)
  15 plot(prior)

To calibrate the model to the data we do:

   1 #Create an Experiments object
   2 experiments = createCaliBayesExperiment(yeast_data, species)
   3 
   4 #Use the calibrate web-service
   5 sid = calibrate(wsdl, LogisticModel, settings, experiments, prior, asText=TRUE)
   6 
   7 #Check CaliBayes to see if calibration has finished
   8 isCaliBayesReady(wsdl, sid)
   9 
  10 #Check CaliBayes every 60 seconds
  11 isCaliBayesFinished(wsdl, sid, 10, 10)
  12 
  13 #Once CaliBayes has finished, retrieve the posterior
  14 posterior = getPosterior(wsdl, sid)
  15 
  16 #Check the output
  17 plot(posterior)
  18 
  19 #Compare the posterior to the prior
  20 compareDistributions(prior, posterior)

Loading your own SBML files into R

We have found that the best way to load an .xml file into R is as follows. This example loads an SBML file called "LogisticModel.xml" from the current working directory:

   1 SBML = toString.XMLNode(xmlInternalTreeParse("LogisticModel.xml"))