Select Page

# Examples using R functions

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

Minimal examples based on demo projects are also included in the documentation of all Simulx functions.

## Load a project and run the simulation

Simple example to load an existing project, run the simulation and look at the results :

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

# Get the project
project <- paste0(getDemoPath(), "/2.models/longitudinal.smlx")

# Run the simulation
runSimulation()

# The results are accessible through the function getSimulationResults()
# The results for the output is the list res and TS is one of the outputs
head(getSimulationResults()$res$TS)


  id time       TS
1  1    0 10.00000
2  1    1 11.04939
3  1    2 12.20862
4  1    3 13.48915
5  1    4 14.90359
6  1    5 16.46585


## Import a Monolix project and simulate a new output variable with the EBEs

In this example, the Monolix demo theophylline_project.mlxtran is imported to Simulx. A new variable AUC is added to the model, and it is simulated on a regular grid using the EBEs estimated in the Monolix project. This requires that the tasks “Population parameters” and “EBEs” have run in the Monolix project.

# Import a Monolix project and simulate a new output variable with the EBEs ####

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

#Load project from Monolix Demos and run it to get EBEs estimated==================================================

MonolixProject <- paste0(getDemoPath(), "/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran") loadProject(MonolixProject) runScenario() initializeLixoftConnectors(software = "simulx", force = TRUE) #IMPORT Monolix project and check that mlx_EBEs element has been created ###=================================================== importProject(MonolixProject) getIndividualElements() #SET ADDITIONAL formula in the model to compute the AUC=================================================== setAddLines("ddt_AUC = Cc") # define a new element output with a regular grid for the variable "AUC" defineOutputElement(name="AUCregular", element = list( data = data.frame( start = 0, interval = 1, final = 120), output = "AUC")) #CHECK the simulation group that already exists=================================================== getGroups() # => we currently have one group, called "simulationGroup1", which uses the following elements:
# - mlx_Pop (population parameter values) => need to be replaced by mlx_EBEs
# - mlx_Adm1 (treatment from monolix data set) => OK
# - mlx_CONC (output CONC with times as in monolix data set) => need to be replaced by AUC with new grid
# - size = 12 (same number of in original data set) => OK

#MODIFY the existing simulation group such that it uses EBEs and the new output element ===================================================

setGroupElement(group = "simulationGroup1",
elements = c("mlx_EBEs", "AUCregular"))

# check the new simulation setup
getGroups()

#Set shared IDs to say the individuals in the treatment element and in the parameter element are the same and should alwyas be sampled together ===================================================

setSharedIds(c("treatment","individual"))

# check the simulation setup
getGroups()
getSharedIds()

# - mlx_EBEs => ok
# - remaining parameters (error model parameters if residual error is simulated) => not needed as we output model predictions
# - AUCregular => ok
# - size=12 => ok

#SAVE and RUN project===================================================

saveProject("simulx_api_example.smlx")

# run simulation
runSimulation()

#GET the results: simulated values for AUC and sampled individual parameters===================================================
simulatedParam <- getSimulationResults()$IndividualParameters$simulationGroup1

simulatedOutput <- getSimulationResults()$res$AUC

#PLOT the results===================================================
library(ggplot2)
ggplot(data = simulatedOutput, aes(x=time, y=AUC))+
geom_line() + aes(color = factor(as.double(original_id)))+
labs(colour="Original ID")



We get the following simulations for the new output variable. “Original_id” is the ID used in the original dataset loaded in the Monolix project.

## Create a Simulx project from scratch and simulate covariate-dependent treatments

This script creates a project similar to the demo project treatment_weight_and_genotype_based.smlx, available in the folder 3.1.treatments of Simulx demos. The model file TMDDmodel.txt can be downloaded here.

While in the interface of Simulx it is possible to directly define distribution laws for covariates and covariate-dependent treatments that are computed at the simulation step, in R it is necessary to sample covariates from distributions and derive the corresponding doses before defining the elements for the simulation.

###############################################################
# Create a Simulx project from scratch and simulate covariate-dependent treatments
###############################################################

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

# Create a new project based on the model
newProject(modelFile = "TMDDmodel.txt")

# Define vector of population parameters
definePopulationElement(name="PopParameters",
element = data.frame(V_pop=3, omega_V=0.1 ,Cl_pop=0.1, omega_Cl=0.1, Q_pop=1, omega_Q=0.2,
V2_pop=3, omega_V2=0.2, beta_V_logtWeight=1, KD_pop=0.01, omega_KD=0.01,
R0_pop=0.2, omega_R0=0.3, kint_pop=50, omega_kint=0.2, kon_pop=10,
omega_kon=0.01, ksyn_pop=10, omega_ksyn=0.2, beta_Cl_logtWeight=0.75,
beta_Q_logtWeight=0.75, beta_V2_logtWeight=0.75,
beta_KD_Genotype_Heterozygous=1.3))

# Define covariates: in the interface it is possible to define distribution laws,
# but in R we have to sample from the distributions prior to the simulation.
# Here we define 3 tables for the 3 simulation groups
for (i in 1:3){
covTable <- data.frame(id=1:250,
Weight = rlnorm(n=250, mean=log(70), sd=0.3),
Genotype=c("Homozygous", "Heterozygous")[sample(1:2, 250, replace=T)])
write.csv(covTable, file = paste0("covTable",i,".csv"), quote=F, row.names=F)
defineCovariateElement(name=paste0("covTable",i),
element = paste0("covTable",i,".csv"))
}

# Define common treatment
amount=1000, tInf=0.208)))

# Define weight-based treatment: in the interface it can be done directly with a scaling formula,
# but in R we have to define a table of individual doses prior to the simulation.
# We use covTable2.csv for simulation group 2.
trt_14nmolPerKg <- data.frame(id=rep(1:250, each=5),
time=rep(seq(0,by=21,length.out = 5), times=250),
amount=14*covTable$Weight, tInf = 0.208) write.csv(trt_14nmolPerKg, file="trtTableWeight.csv", quote=F, row.names=F) defineTreatmentElement(name="14nmolPerKg", element=list(admID=1, data="trtTableWeight.csv")) # Define genotype-based treatment: in the interface it can be done directly with a scaling formula, # but in R we have to define a table of individual doses prior to the simulation. # We use covTable3.csv for simulation group 3. covTable <- read.csv("covTable3.csv") trt_genotype <- data.frame(id=rep(1:250, each=5), time=rep(seq(0,by=21,length.out = 5), times=250), amount=ifelse(covTable$Genotype=="Heterozygous", 1000, 800), tInf = 0.208)
write.csv(trt_genotype, file="trtTableGenotype.csv", quote=F, row.names=F)

# Define outputs
defineOutputElement(name="Concentration", element = list(output="L", data=data.frame(time=seq(0,96,by=1))))
defineOutputElement(name="TargetOccupancy", element = list(output="TO", data=data.frame(time=seq(0,96,by=1))))

# By default there is a single simulation group named simulationGroup1,
# change its name and set all elements for the simulation
setGroupSize("simulationGroup1", 250)
setGroupElement("simulationGroup1", elements = c("PopParameters","14nmolPerKg","Concentration","TargetOccupancy", "covTable1"))
renameGroup("simulationGroup1", "Weight_based")

# Define the two other simulation groups with different covariates and treatments.
# By default the other elements are the same as the first group
setGroupElement("Flat_dose", elements = c("1000nmol","covTable2"))
setGroupElement("Genotype_based", elements = c("1000nmolHomo_800nmolHetero","covTable3"))

# Save project and run the simulation
saveProject(projectFile = "treatment_weight_and_genotype_based.smlx")
runSimulation()

# Get simulation results: a table of individual parameters for each group, and a table of simulated values for each output
sim <- getSimulationResults()
names(sim$IndividualParameters) # "Weight_based", "Flat_dose", "Genotype_based" names(sim$res) # "L", "TO"



## Import a Monolix project and simulate with new output times

Monolix simulates the original dataset with replicates to generate the VPC. It can be useful to perform the same simulations in Simulx with modified observation times, for example to correct the bias in VPC due to dropout. The example below shows this step, and uses the project PDTTE_dropout. The full example to correct a VPC for dropout is detailed on this page.

In this example the Monolix project with the biased VPC is named PD_TTE.mlxtran, the new regular measurement times are similar to the measurements in the original dataset, but are not impacted by dropout, and the number of replicates is 500 (as the default number of simulations for the VPC in Monolix). The random variable capturing dropouts in the model used in Monolix is called survival.

###############################################################
# Import a Monolix project and simulate with new output times
###############################################################

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

importProject(projectFile = "PDTTE_dropout.mlxtran")
defineOutputElement(name="y1_alltimes", element = list(output="y1", data=data.frame(time=seq(-100,1500,by=100))))
setGroupElement("simulationGroup1", elements = c("y1_alltimes","mlx_survival"))
setNbReplicates(500)
saveProject(projectFile = "PD_TTE_alltimes.smlx")
setPreferences(exportsimulationfiles=T)
runSimulation()



## Run a simulation with different sample sizes and plot the power of the study based on an endpoint

This piece of code uses the demo project OutcomeEndpoint_PKPD_changeFromBaseline.smlx (demo folder 6.1), with the goal of comparing two treatment arms (400mg Q4W and 800mg Q8W) with a range of different sample sizes (number of simulated individuals per arm). The power of the clinical trial is defined as the chance of having a significant difference in the mean change from baseline between the two arms. This chance is computed over 200 replicates for each sample size. The script then plots the power of the clinical trial with respect to the sample size, with a target power of 85% displayed as a horizontal line.

Note that Simulx2020R1 has R functions to modify the project and run the simulations, but no function to define the outcomes and endpoints like in the interface. Outcomes and endpoints must be defined directly in R as postprocessing of the simulation results, as shown on this example.

library(ggplot2)
library(dplyr) # for data reformatting
library(lixoftConnectors)
initializeLixoftConnectors(software="simulx")

project <- paste0(simulxDemosPath, "6.outcome_endpoints/6.1.outcome_endpoints/OutcomeEndpoint_PKPD_changeFromBaseline.smlx")
Nrep <- 200
sample_sizes <- c(15, 30, 50, 75, 100)

# Focus on comparing the two groups 400mg_q4w_ and 800mg_q8w_ for the output "Cholesterol_measured_atCtrough"
# -> remove other groups and outputs
removeGroup("400mg_q8w_")
removeGroupElement("400mg_q4w_", "Cholesterol_prediction")
removeGroupElement("800mg_q8w_", "Cholesterol_prediction")
removeGroupElement("400mg_q4w_", "mAbTotal_prediction")
removeGroupElement("800mg_q8w_", "mAbTotal_prediction")

outcomes <- NULL

for(N in sample_sizes){

# change the number of individuals for each arm
setGroupSize("400mg_q4w_", N)
setGroupSize("800mg_q8w_", N)
setNbReplicates(Nrep)

# run the simulation
runSimulation()

# get the results: output is measured cholesterol (variable yLDLc)
sim <- getSimulationResults()$res$yLDLc

# calculate the outcome: remove baseline and take the mean
outcome <- sim %>% group_by(rep,id,group) %>% # for each id
summarise(diff=yLDLc-first(yLDLc)) %>% # relative to baseline with "difference"
slice(2:n()) %>% # remove first record which is zero
summarise(chFromBL=mean(diff)) %>% # mean/max/min per id
ungroup()

outcomes <- rbind(outcomes, cbind(outcome, N=N))

}
endpoint <- outcomes %>% group_by(rep, group, N) %>% summarise(meanChFromBL=mean(chFromBL)) %>% ungroup()

print(ggplot(endpoint, aes(group,meanChFromBL)) +
geom_boxplot(aes(fill = group)) +
facet_grid(. ~ N) +
ylab("Mean change from baseline at Ctroughs")+
xlab("")+
theme(axis.text.x = element_blank(), axis.ticks = element_blank()))


successDF <- NULL # loop over number of individuals per arm
for (n in unique(outcomes$N)){ # loop over the replicates for (i in unique(outcomes$rep)) {

outcome400mg <- outcomes[outcomes$group=="400mg_q4w_" & outcomes$rep==i & outcomes$N==n,'chFromBL'] outcome800mg <- outcomes[outcomes$group=="800mg_q8w_" & outcomes$rep==i & outcomes$N==n,'chFromBL']

tt <- t.test(outcome400mg,outcome800mg, alternative = "two.sided", mu = 0)

if(tt$p.value<0.05){ success <- T }else{ success <- F } successDF <- rbind(successDF, data.frame(SampleSize=n, Replicate=i, Success= success)) } } summary <- as.data.frame(successDF %>% group_by(SampleSize) %>% summarise(Power=sum(Success)/n()*100) %>% ungroup()) ggplot(summary)+geom_line(aes(x=SampleSize,y=Power, color="red"))+ ylab("Power (%)")+xlab("Number of subjects per arm")+geom_hline(yintercept = 85)+ ylim(c(0,100))+ theme(legend.position = "none") ## Handling of warning/error/info messages Error, warning and info messages from Simulx are displayed in the R console when performing actions on a Simulx project. They can be hidden via the R options. Set lixoft_notificationOptions$errors, lixoft_notificationOptions$warnings and lixoft_notificationOptions$info to 1 or 0 to respectively hide or show the messages.

#### Example

op = options()
op$lixoft_notificationOptions$warnings = 1   #hide the warning messages
options(op)

### Force software switch

By default, function initializeLixoftConnectors() prompts users to confirm that they want to proceed with the software switch, in order to avoid losing unsaved changes in the currently loaded project. To override this behavior, force argument can be set to TRUE when calling the function. However, the default behavior can be changed globally as well.

#### Example

op = options()
op\$lixoft_lixoftConnectors_forceSoftwareSwitch <- TRUE
options(op)