Package {CGNM}

CGNM: Cluster Gauss-Newton Method

Description

Cluster Gauss-Newton method (CGNM) finds multiple approximate minimizers of a nonlinear least squares problem argmin_x||f(x)-y*||, without assuming the minimizer is unique. In the context of model fitting, f is the model, x is the parameter vector, and y* is the observed data.

Typical workflow

1. Cluster_Gauss_Newton_method — fit the model; returns a CGNM_result list (fields X, Y, residual_history, initialX, runSetting).

2. Cluster_Gauss_Newton_Bootstrap_method (optional) — residual resampling bootstrap using the fit above, for uncertainty quantification.

3. acceptedApproximateMinimizers / acceptedIndices / bestApproximateMinimizers — select the subset of the found minimizers considered to have converged to (approximately) the same minimum sum of squared residuals.

4. table_parameterSummary — tabulate parameter estimates across the accepted minimizers.

5. plot_Rank_SSR, plot_goodnessOfFit, plot_paraDistribution_byHistogram, plot_paraDistribution_byViolinPlots, plot_parameterValue_scatterPlots, plot_profileLikelihood, plot_2DprofileLikelihood — visualize the fit (ggplot2 must be loaded by the user).

Conventions worth knowing before calling these functions

• "Unset" is NA, not NULL. Arguments such as lowerBound, upperBound, ParameterNames, and initialIterateMatrix default to NA to mean "not supplied". Passing NULL instead of NA may error rather than being treated as unset.

• saveLog=TRUE is the default for Cluster_Gauss_Newton_method() and writes iteration history to a CGNM_log (or CGNM_log_<runName>) folder under the current working directory. plot_profileLikelihood, compare_profileLikelihood, table_profileLikelihoodConfidenceInterval, plot_SSRsurface, suggestInitialLowerRange, and suggestInitialUpperRange all need that per-iteration history (a single final result normally only keeps the first and last iteration); they accept the folder path (a string, or vector of strings), the classic list (looked up from the object's own runName, so the log must still be on disk), or an outputS4 = TRUE object that retains this history in memory (see below) and needs no disk access at all.

• nonlinearFunction may be written either as vector-in/vector-out (one parameter set at a time) or matrix-in/matrix-out (one parameter set per row, for parallelization); see the Parallel computation section of vignette("CGNM-vignette", package = "CGNM").

• The classic return value is a plain list, unless you opt into outputS4 = TRUE. Cluster_Gauss_Newton_method, Cluster_Gauss_Newton_Bootstrap_method, and Cluster_Gauss_Newton_EBE_method all accept this argument (default FALSE) and return a CGNM_result-class S4 object instead, with every field still reachable via $/[[ exactly as on the list, so all postprocessing/plotting functions accept either form unmodified. This lets calling code confirm inherits(x, "CGNM_result") instead of assuming a bare list has the right shape. It also retains the per-iteration history needed for the profile-likelihood functions above (main fit, and the bootstrap run's own history too if Cluster_Gauss_Newton_Bootstrap_method is subsequently run on the object), so those functions work directly on the returned object with no saveLog/disk dependency at all.

See vignette("CGNM-vignette", package = "CGNM") for a full worked example, and the top-level CLAUDE.md / README.md in the package source repository for a condensed quick-start reference.

Author(s)

Maintainer: Yasunori Aoki yaoki@uwaterloo.ca

Authors:

• Yasunori Aoki yaoki@uwaterloo.ca

References

Aoki et al. (2020) Cluster Gauss-Newton method. Optimization and Engineering, 1-31. doi:10.1007/s11081-020-09571-2

Aoki and Sugiyama (2024) Cluster Gauss-Newton method for a quick approximation of profile likelihood: With application to physiologically-based pharmacokinetic models. CPT Pharmacometrics Syst Pharmacol. 13(1):54-67. doi:10.1002/psp4.13055

CGNM_result-class

Description

An optional, opt-in S4 alternative to the classic plain-list output of Cluster_Gauss_Newton_method, Cluster_Gauss_Newton_Bootstrap_method, and Cluster_Gauss_Newton_EBE_method (set outputS4 = TRUE on any of those functions to receive one instead of a list).

Every field that the classic list output ever contains (X, Y, residual_history, initialX, initialY, lambda_history, finalParameterCombinations, runSetting, the bootstrap fields, and the EBE fields) is a slot here, and $, $<-, [[, [[<-, and names() all work exactly as they do on the classic list. This means every existing postprocessing or plotting function in this package (acceptedApproximateMinimizers(), table_parameterSummary(), plot_goodnessOfFit(), etc.) accepts a CGNM_result S4 object as a drop-in replacement for the list, with no change in behavior. The only difference from the classic list is that inherits(x, "CGNM_result")/is(x, "CGNM_result") lets calling code (including automated tooling) confirm the object's identity, which a plain list cannot do.

Existing code that relies on the classic list output is unaffected: outputS4 defaults to FALSE everywhere.

When outputS4 = TRUE, two further slots are populated so that profile-likelihood functions (plot_profileLikelihood, compare_profileLikelihood, table_profileLikelihoodConfidenceInterval, plot_SSRsurface, suggestInitialLowerRange, suggestInitialUpperRange) can be used directly on the object with no dependency on saveLog/disk/working directory: X_history/Y_history hold every intermediate iteration of the main fit (normally discarded once the final result is assembled), and bootstrapIterationHistory holds the equivalent history from Cluster_Gauss_Newton_Bootstrap_method, if one was run. Both are NULL unless outputS4 = TRUE was used, so the memory cost of retaining them is opt-in.

Usage

## S4 method for signature 'CGNM_result'
x$name

## S4 replacement method for signature 'CGNM_result'
x$name <- value

## S4 method for signature 'CGNM_result'
x[[i, j, ...]]

## S4 replacement method for signature 'CGNM_result'
x[[i, j, ...]] <- value

## S4 method for signature 'CGNM_result'
names(x)

Arguments

Cluster_Gauss_Newton_Bootstrap_method

Description

Conduct residual resampling bootstrap analyses using CGNM.

Usage

Cluster_Gauss_Newton_Bootstrap_method(
  CGNM_result,
  nonlinearFunction,
  num_bootstrapSample = 200,
  indicesToUseAsInitialIterates = NA,
  bootstrapType = 1,
  outputS4 = FALSE,
  ...
)

Arguments

Value

list of a matrix X, Y,residual_history, initialX, bootstrapX, bootstrapY as well as a list runSetting.

1. X, Y, residual_history, initialX: identical to what was given as CGNM_result.

2. X: a num_bootstrapSample by n matrix which stores the the X values that was sampled using residual resampling bootstrap analyses (In terms of model fitting this is the parameter combinations with variabilities that represent parameter estimation uncertainties.).

3. Y: a num_bootstrapSample by m matrix which stores the nonlinearFunction evaluated at the corresponding bootstrap analyses results in matrix bootstrapX above. In the context of model fitting each row corresponds to the model simulations.

4. runSetting: identical to what is given as CGNM_result but in addition including num_bootstrapSample and indicesToUseAsInitialIterates.

Examples

##lip-flop kinetics (an example known to have two distinct solutions)

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation, num_iteration = 10, num_minimizersToFind = 100,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
lowerBound=rep(0,3), ParameterNames=c("Ka","V1","CL_2"), saveLog = FALSE)

CGNM_bootstrap=Cluster_Gauss_Newton_Bootstrap_method(CGNM_result,
     nonlinearFunction=model_analytic_function, num_bootstrapSample=100)

plot_paraDistribution_byHistogram(CGNM_bootstrap)

Cluster_Gauss_Newton_EBE_method

Description

obtain ebe empirical bayes estimate (EBE) using CGNM.

Usage

Cluster_Gauss_Newton_EBE_method(
  CGNM_result,
  nonlinearFunction,
  individualIndices_vec,
  numRepeat = 1,
  keepInitialDistribution = NA,
  algorithmParameter_EBEweight = 1,
  outputS4 = FALSE,
  ...
)

Arguments

Value

list of a matrix X, Y,residual_history, initialX, bootstrapX, bootstrapY as well as a list runSetting.

1. X, Y, residual_history, initialX: identical to what was given as CGNM_result.

2. X: a num_bootstrapSample by n matrix which stores the the X values that was sampled using residual resampling bootstrap analyses (In terms of model fitting this is the parameter combinations with variabilities that represent parameter estimation uncertainties.).

3. Y: a num_bootstrapSample by m matrix which stores the nonlinearFunction evaluated at the corresponding bootstrap analyses results in matrix bootstrapX above. In the context of model fitting each row corresponds to the model simulations.

4. runSetting: identical to what is given as CGNM_result but in addition including num_bootstrapSample and indicesToUseAsInitialIterates.

Examples

##lip-flop kinetics (an example known to have two distinct solutions)

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation, num_iteration = 10, num_minimizersToFind = 100,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
lowerBound=rep(0,3), ParameterNames=c("Ka","V1","CL_2"), saveLog = FALSE)

CGNM_EBE=Cluster_Gauss_Newton_EBE_method(CGNM_result,
     nonlinearFunction=model_analytic_function, individualIndices_vec=seq(1,9))

Cluster_Gauss_Newton_method

Description

Find multiple minimisers of the nonlinear least squares problem.

argmin_x ||f(x)-y*||

where

1. f: nonlinear function (e.g., mathematical model)

2. y*: target vector (e.g., observed data to fit the mathematical model)

3. x: variable of the nonlinear function that we aim to find the values that minimize (minimizers) the differences between the nonlinear function and target vector (e.g., model parameter)

Parameter estimation problems of mathematical models can often be formulated as nonlinear least squares problems. In this context f can be thought at a model, x is the parameter, and y* is the observation. CGNM iteratively estimates the minimizer of the nonlinear least squares problem from various initial estimates hence finds multiple minimizers. Full detail of the algorithm and comparison with conventional method is available in the following publication, also please cite this publication when this algorithm is used in your research: Aoki et al. (2020) <doi.org/10.1007/s11081-020-09571-2>. Cluster Gauss–Newton method. Optimization and Engineering, 1-31. As illustrated in this paper, CGNM is faster and more robust compared to repeatedly applying the conventional optimization/nonlinear least squares algorithm from various initial estimates. In addition, CGNM can realize this speed assuming the nonlinear function to be a black-box function (e.g. does not use things like adjoint equation of a system of ODE as the function does not have to be based on a system of ODEs.).

Usage

Cluster_Gauss_Newton_method(
  nonlinearFunction,
  targetVector,
  initial_lowerRange,
  initial_upperRange,
  lowerBound = NA,
  upperBound = NA,
  ParameterNames = NA,
  stayIn_initialRange = FALSE,
  num_minimizersToFind = 250,
  num_iteration = 25,
  saveLog = TRUE,
  runName = "",
  textMemo = "",
  algorithmParameter_initialLambda = 1,
  algorithmParameter_gamma = 2,
  algorithmVersion = 3,
  initialIterateMatrix = NA,
  targetMatrix = NA,
  weightMatrix = NA,
  keepInitialDistribution = NA,
  MO_weights = NA,
  MO_values = NA,
  outputS4 = FALSE,
  varianceExponent = 0,
  ...
)

Arguments

Value

list of a matrix X, Y,residual_history and initialX, as well as a list runSetting

1. X: a num_minimizersToFind by n matrix which stores the approximate minimizers of the nonlinear least squares in each row. In the context of model fitting they are the estimated parameter sets.

2. Y: a num_minimizersToFind by m matrix which stores the nonlinearFunction evaluated at the corresponding approximate minimizers in matrix X above. In the context of model fitting each row corresponds to the model simulations.

3. residual_history: a num_minimizersToFind by (num_iteration+1) matrix storing the sum of squares residual of each approximate minimizer (rows) at each iteration including the initial iterate (columns).

4. initialX: a num_minimizersToFind by n matrix which stores the set of initial iterates.

5. runSetting: a list containing all the input variables to Cluster_Gauss_Newton_method (i.e., nonlinearFunction, targetVector, initial_lowerRange, initial_upperRange ,algorithmParameter_initialLambda, algorithmParameter_gamma, num_minimizersToFind, num_iteration, saveLog, runName, textMemo).

Examples

## Example 1 (start here): flip-flop kinetics defined as a plain R function
## (an example known to have two distinct solutions). This is the minimal,
## recommended pattern for calling Cluster_Gauss_Newton_method().

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation, num_iteration = 10, num_minimizersToFind = 100,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
saveLog = FALSE)

acceptedApproximateMinimizers(CGNM_result)

## Example 2 (optional, advanced): the same flip-flop kinetics model defined
## as an ODE system via the RxODE package, to illustrate fitting simulator-
## based models. Requires the RxODE package (not a dependency of CGNM), so
## this example is not run automatically.
## Not run: 
library(RxODE)

model_text="
d/dt(X_1)=-ka*X_1
d/dt(C_2)=(ka*X_1-CL_2*C_2)/V1"

model=RxODE(model_text)
#define nonlinearFunction
model_function=function(x){

observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)

theta <- c(ka=x[1],V1=x[2],CL_2=x[3])
ev <- eventTable()
ev$add.dosing(dose = 1000, start.time =0)
ev$add.sampling(observation_time)
odeSol=model$solve(theta, ev)
log10(odeSol[,"C_2"])

}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(nonlinearFunction=model_function,
targetVector = observation, saveLog = FALSE,
initial_lowerRange = c(0.1,0.1,0.1),initial_upperRange =  c(10,10,10))
## End(Not run)

acceptedApproximateMinimizers

Description

CGNM find multiple sets of minimizers of the nonlinear least squares (nls) problem by solving nls from various initial iterates. Although CGNM is shown to be robust compared to other conventional multi-start algorithms, not all initial iterates minimizes successfully. By assuming sum of squares residual (SSR) follows the chai-square distribution we first reject the approximated minimiser who SSR is statistically significantly worse than the minimum SSR found by the CGNM. Then use elbow-method (a heuristic often used in mathematical optimisation to balance the quality and the quantity of the solution found) to find the "acceptable" maximum SSR. This function outputs the acceptable approximate minimizers of the nonlinear least squares problem found by the CGNM.

Usage

acceptedApproximateMinimizers(
  CGNM_result,
  cutoff_pvalue = 0.05,
  numParametersIncluded = NA,
  useAcceptedApproximateMinimizers = TRUE,
  algorithm = 2,
  ParameterNames = NA,
  ReparameterizationDef = NA
)

Arguments

Value

A dataframe that each row stores the accepted approximate minimizers found by CGNM.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

acceptedApproximateMinimizers(CGNM_result)

acceptedIndices

Description

CGNM find multiple sets of minimizers of the nonlinear least squares (nls) problem by solving nls from various initial iterates. Although CGNM is shown to be robust compared to other conventional multi-start algorithms, not all initial iterates minimizes successfully. By assuming sum of squares residual (SSR) follows the chai-square distribution we first reject the approximated minimiser who SSR is statistically significantly worse than the minimum SSR found by the CGNM. Then use elbow-method (a heuristic often used in mathematical optimisation to balance the quality and the quantity of the solution found) to find the "acceptable" maximum SSR. This function outputs the indices of acceptable approximate minimizers of the nonlinear least squares problem found by the CGNM.

Usage

acceptedIndices(
  CGNM_result,
  cutoff_pvalue = 0.05,
  numParametersIncluded = NA,
  useAcceptedApproximateMinimizers = TRUE,
  algorithm = 2
)

Arguments

Value

A vector of natural number that contains the indices of accepted approximate minimizers found by CGNM.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

acceptedIndices(CGNM_result)

acceptedIndices_binary

Description

CGNM find multiple sets of minimizers of the nonlinear least squares (nls) problem by solving nls from various initial iterates. Although CGNM is shown to be robust compared to other conventional multi-start algorithms, not all initial iterates minimizes successfully. By assuming sum of squares residual (SSR) follows the chai-square distribution we first reject the approximated minimiser who SSR is statistically significantly worse than the minimum SSR found by the CGNM. Then use elbow-method (a heuristic often used in mathematical optimisation to balance the quality and the quantity of the solution found) to find the "acceptable" maximum SSR. This function outputs the indices of acceptable approximate minimizers of the nonlinear least squares problem found by the CGNM. (note that acceptedIndices(CGNM_result) is equal to seq(1,length(acceptedIndices_binary(CGNM_result)))[acceptedIndices_binary(CGNM_result)])

Usage

acceptedIndices_binary(
  CGNM_result,
  cutoff_pvalue = 0.05,
  numParametersIncluded = NA,
  useAcceptedApproximateMinimizers = TRUE,
  algorithm = 2
)

Arguments

Value

A vector of TRUE and FALSE that indicate if the each of the approximate minimizer found by CGNM is acceptable or not.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

acceptedIndices_binary(CGNM_result)

acceptedMaxSSR

Description

CGNM find multiple sets of minimizers of the nonlinear least squares (nls) problem by solving nls from various initial iterates. Although CGNM is shown to be robust compared to other conventional multi-start algorithms, not all initial iterates minimizes successfully. By assuming sum of squares residual (SSR) follows the chai-square distribution we first reject the approximated minimiser who SSR is statistically significantly worse than the minimum SSR found by the CGNM. Then use elbow-method (a heuristic often used in mathematical optimisation to balance the quality and the quantity of the solution found) to find the "acceptable" maximum SSR.

Usage

acceptedMaxSSR(
  CGNM_result,
  cutoff_pvalue = 0.05,
  numParametersIncluded = NA,
  useAcceptedApproximateMinimizers = TRUE,
  algorithm = 2
)

Arguments

Value

A positive real number that is the maximum sum of squares residual (SSR) the algorithm has selected to accept.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

acceptedMaxSSR(CGNM_result)

Convert a classic CGNM result list into a CGNM_result S4 object

Description

Cluster_Gauss_Newton_method(), Cluster_Gauss_Newton_Bootstrap_method(), and Cluster_Gauss_Newton_EBE_method() all accept an outputS4 = TRUE argument that does this conversion automatically. Use as_CGNM_result_S4() directly when you have an existing classic list result (for example one loaded from an older saved .RDATA file) that you want to convert after the fact. Calling it on an object that is already a CGNM_result S4 object returns that object unchanged.

Usage

as_CGNM_result_S4(CGNM_result_list)

Arguments

Value

a CGNM_result S4 object exposing the same fields via $/[[.

Examples

model_function <- function(x) x
CGNM_result <- Cluster_Gauss_Newton_method(
  nonlinearFunction = model_function, targetVector = c(1, 2, 3),
  initial_lowerRange = rep(0.1, 3), initial_upperRange = rep(10, 3),
  num_iteration = 2, num_minimizersToFind = 5, saveLog = FALSE)
CGNM_result_S4 <- as_CGNM_result_S4(CGNM_result)
CGNM_result_S4$X

bestApproximateMinimizers

Description

Returns the approximate minimizers with minimum SSR found by CGNM.

Usage

bestApproximateMinimizers(
  CGNM_result,
  numParameterSet = 1,
  ParameterNames = NA,
  ReparameterizationDef = NA
)

Arguments

Value

A vector a vector of accepted approximate minimizers with minimum SSR found by CGNM.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

bestApproximateMinimizers(CGNM_result,10)

col_quantile

Description

Obtain columb wise quantile

Usage

col_quantile(data_in, prob)

Arguments

Value

a vector of number dim(data_in)[2] containing: quantile of the each column where the probability is specified as "prob"

Examples

A=matrix(seq(1,100),nrow = 25)
col_quantile(A, 0.5)

compare_profileLikelihood

Description

Draw profile likelihood surface using the function evaluations conducted during CGNM computation. Note plot_SSRsurface can only be used when log is saved by setting saveLog=TRUE option when running Cluster_Gauss_Newton_method(). The grey horizontal line is the threshold for 95% pointwise confidence interval.

Usage

compare_profileLikelihood(
  logLocation,
  alpha = 0.25,
  numBins = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  showInitialRange = TRUE,
  Likelihood_function = Residual_function_def
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

plot_profileLikelihood("CGNM_log")

## End(Not run)

generateCGNM_script

Description

Programmatic (non-GUI) equivalent of the shinyCGNM app's "make CGNM script" workflow. Takes the ODE model text and the same three tables the shinyCGNM app uses (parameter info, observed data, and dose data), runs the same validation checks the app performs on those tables (ODE compiles, dosing targets a real compartment, every referenced parameter has a range, observed values are numeric, residual error model is 0/1, ...), and returns the generated R script (an rxode2-based model function plus a Cluster_Gauss_Newton_method/Cluster_Gauss_Newton_Bootstrap_method call) as a character string, ready to be written to a file, cat()'d, or eval(parse(text = .))'d directly.

Usage

generateCGNM_script(
  ODE_text,
  parameterInfo_table,
  observedData_table,
  doseData_table = NA,
  initialConditionData_table = NA,
  runName = "",
  num_minimizersToFind = 250,
  num_iteration = 25,
  bootstrap = TRUE,
  parallel = "none",
  includeMiddleOutCode = FALSE,
  includeSimulationCode = FALSE,
  simulationDoseData_table = NA,
  simulationTimepoints_table = NA,
  scriptFileName = NA
)

Arguments

Value

string the generated R script.

Examples

## Not run: 
ODE_text="
d/dt(depot) = -ka*depot
d/dt(central) = ka*depot - (CL/V1)*central
C_central = central/V1
"

parameterInfo_table=make_ShinyCGNM_parameterInfo(
  ParameterName=c("ka","CL","V1"),
  Initial_lower_range=c(0.01,0.01,0.01),
  Initial_upper_range=c(100,100,100)
)

doseData_table=make_ShinyCGNM_doseData(
  ID="1", dose=1000, dosing.to="depot", start.time=0, rate=NA,
  nbr.doses=1, dosing.interval=NA
)

observedData_table=make_ShinyCGNM_observationData(
  ID="1",
  time=c(0.1, 0.2, 0.4, 0.6, 1, 2, 3, 6, 12),
  Observation_expression="log10(C_central)",
  Observed_value=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238)),
  ResidualError_model=0
)

script_text=generateCGNM_script(
  ODE_text=ODE_text,
  parameterInfo_table=parameterInfo_table,
  observedData_table=observedData_table,
  doseData_table=doseData_table,
  runName="oral1cpt"
)

cat(script_text)
eval(parse(text=script_text))

## Same model, but with an explicit initial condition instead of (or in addition
## to) a dose: e.g. central starts at 5 instead of implicitly at 0.
initialConditionData_table=make_ShinyCGNM_initialCondition(
  ID="1", state="central", value=5
)
script_text_withInit=generateCGNM_script(
  ODE_text=ODE_text,
  parameterInfo_table=parameterInfo_table,
  observedData_table=observedData_table,
  doseData_table=doseData_table,
  initialConditionData_table=initialConditionData_table,
  runName="oral1cpt_withInit"
)

## Same fit, plus a post-hoc middle-out template and a simulation of a denser time
## grid at a higher dose (bootstrap = TRUE is required for the simulation part).
## note: unlike observedData_table, Observation_expression here must be an exact
## compartment/derived-variable name from ODE_text (e.g. "C_central"), not an
## expression like "log10(C_central)" - it is not evaluated, only matched by name.
simulationTimepoints_table=make_ShinyCGNM_simulationTimepoints(
  ID="1", time=seq(0.1,12,by=0.1), Observation_expression="C_central"
)
simulationDoseData_table=make_ShinyCGNM_doseData(
  ID="1", dose=2000, dosing.to="depot", start.time=0, rate=NA,
  nbr.doses=1, dosing.interval=NA
)
script_text_full=generateCGNM_script(
  ODE_text=ODE_text,
  parameterInfo_table=parameterInfo_table,
  observedData_table=observedData_table,
  doseData_table=doseData_table,
  runName="oral1cpt_full",
  bootstrap=TRUE,
  includeMiddleOutCode=TRUE,
  includeSimulationCode=TRUE,
  simulationDoseData_table=simulationDoseData_table,
  simulationTimepoints_table=simulationTimepoints_table
)

## End(Not run)

make_ShinyCGNM_doseData

Description

A helper function to write out the csv file that can be read in as the dose file in shinyCGNM

Usage

make_ShinyCGNM_doseData(
  ID,
  dose,
  dosing.to,
  start.time,
  rate = NA,
  nbr.doses = 1,
  dosing.interval = NA,
  fileName = NA
)

Arguments

Value

data.frame if fileName is NA Null if fileName is not NA but instead write out the csv file

Examples

make_ShinyCGNM_doseData(
ID=seq(1,5),
dose=10,
dosing.to="A_admin",
start.time=0,
rate=NA
)

make_ShinyCGNM_initialCondition

Description

A helper function to build the (optional) initial condition table accepted by generateCGNM_script's initialConditionData_table. This lets you set the value of an ODE compartment at time 0 explicitly, independently of (or together with) any dosing into that compartment from a dose table.

Usage

make_ShinyCGNM_initialCondition(ID, state, value, fileName = NA)

Arguments

Value

data.frame if fileName is NA Null if fileName is not NA but instead write out the csv file

Examples

make_ShinyCGNM_initialCondition(
ID="1",
state="central",
value=5
)

make_ShinyCGNM_observationData

Description

A helper function to write out the csv file that can be read in as the observation file in shinyCGNM

Usage

make_ShinyCGNM_observationData(
  ID,
  time,
  Observation_expression,
  Observed_value,
  ResidualError_model,
  Memo = NA,
  fileName = NA
)

Arguments

Value

data.frame if fileName is NA Null if fileName is not NA but instead write out the csv file

Examples

make_ShinyCGNM_observationData(
ID=1,
time=c(1,2,3,6,12,24),
Observation_expression="C_central",
Observed_value=c(0.1, 0.3, 0.6, 0.1, 0.05, 0.01),
ResidualError_model=1
)

make_ShinyCGNM_parameterInfo

Description

A helper function to write out the csv file that can be read in as the parameter info file in shinyCGNM

Usage

make_ShinyCGNM_parameterInfo(
  ParameterName,
  Initial_lower_range,
  Initial_upper_range,
  Lower_bound = 0,
  Upper_bound = NA,
  VaryByID = 0,
  MO_weight = 0,
  MO_value = NA,
  Unit = NA,
  fileName = NA
)

Arguments

Value

data.frame if fileName is NA Null if fileName is not NA but instead write out the csv file

Examples

make_ShinyCGNM_parameterInfo(
ParameterName=c("ka","CL","V1"),
Initial_lower_range=c(0.01,0.01,0.01),
Initial_upper_range=c(100,100,100)
)

make_ShinyCGNM_simulationTimepoints

Description

A helper function to build the simulation time point table accepted by generateCGNM_script's simulationTimepoints_table (used when includeSimulationCode = TRUE). Unlike make_ShinyCGNM_observationData, there is no Observed_value/ResidualError_model, since nothing is being fit here - these are just the time points and expressions to simulate and plot.

Usage

make_ShinyCGNM_simulationTimepoints(
  ID,
  time,
  Observation_expression,
  Memo = NA,
  fileName = NA
)

Arguments

Value

data.frame if fileName is NA Null if fileName is not NA but instead write out the csv file

Examples

make_ShinyCGNM_simulationTimepoints(
ID=1,
time=seq(0,24,by=0.5),
Observation_expression="C_central"
)

plot_2DprofileLikelihood

Description

Make likelihood related values v.s. parameterValues plot using the function evaluations used during CGNM computation. Note plot_SSRsurface can only be used when log is saved by setting saveLog=TRUE option when running Cluster_Gauss_Newton_method().

Usage

plot_2DprofileLikelihood(
  logLocation,
  index_x = NA,
  index_y = NA,
  plotType = 2,
  plotMax = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  numBins = NA,
  showInitialRange = TRUE,
  alpha = 0.25,
  Likelihood_function = Residual_function_def,
  dependentVariable = NA
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=10^x[1]
 V1=10^x[2]
 CL_2=10^x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(-1,-1,-1), initial_upperRange =  c(1,1,1),
num_iter = 10, num_minimizersToFind = 500, saveLog=TRUE)

## the minimum example
plot_2DprofileLikelihood("CGNM_log")

## we can draw profilelikelihood also including bootstrap result
CGNM_result=Cluster_Gauss_Newton_Bootstrap_method(CGNM_result,
                      nonlinearFunction = model_analytic_function)

## example with various options
plot_2DprofileLikelihood(c("CGNM_log","CGNM_log_bootstrap"),
 showInitialRange = TRUE,index_x = c("ka","V1"))
 
## End(Not run)

plot_Rank_SSR

Description

Make SSR v.s. rank plot. This plot is often used to visualize the maximum accepted SSR.

Usage

plot_Rank_SSR(CGNM_result, indicesToInclude = NA)

Arguments

Value

A ggplot object of SSR v.s. rank.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_Rank_SSR(CGNM_result)

plot_SSR_parameterValue

Description

Make SSR v.s. parameterValue plot of the accepted approximate minimizers found by the CGNM. Bars in the violin plots indicates the interquartile range.

Usage

plot_SSR_parameterValue(
  CGNM_result,
  indicesToInclude = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  showInitialRange = TRUE
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_SSR_parameterValue(CGNM_result)

plot_SSRsurface

Description

Make minimum SSR v.s. parameterValue plot using the function evaluations used during CGNM computation. Note plot_SSRsurface can only be used when log is saved by setting saveLog=TRUE option when running Cluster_Gauss_Newton_method().

Usage

plot_SSRsurface(
  logLocation,
  alpha = 0.25,
  profile_likelihood = FALSE,
  numBins = NA,
  maxSSR = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  showInitialRange = FALSE,
  Residual_function = Residual_function_def,
  compareBetweenModels = FALSE
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

plot_SSRsurface("CGNM_log") + scale_y_continuous(trans='log10')

## End(Not run)

plot_goodnessOfFit

Description

Make goodness of fit plots to assess the model-fit and bias in residual distribution. The linear model is fit to the residual and plotted using geom_smooth(method=lm) in ggplot.

Explanation of the terminologies in terms of PBPK model fitting to the time-course drug concentration measurements:
"independent variable" is time
"dependent variable" is the concentration.
"Residual" is the difference between the measured concentration and the model simulation with the parameter fond by the CGNM.
"m" is number of observations

Usage

plot_goodnessOfFit(
  CGNM_result,
  plotType = 1,
  plotRank = c(1),
  independentVariableVector = NA,
  dependentVariableTypeVector = NA,
  absResidual = FALSE
)

Arguments

Value

A ggplot object of the goodness of fit plot.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=10^x[1]
 V1=10^x[2]
 CL_2=10^x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = rep(0.01,3), initial_upperRange =  rep(100,3),
lowerBound=rep(0,3), ParameterNames = c("Ka","V1","CL"),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_goodnessOfFit(CGNM_result)
plot_goodnessOfFit(CGNM_result,
     independentVariableVector=c(0.1,0.2,0.4,0.6,1,2,3,6,12))

plot_paraDistribution_byHistogram

Description

Make histograms to visualize the initial distribution and distribition of the accepted approximate minimizers found by the CGNM.

Usage

plot_paraDistribution_byHistogram(
  CGNM_result,
  indicesToInclude = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  bins = 30
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))



CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = rep(0.01,3), initial_upperRange =  rep(100,3),
lowerBound=rep(0,3), ParameterNames = c("Ka","V1","CL"),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_paraDistribution_byHistogram(CGNM_result)
plot_paraDistribution_byHistogram(CGNM_result,
     ReparameterizationDef=c("log10(Ka)","log10(V1)","log10(CL)"))

plot_paraDistribution_byViolinPlots

Description

Make violin plot to compare the initial distribution and distribition of the accepted approximate minimizers found by the CGNM. Bars in the violin plots indicates the interquartile range. The solid line connects the interquartile ranges of the initial distribution and the distribution of the accepted approximate minimizer at the final iterate. The blacklines connets the minimums and maximums of the initial distribution and the distribution of the accepted approximate minimizer at the final iterate. The black dots indicate the median.

Usage

plot_paraDistribution_byViolinPlots(
  CGNM_result,
  indicesToInclude = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))



CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = rep(0.01,3), initial_upperRange =  rep(100,3),
lowerBound=rep(0,3), ParameterNames = c("Ka","V1","CL"),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_paraDistribution_byViolinPlots(CGNM_result)
plot_paraDistribution_byViolinPlots(CGNM_result,
     ReparameterizationDef=c("log10(Ka)","log10(V1)","log10(CL)"))

plot_parameterValue_scatterPlots

Description

Make scatter plots of the accepted approximate minimizers found by the CGNM. Bars in the violin plots indicates the interquartile range.

Usage

plot_parameterValue_scatterPlots(CGNM_result, indicesToInclude = NA)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

plot_parameterValue_scatterPlots(CGNM_result)

plot_profileLikelihood

Description

Draw profile likelihood surface using the function evaluations conducted during CGNM computation. Note plot_SSRsurface can only be used when log is saved by setting saveLog=TRUE option when running Cluster_Gauss_Newton_method(). The grey horizontal line is the threshold for 95% pointwise confidence interval.

Usage

plot_profileLikelihood(
  logLocation,
  alpha = 0.25,
  numBins = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  showInitialRange = TRUE,
  Likelihood_function = Residual_function_def
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

plot_profileLikelihood("CGNM_log")

## End(Not run)

plot_simulationMatrixWithCI

Description

Plot simulation that are provided to plot confidence interval (or more like a confidence region).

Usage

plot_simulationMatrixWithCI(
  simulationMatrix,
  independentVariableVector = NA,
  dependentVariableTypeVector = NA,
  confidenceLevels = c(0.25, 0.75),
  observationVector = NA,
  observationIndpendentVariableVector = NA,
  observationDependentVariableTypeVector = NA,
  overLay = FALSE
)

Arguments

Value

A ggplot object.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 (Cp)
}

observation=(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation, num_iteration = 10, num_minimizersToFind = 100,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
lowerBound=rep(0,3), ParameterNames=c("Ka","V1","CL_2"), saveLog = FALSE)

CGNM_bootstrap=Cluster_Gauss_Newton_Bootstrap_method(CGNM_result,
     nonlinearFunction=model_analytic_function, num_bootstrapSample=100)


plot_simulationMatrixWithCI(CGNM_result$bootstrapY,
independentVariableVector=observation_time, observationVector=observation)

## End(Not run)

plot_simulationWithCI

Description

Plot model simulation where the various parameter combinations are provided and conduct simulations and then the confidence interval (or more like a confidence region) is plotted.

Usage

plot_simulationWithCI(
  simulationFunction,
  parameter_matrix,
  independentVariableVector = NA,
  dependentVariableTypeVector = NA,
  confidenceLevels = c(0.25, 0.75),
  observationVector = NA,
  observationIndpendentVariableVector = NA,
  observationDependentVariableTypeVector = NA,
  overLay = FALSE
)

Arguments

Value

A list including ggplot object ($plot) and simulated data matrix ($plotData_matrix).

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 (Cp)
}

observation=(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation, num_iteration = 10, num_minimizersToFind = 100,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
lowerBound=rep(0,3), ParameterNames=c("Ka","V1","CL_2"), saveLog = FALSE)

CGNM_bootstrap=Cluster_Gauss_Newton_Bootstrap_method(CGNM_result,
     nonlinearFunction=model_analytic_function, num_bootstrapSample=100)

plot_simulationWithCI(model_analytic_function, as.matrix(CGNM_result$bootstrapTheta),
independentVariableVector=observation_time, observationVector=observation)

## End(Not run)

shinyCGNM

Description

Start Shiny app that assist the user to make R-script to conduct PBPK model fitting using CGNM.

Usage

shinyCGNM()

Value

NULL and start graphifcal user interface.

Examples

## Not run: 
shinyCGNM()

## End(Not run)

suggestInitialLowerRange

Description

Suggest initial lower range based on the profile likelihood. The user can re-run CGNM with this suggested initial range so that to improve the convergence.

Usage

suggestInitialLowerRange(logLocation, alpha = 0.25, numBins = NA)

Arguments

Value

A numerical vector of suggested initial lower range based on profile likelihood.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

suggestInitialLowerRange("CGNM_log")

## End(Not run)

suggestInitialUpperRange

Description

Suggest initial upper range based on the profile likelihood. The user can re-run CGNM with this suggested initial range so that to improve the convergence.

Usage

suggestInitialUpperRange(logLocation, alpha = 0.25, numBins = NA)

Arguments

Value

A numerical vector of suggested initial upper range based on profile likelihood.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

suggestInitialLowerRange("CGNM_log")

## End(Not run)

table_parameterSummary

Description

Make summary table of the approximate local minimizers found by CGNM. If bootstrap analysis result is available, relative standard error (RSE: standard deviation/mean) will also be included in the table.

Usage

table_parameterSummary(
  CGNM_result,
  indicesToInclude = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  pretty = FALSE
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))



CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = rep(0.01,3), initial_upperRange =  rep(100,3),
lowerBound=rep(0,3), ParameterNames = c("Ka","V1","CL"),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

table_parameterSummary(CGNM_result)
table_parameterSummary(CGNM_result,
     ReparameterizationDef=c("log10(Ka)","log10(V1)","log10(CL)"))
table_parameterSummary(CGNM_result, pretty=TRUE)

table_profileLikelihoodConfidenceInterval

Description

Make table of confidence intervals that are approximated from the profile likelihood. First inspect profile likelihood plot and make sure the plot is smooth and has good enough resolution and the initial range is appropriate. Do not report this table without checking the profile likelihood plot.

Usage

table_profileLikelihoodConfidenceInterval(
  logLocation,
  alpha = 0.25,
  numBins = NA,
  ParameterNames = NA,
  ReparameterizationDef = NA,
  pretty = FALSE,
  Likelihood_function = Residual_function_def,
  silent = FALSE
)

Arguments

Value

A ggplot object including the violin plot, interquartile range and median, minimum and maximum.

Examples

## Not run: 
model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog=TRUE)

table_profileLikelihoodConfidenceInterval("CGNM_log")

## End(Not run)

topIndices

Description

CGNM find multiple sets of minimizers of the nonlinear least squares (nls) problem by solving nls from various initial iterates. Although CGNM is shown to be robust compared to other conventional multi-start algorithms, not all initial iterates minimizes successfully. One can visually inspect rank v.s. SSR plot and manually choose number of best fit acceptable parameters. By using this function "topIndices", we can obtain the indices of the "numTopIndices" best fit parameter combinations.

Usage

topIndices(CGNM_result, numTopIndices)

Arguments

Value

A vector of natural number that contains the indices of accepted approximate minimizers found by CGNM.

Examples

model_analytic_function=function(x){

 observation_time=c(0.1,0.2,0.4,0.6,1,2,3,6,12)
 Dose=1000
 F=1

 ka=x[1]
 V1=x[2]
 CL_2=x[3]
 t=observation_time

 Cp=ka*F*Dose/(V1*(ka-CL_2/V1))*(exp(-CL_2/V1*t)-exp(-ka*t))

 log10(Cp)
}

observation=log10(c(4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238))

CGNM_result=Cluster_Gauss_Newton_method(
nonlinearFunction=model_analytic_function,
targetVector = observation,
initial_lowerRange = c(0.1,0.1,0.1), initial_upperRange =  c(10,10,10),
num_iter = 10, num_minimizersToFind = 100, saveLog = FALSE)

topInd=topIndices(CGNM_result, 10)

## This gives top 10 approximate minimizers
CGNM_result$X[topInd,]