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
- 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.