Bayesian individual dynamics predictions
This partial example corresponds to Bayesian individual dynamic predictions with uncertainty. This is what it does:
estimates population parameters for a model,
fixes the population parameters to the estimates,
changes the dataset for another, that can be filtered up to some landmark time,
runs the task Conditional distribution that samples sets of parameters from the individual posterior distributions,
simulates the model based on the sampled parameters with Simulx.
library("lixoftConnectors")
initializeLixoftConnectors(software = "monolix")
project <- paste0(getDemoPath(), "/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran")
loadProject(projectFile = project)
# run parameter estimation (if needed)
runPopulationParameterEstimation()
# set pop params to estimates and fix pop params
setInitialEstimatesToLastEstimates()
popparams <- getPopulationParameterInformation()
popparams$method <- "FIXED"
setPopulationParameterInformation(popparams)
# the original dataset is truncated at TIME=8 and used as new dataset in the monolix project
baseDataset <- read.table(getData()$dataFile, header=T)
filteredDataset <- baseDataset[baseDataset$TIME<8,]
write.csv(filteredDataset, "filteredDataset.csv", row.names = F, quote=F)
setData("filteredDataset.csv", getData()$headerTypes, getData()$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)
# save the project to avoid overwriting previous results
saveProject("bayesian_forecasting.mlxtran")
# run SAEM (needed even if pop param are fixed) 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
# retrieve the samples from the conditional distribution
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)
# export project from Monolix to Simulx
exportProject(settings=list(targetSoftware="simulx"), force=T)
# define an element from an external file containing the samples from the conditional distribution
defineIndividualElement(name="simulatedParameters", element="table_simulated_parameters.csv")
# use the defined element for the simulation
setGroupElement(group = "simulationGroup1", elements = c("simulatedParameters"))
# set the number of "individuals" (actually nb indiv * nb rep)
setGroupSize(group = "simulationGroup1", size = nrow(simparams))
# run the simulation
runSimulation()
# retrieve the results
simresults <- getSimulationResults()