backup_options <- options() options(width = 3000) knitr::opts_chunk$set(collapse = TRUE, comment = "#>",echo = FALSE, warning = FALSE, message = FALSE, cache = FALSE, tidy = FALSE, size = "small") library(PFIM)
In this example, we simulate an 1-compartment model with linear elimination for IV infusion over 1 hour (inspired by [@Sukeishi2021]). One hundred and fifty (150) subjects receive a 400mg loading dose on the first day, followed by 4 daily doses of 200mg. Blood samples are taken at the end of the $1^{st}$ infusion (H1), H20, H44, H66 and H120. By evaluating this design, we will then select 4 sampling times on intervals (0,48) and (72,120) for an optimal design using PSO (Particle Swarm Optimization) algorithm. The same optimization problem will also be run with PGBO (Population Based Genetic Optimization) and Simplex algorithm.
Reports of the design evaluation and optimization are available at https://github.com/iame-researchCenter/PFIM
Linear1InfusionSingleDose_ClV
from the library of modelsmodelEquations = list( outcomes = list( "RespPK"), modelFromLibrary = list("PKModel" = "Linear1InfusionSingleDose_ClV") )
modelParameters = list( ModelParameter( name = "V", distribution = LogNormal( mu = 50, omega = sqrt( .26 ) ) ), ModelParameter( name = "Cl", distribution = LogNormal( mu = 5, omega = sqrt( .34 ) ) ) )
RespPK
errorModelRespPK = Combined1( outcome = "RespPK", sigmaInter = 0.5, sigmaSlope = sqrt( 0.15 ) ) modelError = list( errorModelRespPK )
administrationRespPK = Administration( outcome = "RespPK", Tinf = rep( 1, 5 ), timeDose = seq( 0, 96, 24 ), dose = c( 400, rep( 200, 4 ) ) )
RespPK
samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 1,12,24,44,72,120 ) )
arm1
of size 150 with administration administrationRespPK
and samplings samplingTimesRespPK
arm1 = Arm( name = "arm1", size = 150, administrations = list( administrationRespPK ) , samplingTimes = list( samplingTimesRespPK ) )
arm1
to the design design1
design1 = Design( name = "design1", arms = list( arm1 ) )
evaluationPop = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" ), designs = list( design1 ), fim = "population", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationFIMPop = run( evaluationPop ) evaluationInd = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" ), designs = list( design1 ), fim = "individual", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationFIMInd = run( evaluationInd ) evaluationBay = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" ), designs = list( design1 ), fim = "Bayesian", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationFIMBay = run( evaluationBay )
show( evaluationFIMPop ) show( evaluationFIMInd ) show( evaluationFIMBay ) fisherMatrix = getFisherMatrix( evaluationFIMPop ) getCorrelationMatrix( evaluationFIMPop ) getSE( evaluationFIMPop ) getRSE( evaluationFIMPop ) getShrinkage( evaluationFIMPop ) getDeterminant( evaluationFIMPop ) getDcriterion( evaluationFIMPop ) plotOptions = list( unitTime = c("hour"), unitOutcomes = c("mcg/mL") ) plotEvaluation = plotEvaluation( evaluationFIMPop, plotOptions ) plotSensitivityIndice = plotSensitivityIndice( evaluationFIMPop, plotOptions ) SE = plotSE( evaluationFIMPop ) RSE = plotRSE( evaluationFIMPop ) # examples plotOutcomesEvaluationRespPK = plotEvaluation[["design1"]][["arm1"]][["RespPK"]] plotSensitivityIndice_RespPK_Cl = plotSensitivityIndice[["design1"]][["arm1"]][["RespPK"]][["V"]] SE = SE[["design1"]] RSE = RSE[["design1"]]
outputPath = "C:/..." plotOptions = list( unitTime=c("hour"), unitResponses= c("mcg/mL") ) outputFile = "Example02_EvaluationPopFIM.html" Report( evaluationFIMPop, outputPath, outputFile, plotOptions ) outputFile = "Example02_EvaluationIndFIM.html" Report( evaluationFIMInd, outputPath, outputFile, plotOptions ) outputFile = "Example02_EvaluationBayFIM.html" Report( evaluationFIMBay, outputPath, outputFile, plotOptions )
We create sampling times that will be used in the initial design for comparison during the optimization process. As PSO does not optimize the dose regimens, we keep the same administration.
samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 1, 48, 72, 120 ) )
We define sampling times constraints that aim to select 4 sampling times: 2 within the interval [1,48], and 2 within [72, 120], with at least a delay of 5 between two points.
samplingConstraintsRespPK = SamplingTimeConstraints( outcome = "RespPK", initialSamplings = c( 1, 48, 72, 120 ), numberOfTimesByWindows = c( 2, 2 ), samplingsWindows = list( c( 1, 48 ), c( 72, 120 ) ), minSampling = 5 )
arm2 = Arm( name = "arm2", size = 150, administrations = list( administrationRespPK ) , samplingTimes = list( samplingTimesRespPK ), samplingTimesConstraints = list( samplingConstraintsRespPK ) ) design2 = Design( name = "design2", arms = list( arm2 ), numberOfArms = 150 )
optimization = Optimization( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "PSOAlgorithm", optimizerParameters = list( maxIteration = 100, populationSize = 10, personalLearningCoefficient = 2.05, globalLearningCoefficient = 2.05, seed = 42, showProcess = T ), designs = list( design2 ), fim = "population", outcomes = list( "RespPK") )
optimizationPSO = run( optimization )
show( optimizationPSO ) fisherMatrix = getFisherMatrix( optimizationPSO ) getCorrelationMatrix( optimizationPSO ) getSE( optimizationPSO ) getRSE( optimizationPSO ) getShrinkage( optimizationPSO ) getDeterminant( optimizationPSO ) getDcriterion( optimizationPSO )
outputPath = "C:/..." outputFile = "Example02_OptimizationPSOPopFIM.html" Report( optimizationPSO, outputPath, outputFile, plotOptions )
optimizationPGBO = Optimization( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "PGBOAlgorithm", optimizerParameters = list( N = 30, muteEffect = 0.65, maxIteration = 1000, purgeIteration = 200, seed = 42, showProcess = TRUE ), designs = list( design2 ), fim = "population", outcomes = list( "RespPK") )
optimizationPGBO = run( optimizationPGBO ) fisherMatrix = getFisherMatrix( optimizationPGBO ) getCorrelationMatrix( optimizationPGBO ) getSE( optimizationPGBO ) getRSE( optimizationPGBO ) getShrinkage( optimizationPGBO ) getDeterminant( optimizationPGBO ) getDcriterion( optimizationPGBO )
show( optimizationPGBO )
outputPath = "C:/..." outputFile = "Example02_OptimizationPGBOPopFIM.html" Report( optimizationPGBO, outputPath, outputFile, plotOptions )
optimizationSimplex= Optimization( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "SimplexAlgorithm", optimizerParameters = list( pctInitialSimplexBuilding = 20, maxIteration = 1000, tolerance = 1e-6, showProcess = TRUE ), designs = list( design2 ), fim = "population", outcomes = list( "RespPK") ) optimizationSimplex = run( optimizationSimplex )
show( optimizationSimplex ) fisherMatrix = getFisherMatrix( optimizationSimplex ) getCorrelationMatrix( optimizationSimplex ) getSE( optimizationSimplex ) getRSE( optimizationSimplex ) getShrinkage( optimizationSimplex ) getDeterminant( optimizationSimplex ) getDcriterion( optimizationSimplex )
outputPath = "C:/..." outputFile = "Example02_OptimizationSimplexPopFIM.html" Report( optimizationSimplex, outputPath, outputFile, plotOptions ) saveRDS( optimizationSimplex, file = paste0( outputPath, "Example02_OptimizationSimplexPopFIM.rds" ) )
options(backup_options)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.