Minimization Tool for Pharmacokinetic-Pharmacodynamic Data Analysis

Description

This is a minimization tool to solve the examples in the book Gabrielsson J, Weiner D. 'Pharmacokinetic and Pharmacodynamic Data Analysis - Concepts and Applications' 5th ed. 2016. (ISBN:9198299107).

Details

This is a set of minimization tools to solve all the examples in the book 'Pharmacokinetic and Pharmacodynamic Data Analysis - Concepts and Applications' 5th ed. 2016.

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Gabrielsson J, Weiner D. Pharmacokinetic and Pharmacodynamic Data Analysis - Concepts and Applications. 5th ed. 2016.

Examples

tData = Theoph
colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV")

fPK = function(THETA)     # Prediction function
{
  DOSE = 320000           # in microgram
  TIME = e$DATA[,"TIME"]  # use data in e$DATA

  K    = THETA[1]
  Ka   = THETA[2]
  V    = THETA[3]

  Cp   = DOSE/V*Ka/(Ka - K)*(exp(-K*TIME) - exp(-Ka*TIME))
  return(Cp)
}

IDs = unique(tData[,"ID"])
nID = length(IDs)
for (i in 1:nID) {
  Data = tData[tData$ID == IDs[i],]
  Res = nlr(fPK, Data, pNames=c("k", "ka", "V"), IE=c(0.1, 3, 500), 
            SecNames=c("CL", "Thalf", "MRT"), SecForms=c(~V*k, ~log(2)/k, ~1/k))
  print(paste("## ID =", i, "##"))
  print(Res)
}

Internal Functions

Description

Internal functions

Details

These are not to be called by the user.

One compartment model - analytical

Description

It calculates using one compartment model.

Usage

Comp1(Ke, Ka=0, DH)

Arguments

Details

First compartment is the gut compartment for oral dosing. IV bolus and infusion dosing should be done at the second compartment.

Value

This returns a table with the gut and the central compartment columns

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DAT
DAT2 = ExpandDH(DAT)
X1 = Comp1(Ke=0.1, Ka=1, DAT2)
X1
matplot(DAT2[, "TIME"], X1, type="l")

An Example of Dosing History Table

Description

This is a conventional NONMEM input data format.

Usage

DAT

Format

This data frame has 5 columns with 18 time-points for the simulation.

TIME

Time

AMT

Amount given for the compartment of CMT column

RATE

Infusion rate

CMT

Compartment number, 1=gut, 2=central, 3=peripheral, etc.

DV

Currently blank and not used.

Details

To be used at Comp1 or nComp, expand dosing history with ExpandDH function.

Environment's Objects

Description

Get an environment's visible objects as a list.

Usage

  EnvObj(envir = e)

Arguments

Details

All the visible objects in the environment including functions and data will be returned.

Value

All visible objects as a list

Author(s)

Kyun-Seop Bae <k@acr.kr>

Expand Dosing History Table

Description

It expands dosing history table.

Usage

ExpandDH(DH, Fo = 1)

Arguments

Details

It expands dosing history table of conventional NONMEM data format. It calculate bioavailable amount, then add time points of non-differentiable, e.g. stopping points of infusion.

Value

Returns expanded dosing history table.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DAT
ExpandDH(DAT) # One observation point is increased at the time of 27.

Internal Obj Functions

Description

Internal Obj functions

Details

These are not to be called by the user.

Get Secondary Parameter Estimates

Description

Get standard error and relative standard error (cv) of the secondary parameter estimate

Usage

Secondary(Formula, PE, COV)

Arguments

Details

Variables within Formula should exist in the names of PE vector.

Value

This returns point estimate, standard error, relative standard error of the secondary parameter estimate.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

tData = Theoph
colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV") # Table requires DV column

fPK = function(THETA) # Prediction function
{
  AMT  = 320000 # in microgram
  TIME = e$DATA[,"TIME"]
  V  = THETA[1]
  K  = THETA[2]
  Ka = THETA[3]
  Cp = AMT/V*Ka/(Ka - K)*(exp(-K*TIME) - exp(-Ka*TIME))
  return(Cp)
}
Data = tData[tData$ID == 1,]
Res = nlr(fPK, Data, pNames=c("V", "K", "Ka"), IE=c(30000, 0.1, 2))
Secondary(~V*K, Res$Est["PE",1:e$nPara], Res$Cov)

Get Lambdas and Coefficients of Two-compartment Model

Description

It calculates lambdas and coefficients for two-compartment model from K10, K12, and K21.

Usage

SolComp2(K10, K12, K21)

Arguments

Details

It calculates lambdas and coefficients of two-compartment model from K10, K12, and K21. Lambdas should have no identical values.

Value

This returns a list of lambdas and coefficients.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DAT
DAT2 = ExpandDH(DAT)
Sol = SolComp2(K10=0.1, K12=3, K21=1)
X2 = nComp(Sol, Ka=1, DAT2)
X2
matplot(DAT2[, "TIME"], X2, type="l")

Get Lambdas and Coefficients of Three-compartment Model

Description

It calculates lambdas and coefficients for three-compartment model from K10, K12, K21, K13, and K31.

Usage

SolComp3(K10, K12, K21, K13, K31)

Arguments

Details

It calculates lambdas and coefficients of two-compartment model from K10, K12, and K21. Lambdas should have no identical values.

Value

This returns a list of lambdas and coefficients.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DAT
DAT2 = ExpandDH(DAT)
Sol = SolComp3(K10=0.1, K12=3, K21=1, K13=2, K31=0.5)
X3 = nComp(Sol, Ka=1, DAT2)
X3
matplot(DAT2[, "TIME"], X3, type="l")

Compare model with Chi-square test

Description

It performs chi-square test for two models comparison.

Usage

cmpChi(r1, r2)

Arguments

Details

One model should include the other model.

Value

Returns a p-value from pchisq

Author(s)

Kyun-Seop Bae <k@acr.kr>

Simplest diagnostic plot for minimization result

Description

It performs a simple diagnostic plot from the result of nlr.

Usage

dx(r)

Arguments

Details

This plots 'Observation vs. Prediction' and 'Normalized Redisual vs. Prediction' only. Normalized residual are meant to be distributed as standard normal distribution, N(0, 1).

Value

This just draws a plot.

Author(s)

Kyun-Seop Bae <k@acr.kr>

environment for internal data

Description

This is for the storage of intermediate data.

Hougaard Measure of Skewness

Description

Hougaard measure of skewness with nonlinear regression

Usage

  hSkew(rx)

Arguments

Details

Hougaard measure of skewness can be used to check if the parameters of nonlinear regression behavior in linear fashion, i.e. symmetric confidence interval. Be cautious on the variable name conflict. All the variables in the nonlinear function should be able to be accessed by the function.

Value

Hougaard estimate of skewness for each parameter

Author(s)

Kyun-Seop Bae k@acr.kr

References

EL-Shehawy SA. On Calculating the Hougaard Measure of Skewness in a Nonlinear Regression Model with Two Parameters. J Math & Stat. 2009;5(4):360-364.

Examples

  r1 = nls(density ~ b1*conc/(b2 + conc), DNase[DNase$Run == 1, ], start=list(b1=3, b2=5))
  hSkew(r1)

Get Amounts of Each Compartments using Lambdas and Coefficients of Multi-compartment Model

Description

It calculates using multi-compartment model.

Usage

nComp(Sol, Ka=0, DH)

Arguments

Details

First compartment is the gut compartment for oral dosing. IV bolus and infusion dosing should be done at the second compartment. If a bolus dose was given at time T, it is reflected at times of larger than T. This is more close to real observation. ADAPT does like this, but NONMEM does not.

Value

This returns a table with the gut and the other compartment columns

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DAT
DAT2 = ExpandDH(DAT)
Sol = SolComp2(K10=0.1, K12=3, K21=1)
X2 = nComp(Sol, Ka=1, DAT2)
X2
matplot(DAT2[, "TIME"], X2, type="l")

Nonlinear Regression in R

Description

It performs nonlinear regression usually for pharmacokinetic and pharmacodynamic models.

Usage

nlr(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms, 
    Method="L-BFGS-B", Sx, conf.level=0.95, k, fix=0)

Arguments

Details

This uses scaled transformed parameters and environment e internally.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

  tData = Theoph
  colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV")

  fPK = function(THETA) # Prediction function
  {
    DOSE = 320000 # in microgram
    TIME = e$DATA[, "TIME"] # use data in e$DATA

    K    = THETA[1]
    Ka   = THETA[2]
    V    = THETA[3]

    P  = DOSE/V*Ka/(Ka - K) * (exp(-K*TIME) - exp(-Ka*TIME))
    return(P)
  }

  IDs = unique(tData[,"ID"])
  nID = length(IDs)
  for (i in 1:nID) {
    Data = tData[tData$ID == IDs[i],]
    Res = nlr(fPK, Data, pNames=c("k", "ka", "V"), IE=c(0.1, 3, 500), 
              SecNames=c("CL", "Thalf", "MRT"), SecForms=c(~V*k, ~log(2)/k, ~1/k))
    print(paste("## ID =", i, "##"))
    print(Res)
  }

  # Another example from radioimmunoassay(RIA)
  d1 = data.frame(conc = c(200, 100, 50, 25, 12.5, 6.25, 3.125, 0),
                  DV = c(1.78, 1.5, 1.17, 0.74, 0.51, 0.31, 0.19, 0.04))

  PRED = function(TH) TH[1] + TH[2]*d1$conc^TH[4]/(TH[3]^TH[4] + d1$conc^TH[4])
  Scale = function(TH) 1/(PRED(TH) - (TH[1] + TH[2])/2)^2

  nlr(PRED, d1, pNames=c("R0", "Rmax", "RC50", "Hill"), IE=c(0.1, 3, 50, 1), 
      Error="S", Sx=Scale)

Plot Compartment Model Diagram

Description

It plots the diagrom of a comparment model.

Usage

pComp(dComp, dRate, Shape="rect", Col=NA, Bx=0.3, By=0.2, Cex=1.0, Lwd=3, 
      Radius=0.3, thIn=pi/2, thOut=pi/2, ...)

Arguments

Details

Flow direction is from the top to bottom.

Value

It plots.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

dA = data.frame(No = c(1, 2, 3, 4), Name=c("Gut Depot", "Skin Depot", "Central", "Peripheral"), 
                Level=c(1, 1, 2, 2), xPos=c(-0.5, 0.5, 0, 1))
dB = data.frame(From = c(1, 2, 3, 4, 3, 0, 0), To=c(3, 3, 4, 3, 5, 1, 2), 
                Name=c("KA", "KA2", "K12", "K21", "CL", "F1", "F2"))

pComp(dA, dB)
#par(oma=c(0, 0, 0, 0), mar=c(0, 0, 1, 0)) # If need, adjust margin before calling
pComp(dA, dB, "circ", main="Compartmental Model Diagram")

pComp(dA, dB, "circ", main="Compartmental Model Diagram", Col="#DDEEFF", asp=1)

Plot Likelihood or Objective Fnnction Value Profile

Description

It plots estimated likelihood profile. This is not profile likelihood profile.

Usage

  pProf(Bag = e, Title = "", ...)

Arguments

Details

This plots likelihood profile from the result of nlr() function. Bag should contain the results of nlr().

Value

No values but a plot.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Old type WinNonlin - Least Square not MLE

Description

It performs old type Winnonlin regression.

Usage

wnl5(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjLS)

Arguments

Details

This uses scaled transformed parameters and environment e internally. Here we do not provide standard error. If you want standard error, use nlr.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

tData = Theoph
colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV")

fPK = function(THETA) # Prediction function
{
  DOSE = 320000 # in microgram
  TIME = e$DATA[,"TIME"] # use data in e$DATA

  K  = THETA[1]
  Ka = THETA[2]
  V  = THETA[3]
  Cp  = DOSE/V*Ka/(Ka - K)*(exp(-K*TIME) - exp(-Ka*TIME))
  return(Cp)
}

IDs = unique(tData[,"ID"])
nID = length(IDs)
for (i in 1:nID) {
  Data = tData[tData$ID == IDs[i],]
  Res = wnl5(fPK, Data, pNames=c("k", "ka", "V"), IE=c(0.1, 3, 500))
  print(paste("## ID =", i, "##"))
  print(Res)
}