R Speaks Non Linear Mixed Effects Modeling, RsNLME, is a suite of R packages and supplementary Shiny apps developed by Certara that supports pharmacometric modeling inside R.

When used together, these packages can add a level of speed and flexibility to your pharmacometrics workflow that cannot be achieved with point-and-click type tools.

Whether you’re still learning R or a seasoned expert, efficiently build, execute, and analyze your models from the Shiny GUI, then generate the corresponding R code to reproduce and expand your workflow from the R command line.


Certara.RsNLME uses tidyverse syntax to generate and update pharmacometric models in R. All functions can easily be chained together using the %>% operator from magrittr.


# Create two-compartment pharmacokinetic model
model <- pkmodel(numCompartments = 2, data = pkData, 
                ID = "Subject", Time = "Act_Time", A1 = "Amount", CObs = "Conc",
                modelName = "TwCpt_IVBolus_FOCE_ELS") 

# Update initial estimates              
model <- model %>%
          fixedEffect(effect = c("tvV", "tvCl", "tvV2", "tvCl2"), value = c(15, 5, 40, 15)) 

Model Builder

If you are still learning the command line syntax, use the Certara.RsNLME.ModelBuilder Shiny application from either RStudio or Pirana to build your NLME model from the GUI and generate the corresponding R and PML code to reproduce your model across multiple environments.


model <- modelBuilderUI(data = pkData)

Learn more about Certara.RsNLME.ModelBuilder here

Fit model

Next, we can execute the above model we created from the Shiny GUI inside R using the command fitmodel():

job <- fitmodel(model)

   Scenario RetCode    LogLik     -2LL      AIC      BIC nParm nObs nSub EpsShrinkage Condition
1: WorkFlow       1 -632.7953 1265.591 1283.591 1308.057     9  112   16      0.17297   3.34287

   #Scenario Parameter   Estimate Units    Stderr       CV%     2.5%CI    97.5%CI Var.Inf.factor
1:  WorkFlow       tvV 15.4073105    NA 1.1100322  7.204581 13.2058234 17.6087976             NA
2:  WorkFlow      tvCl  6.6015335    NA 0.8196260 12.415691  4.9759987  8.2270683             NA
3:  WorkFlow      tvV2 41.1964591    NA 1.0869697  2.638503 39.0407112 43.3522071             NA
4:  WorkFlow     tvCl2 14.0376052    NA 1.0477128  7.463615 11.9597141 16.1154964             NA
5:  WorkFlow      CEps  0.1611031    NA 0.0199441 12.379710  0.1215487  0.2006575             NA

   Label         nV       nCl          nV2       nCl2
1:    nV 0.06927885 0.0000000 0.000000e+00 0.00000000
2:   nCl 0.00000000 0.1809428 0.000000e+00 0.00000000
3:   nV2 0.00000000 0.0000000 4.357874e-08 0.00000000
4:  nCl2 0.00000000 0.0000000 0.000000e+00 0.04265183

   Label nV nCl nV2 nCl2
1:    nV  1   0   0    0
2:   nCl  0   1   0    0
3:   nV2  0   0   1    0
4:  nCl2  0   0   0    1

       Label        nV        nCl       nV2      nCl2
1: Shrinkage 0.1198519 0.01591489 0.9988085 0.2864404

   Label        nV        nCl        nV2       nCl2
1:    nV 0.0394158 0.00000000 0.0000e+00 0.00000000
2:   nCl 0.0000000 0.05395801 0.0000e+00 0.00000000
3:   nV2 0.0000000 0.00000000 1.7392e-05 0.00000000
4:  nCl2 0.0000000 0.00000000 0.0000e+00 0.05324145

Model Executor

Or alternatively, use the Certara.RsNLME.ModelExecutor Shiny application to specify additional engine arguments, change the estimation algorithm, add output tables, and more - all from the Shiny GUI!



Learn more about Certara.RsNLME.ModelExecutor here

Xpose NLME

After executing the model, we use the Certara.Xpose.NLME, xpose, ggplot2, and flextable packages to generate our model diagnostic plots.


xpdb <- xposeNlme(dir = "~/TwCpt_IVBolus_FOCE_ELS")

res_vs_idv(xpdb) + 

Learn more about Certara.Xpose.NLME here

Model Results

Or alternatively, use the Certara.ModelResults Shiny application to easily preview, customize, and report model diagnostics plots and tables from the Shiny GUI. Furthermore, the application will generate the corresponding .R and .Rmd code to reproduce your model diagnostics.

Learn more about Certara.ModelResults here

Fit VPC model

Lastly, users can execute a vpcmodel() to generate a Visual Predictive Check (VPC) plot and assess model fit.

Using the Certara.RsNLME package, we will execute the function vpcmodel() to return our observed and simulated data used to generate our VPC.


vpcJob <- vpcmodel(model)

Next we’ll extract our observed and simulated data from the return value of vpcmodel().

obs_data <- vpcJob$predcheck0

sim_data <- vpcJob$predout


Then we can use the tidyvpc package to parameterize our VPC.


 vpc <- observed(obs_data, y = DV, x = IVAR) %>%
      simulated(sim_data, y = DV) %>%
      binless() %>%

Learn more about tidyvpc here

VPC Results

Or alternatively, use the Certara.VPCResults Shiny application to easily parameterize, customize, and report VPC plots from the Shiny GUI. Furthermore, the application will generate the corresponding .R and .Rmd code to reproduce your VPC’s in R.

vpcResultsUI(obs_data, sim_data)

Learn more about Certara.VPCResults here