Examples using R functions

This page provides pieces of code to perform some procedures with the R functions for Monolix. The code has to be adapted to use your projects and use the results as you wish in R.

Load and run a project
Convergence assessment
Profile likelihood
Bayesian individual dynamic predictions
Covariate search
Generate plots

Load and run a project

Simple example to load an existing project, run SAEM and print the number of iterations used in the exploratory and smoothing phases:

# load and initialize the API
library(lixoftConnectors) 
initializeLixoftConnectors(software="monolix")

demoPath = 'C:/Users/username/lixoft/monolix/monolix2020R1/demos/1.creating_and_using_models/1.1.libraries_of_models/'
project <- paste0(demoPath, "theophylline_project.mlxtran"
loadProject(projectFile = project)

runPopulationParameterEstimation()
iter <- getSAEMiterations()
print(paste0("Iterations in exploratory phase: ",iter$iterationNumbers[1]))
print(paste0("Iterations in smoothing phase: ",iter$iterationNumbers[2]))

Convergence assessment

Example of the functions to call to build a project from scratch using one of the demo data sets and a model from the libraries, and run a convergence assessment to evaluate the robustness of the convergence. Running the convergence assessment in R gives the flexibility to choose a different strategy than the one done in the interface, in this case it samples new initial values for all parameters instead of fixed effects only, and it does not change the random seed:

# load and initialize the API
library(lixoftConnectors)
initializeLixoftConnectors(software="monolix")

# create a new project by setting a data set and a structural model
# replace  by the path to your home directory
demoPath = paste0(getDemoPath(), '/1.creating_and_using_models/1.1.libraries_of_models/')
librariesPath = 'C:/ProgramData/Lixoft/MonolixSuite2021R2/factory/library/pk'
newProject(data = list(dataFile = paste0(demoPath,'data/warfarin_data.txt'),
                       headerTypes =c("id", "time", "amount", "observation", "obsid", "contcov", "catcov", "ignore"),
                       mapping = list("1" = "y1")),
           modelFile = paste0(librariesPath,'/oral1_1cpt_TlagkaVCl.txt'))

# set tasks in scenario
scenario <- getScenario()
scenario$tasks = c(populationParameterEstimation = T, 
                   conditionalModeEstimation = T, 
                   conditionalDistributionSampling = T, 
                   standardErrorEstimation=T, 
                   logLikelihoodEstimation=T)
scenario$linearization = TRUE
setScenario(scenario)

# ----------------------------------------------------------------------------
# convergence assessment: run 5 estimations with different initial estimates,
# store the results in tabestimates
# ----------------------------------------------------------------------------
popparams <- getPopulationParameterInformation()
tabestimates <- NULL; tabiters <- NULL
for(i in 1:5){
  # sample new initial estimates
  popini <- sapply(1:nrow(popparams), function(j){runif(n=1, min=popparams$initialValue[j]/2, max=popparams$initialValue[j]*2)})
  
  # set sampled values as new initial estimates
  newpopparams <- popparams
  newpopparams$initialValue <- popini
  setPopulationParameterInformation(newpopparams)
  
  # run the estimation
  runScenario()
  
  # store the estimates and s.e. in a table
  estimates <- as.data.frame(getEstimatedPopulationParameters())
  names(estimates) <- "estimate"
  rses <- getEstimatedStandardErrors()$linearization$rse
  names(rses) <- getEstimatedStandardErrors()$linearization$parameter
  rses <- as.data.frame(rses)
  estimates <- merge(estimates, rses, by = "row.names")
  estimates$run <- i
  names(estimates)[names(estimates) == "Row.names"] <- "param"
  tabestimates <- rbind(tabestimates, estimates)
  
  # store the iterations
  iters <- getChartsData("plotSaem")
  iters$run <- i
  tabiters <- rbind(tabiters, iters)
}

# load plotting libraries
library(ggplot2)
library(gridExtra)

# plot SAEM iterations
plotList <- list()
i <- 1
for (param in popparams$name) {
  if (popparams[popparams$name == param, ]$method == "FIXED") next
  changePhase <- tabiters$iteration[which(diff(tabiters$phase) == 1) + 1]
  plotList[[i]] <- ggplot(tabiters, aes_string(x = "iteration", y = param)) +
    geom_line(aes(group = run, color = factor(run))) +
    theme(legend.position = "none", plot.title = element_text(hjust = .5)) +
    geom_vline(xintercept = changePhase, color = 1:length(changePhase)) +
    labs(title = param, x = NULL, y = NULL)
  i <- i + 1
}
grid.arrange(grobs = plotList, ncol = 3)

# plot population parameters
plotList <- list()
i <- 1
for (param in popparams$name) {
  if (popparams[popparams$name == param, ]$method == "FIXED") next
  estimates <- tabestimates[tabestimates$param == param, ]
  plotList[[i]] <- ggplot(estimates, aes(x = run, y = estimate)) +
    geom_point(aes(color = factor(run))) +
    geom_errorbar(aes(ymax = estimate * (1+rses/100), ymin = estimate * (1-rses/100), color = factor(run))) +
    theme(legend.position = "none", plot.title = element_text(hjust = .5)) +
    labs(title = param, x = NULL, y = NULL)
  i <- i + 1
}
grid.arrange(grobs = plotList, ncol = 3)

Profile likelihood

This code performs a profile likelihood for the parameter f_LDLc_pop. This code which can be modified allows more flexibility than using the function confintmlx predefined in the Rsmlx package.

library(lixoftConnectors)
initializeLixoftConnectors()

loadProject("path_to_project.mlxtran")

tested_values=rev(c(seq(0.5,0.975,by=0.025),0.99))

for(i in tested_values){
  setPopulationParameterInformation(f_LDLc_pop = list(initialValue = i, method = "FIXED"))
  runScenario()
  LL=getEstimatedLogLikelihood()$importanceSampling[1]
  df <- data.frame(param='f_LDLc_pop’, value=i, LL=LL)
}

Bayesian individual dynamic predictions

This partial example corresponds to bayesian individual dynamic predictions with uncertainty. This is what it does:

estimate population parameters for a model,
fix the population parameters to the estimates,
change the dataset for another, that can be filtered up to some landmark time,
run the task Conditional distribution that samples sets of parameters from the individual posterior distributions,
simulate the model based on the sampled parameters with Simulx.

library("lixoftConnectors")
initializeLixoftConnectors(software = "monolix")

loadProject(projectFile = "path_to_projet_file.mlxtran")

# set pop params to estimates and fix pop params
runPopulationParameterEstimation()
setInitialEstimatesToLastEstimates()
popparams <- getPopulationParameterInformation()
popparams$method <- "FIXED"
setPopulationParameterInformation(popparams)

# here the data can be modified in the project to take into account data until landmark time
filtereddatfile <- "path to data filtered up to landmark time.csv"
BaseData <- getData()
setData(filtereddatfile, BaseData$headerTypes, BaseData$observationTypes)

# set MCMC settings: min number iter=1000, relative interval width=1, number of param per indiv=200
# this is to sample 200 parameters from 1000 iterations. Samples are uniformly spread on the iterations.
CondDistSettings <- getConditionalDistributionSamplingSettings()
setConditionalDistributionSamplingSettings(nbminiterations=500, ratio=1, nbsimulatedparameters=200)

# run SAEM and conditional distribution
runPopulationParameterEstimation() # mandatory before other tasks, but nb of iterations is null as all parameters are fixed
runConditionalDistributionSampling() # this is sampling from the posterior distribution for each individual

simparams <- getSimulatedIndividualParameters()
simparams$id <- row.names(simparams) # replace id by rank of {rep, id}
simparams$rep <- NULL
write.csv(simparams, file="table_simulated_parameters.csv", row.names = F)

# simulate based on simulated parameters using Simulx

initializeLixoftConnectors(software = "simulx", force=TRUE)

importMonolixProject(projectFile = "path_to_projet_file.mlxtran")

defineIndividualElement(name="simulatedParameters", element="table_simulated_parameters.csv")
setGroupElement(group = "simulationGroup1", elements = c("simulatedParameters"))
setGroupSize(group = "simulationGroup1", size = nrow(simparams))

runSimulation()
simresults <- getSimulationResults()

Covariate search

The automatic covariate search procedures can be launched using the lixoftConnectors. Note that when reloading the project (in the GUI or via the lixoftConnectors) with the 2020R1, the covariate search results are displayed but the settings are lost. From the 2021R1 version on, the settings will be saved and reloaded.

In the example below, we do the following steps:

load a base monolix project
add covariate transformations
get the default settings of the covariate search and modify them
specify which covariates and parameters to test along with relationships which are lock-in or out.
run the covariate search
get the results

As an example, we are using the demo 5.2 warfarin_covariate1_project.mlxtran as base project.

library(lixoftConnectors)
initializeLixoftConnectors(software="monolix") # makes the link to the MonolixSuite installation

# load base project (from the demos)
loadProject("warfarin_covariate1_project.mlxtran")

# wt and sex are defined as covariates. Add logtWT=log(wt/70)
addContinuousTransformedCovariate(logtWt = "log(wt/70)" )

# save the project with the additional covariate under a new name
saveProject("warfarin_covariate1_project_covsearch.mlxtran")

# get the default covariate search settings
set <- getModelBuildingSettings()

# modify the settings
set$strategy <- 'cossac' # "cossac" or "scm"
set$useSambaBeforeCossac <- FALSE # if true, corresponds to "covSAMBA-COSSAC" in the GUI
set$criterion <- "LRT" # "BIC" or "LRT"
set$threshold$lrt[1] <- 0.01 # forward p-value threshold for LRT
set$threshold$lrt[2] <- 0.001 # backward p-value threshold for LRT

# set logtWt and sex as covariates to test on all parameters except Tlag 
set$covariates <- c("sex","logtWt")
set$parameters <- c("ka","V","Cl")

# in addition lock-out (never test it) "sex on ka" and lock-in (have it in all tested models) "logtWt on Cl"
set$relationships[1,] <- c("ka", "sex", FALSE) # sex never included on ka
set$relationships[2,] <- c("Cl", "logtWt", TRUE) # logtWt always included on Cl
# all other relationships corresponding to set$covariates and set$parameters will be tested

# run the covariate search
runModelBuilding(settings=set)

# get the results. The best model is indicated with $bestModel = T
getModelBuildingResults()

Generate plots

Starting with the version 2021R1, the plots of Monolix and PKanalix are now available as ggplot objects in R, thanks to new functions in the package lixoftConnectors.