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.


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

# Get the project
demoPath = '/lixoft/simulx/simulx2020R1/demos/'
project <- paste0(demoPath, "2.models/longitudinal.smlx")
loadProject(projectFile = project)

# Run the simulation

# The results are accessible through the function getSimulationResults()
# The results for the output is the list res and TS is one of the outputs
  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

where “/” is the user’s home folder (on windows C:/Users/toto if toto is your username).


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.

initializeLixoftConnectors(software = "simulx")

importMonolixProject(projectFile = "C:/Users/username/lixoft/Monolix/monolix2020R1/demos/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran")

# check that mlx_EBEs element has properly been created when importing the project

# 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 = 1, interval = 0.1, final = 100), output = "AUC"))

# check the simulation group that already exists
# => 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 output for AUC
setGroupElement(group = "simulationGroup1", elements = c("mlx_EBEs", "AUCregular"))

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

# run simulation

# get the results: simulated values for AUC on grid 1:0.1:100
# for the 12 individuals of the original data set using the EBEs and original treatment
sim <- getSimulationResults()$res$AUC

# plot the results
ggplot(data=sim, aes(x=time, y=AUC))+
geom_line() + facet_wrap(~id)


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.

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

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

# Define vector of population parameters
                        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, 

# 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)
  element = paste0("covTable",i,".csv"))

# Define common treatment
defineTreatmentElement(name="1000nmol", element=list(admID=1, data=data.frame(time=seq(0,by=21,length.out = 5), 
                                                                              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.
covTable <- read.csv("covTable1.csv")
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)
defineTreatmentElement(name="1000nmolHomo_800nmolHetero", element=list(admID=1, data="trtTableGenotype.csv"))

# 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") 

# 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 the full example to correct a VPC for dropout is given 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.

initializeLixoftConnectors(software = "simulx")

importMonolixProject(projectFile = "PD_TTE.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"))
saveProject(projectFile = "PD_TTE_alltimes.smlx")



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(dplyr) # for data reformatting

simulxDemosPath <- "C:/Users/username/lixoft/simulx/simulx2020R1/demos/"
project <- paste0(simulxDemosPath, "6.outcome_endpoints/6.1.outcome_endpoints/OutcomeEndpoint_PKPD_changeFromBaseline.smlx")
Nrep <- 200
sample_sizes <- c(15, 30, 50, 75, 100)

loadProject(projectFile = project)

# Focus on comparing the two groups 400mg_q4w_ and 800mg_q8w_ for the output "Cholesterol_measured_atCtrough" 
# -> remove other groups and outputs
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)

   # run the simulation

   # 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

 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")+
 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)

   success <- T
   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")

modal close image