R/nca_utils.R
In IntiQuan/NonCompartIQR: Supportive functions for NCA
Defines functions
singleNCA_IQRnca
getSlope_IQRnca
LogAUC_IQRnca
LinAUC_IQRnca
Interpol_IQRnca
IntAUC_IQRnca
BestSlope_IQRnca
AUC_IQRnca
UT_IQRnca
UnitUrine_IQRnca
Unit_IQRnca
# Original Author: Kyun-Seop Bae k@acr.kr
# Adapted by: Henning Schmidt henning.schmidt@intiquan.com
# Based on NonCompart package on CRAN
########################################################################### ###
# singleNCA_IQRnca function
########################################################################### ###
singleNCA_IQRnca <- function(
x,
y,
dose=0,
adm="Extravascular",
dur=0,
doseUnit="mg",
timeUnit="h",
concUnit="ug/L",
NCAinfo=0, # Dummy for licensing check
iAUC="",
down="Log",
MW=0,
slopePoints = NULL,
NOSLOPE = NULL,
TOL = 1e-4,
FLAGsteadystate = FALSE,
tau = NULL,
IQRversion = NULL,
FLAGivSlopeCmax = FALSE) {
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.8.14
# Called by: tblNCA
# Calls: BestSlope, IntAUC, Unit, UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# dose: dose, actually given dose, not per kg or per square meter dose!
# adm: administration method, "Extravascular", "Bolus", or "Infusion"
# dur: duration of infusion
# doseUnit: dose unit (not per kg or per m2)
# timeUnit: time unit
# concUnit: concentration unit
# iAUC: data frame for interval AUC. See the example for the detail
# down: trapezoidal downward caculation. "Linear" or "Log"
# MW: molecular weight
# slopePoints: vector with indices of observations to use for slope calculation
# NOSLOPE: logical vector TRUE=>do not consider point for bestslope, FALSE: do consider
# IQRversion: version number for handling compatibility outside
# FLAGivSlopeCmax: if set to TRUE then Cmax obs will be allowed to be considered
# for bestslope calc. useful for 1 cpt dynamics and availability of very few points for NCA
# cannot be implemented by default as then not matching with other examples with richer data and
# comparison to WinNonLin
# RETURNS
# NCA result
# For steady-state data: Require presence of pre-dose sample at time=0 if infusion or extravascular administration
if (FLAGsteadystate) {
if (x[1] != 0) {
stop("Analysis of steady-state data requires availability of pre-dose sample\nin the data for infusion and extravascular administration")
}
if (is.na(y[1])) {
stop("Analysis of steady-state data requires availability of pre-dose sample (non-NA)\nin the data for infusion and extravascular administration")
}
# Check if last observation time point larger or equal to tau - otherwise error
if (tau > max(x)) {
stop("For steady-state profile analysis the last sample time needs to be larger than tau")
}
}
# Other checks
if (!(is.numeric(x) & is.numeric(y) & is.numeric(dose) & is.numeric(dur) & is.character(adm) & is.character(down))) {
stop("Check input types!")
}
if (UT_IQRnca(adm) == "INFUSION" & !all(dur > 0)) stop("Infusion mode should have dur larger than 0!")
###################################### ###
# Special Handling for slope calculation
###################################### ###
# Save original x and y values for determination of the right slope points after
# slope calculation (done via x1 to xorig__ matching). In addition xorig__ and yorig__
# are also used when manually points are selected for slope calculation
xorig__ <- x
yorig__ <- y
# Get xslope__, yslope__ to be used for bestslope algorithm
xslope__ <- x[!NOSLOPE] # Remove the points to be avoided in bestslope algorithm
yslope__ <- y[!NOSLOPE] # Remove the points to be avoided in bestslope algorithm
xslope__ = xslope__[yslope__ != 0] # Remove all points with zeros in y (including mid) for LAMZ
yslope__ = yslope__[yslope__ != 0] # Remove all points with zeros in y
NApoints__ <- is.na(xslope__) | is.na(yslope__)
xslope__ = xslope__[!NApoints__] # remove NA points in x
yslope__ = yslope__[!NApoints__] # remove NA points in y
###################################### ###
# Standard handling
###################################### ###
# Remove NA points
NApoints__ <- is.na(x) | is.na(y)
x <- x[!NApoints__] # remove NA points in x
y <- y[!NApoints__] # remove NA points in y
n <- length(x)
RetNames1__ = c("b0", "CMAX", "CMAXD", "TMAX", "TLAG", "CLST", "CLSTP", "TLST", "LAMZHL", "LAMZ",
"LAMZLL", "LAMZUL", "LAMZNPT", "CORRXY", "R2", "R2ADJ", "AUCLST", "AUCALL",
"AUCIFO", "AUCIFOD", "AUCIFP", "AUCIFPD", "AUCPEO", "AUCPEP",
"AUMCLST", "AUMCIFO", "AUMCIFP", "AUMCPEO", "AUMCPEP")
# If steady-state profile then add AUCTAU
if (FLAGsteadystate) {
RetNames1__ = c(RetNames1__, "AUCTAU","AUCTAUD","TAU")
}
if (UT_IQRnca(adm) == "BOLUS") {
RetNames1__ = union(RetNames1__, c("C0", "AUCPBEO", "AUCPBEP"))
}
if (UT_IQRnca(adm) == "EXTRAVASCULAR") {
# Some of them are dependent of SS and SD ... removed later what is not needed
RetNames1__ = union(RetNames1__, c("VZFO", "VZFP", "CLFO", "CLFP", "CLFSS", "VZFSS", "MRTEVLST", "MRTEVIFO", "MRTEVIFP"))
} else {
# Some of them are dependent of SS and SD ... removed later what is not needed
RetNames1__ = union(RetNames1__, c("VZO", "VZP", "CLO", "CLP", "CLSS", "VZSS", "MRTIVLST", "MRTIVIFO", "MRTIVIFP", "VSSO", "VSSP"))
}
Res__ = rep(NA_real_, length(RetNames1__))
names(Res__) = RetNames1__
if (n == 0 | n != length(y) | length(y[y < 0]) > 0) {
Res__["LAMZNPT"] = 0
return(Res__)
}
uY__ = unique(y)
if (length(uY__) == 1) { # Case of all the same values
Res__["CMAX"] = uY__
if (dose > 0) Res__["CMAXD"] = uY__ / dose
Res__["TMAX"] = x[y==uY__][1] # First Tmax
if (which(y==uY__)[1] > 1) {
Res__["TLAG"] = x[which(y==uY__) - 1]
} else {
Res__["TLAG"] = 0
}
Res__["CLST"] = uY__
Res__["TLST"] = x[y==uY__][1]
Res__["LAMZNPT"] = 0
Res__["b0"] = uY__
return(Res__)
}
iLastNonZero__ = max(which(y > 0)) # Index of last non-zero y
x0 = x[1:iLastNonZero__] # Till Non-zero concentration. i.e. removing trailing zeros
y0 = y[1:iLastNonZero__] # Till Non-zero concentration. i.e. removing trailing zeros
if (UT_IQRnca(adm) == "BOLUS") {
# Bolus - add an initial 0 at time 0 if not included in vector
# Works for steady-state and single dose data
if (is.na(x[x==0][1])) {
# Determine concentration at time 0
if (y[1] > y[2] & y[2] > 0) {
# Use extrapolation
C0 = exp(-x[1]*(log(y[2]) - log(y[1]))/(x[2] - x[1]) + log(y[1]))
} else {
# Use first non-zero value
C0 = y[x==min(x[y > 0])]
}
x2 = c(0, x)
y2 = c(C0, y)
x3 = c(0, x0)
y3 = c(C0, y0)
} else {
# 0 in x as first element. Check that y[1] > 0 otherwise error
if (y[1] <= 0) stop("For bolus administration the concentration at time 0 is not allowed to be 0")
x2 = x
y2 = y
x3 = x0
y3 = y0
# Set initial concentration at time 0 to observed one
C0 = y[1]
}
} else {
# Infusion & extravascular - add an initial 0 at time 0 if not included in vector
# Works for steady-state and single dose data
if (is.na(x[x==0][1])) {
# This case cannot happen for steady-state data as the 0 time point is required to be present and
# populated from the data.
x2 = c(0, x) # for AUCALL including trailing zero y values
y2 = c(0, y)
x3 = c(0, x0) # for AUCLST without trailing zero y values
y3 = c(0, y0)
} else {
# Time 0 included in x ... no need to add
# Works for single dose and steady-state dose!
x2 = x # for AUCALL including trailing zero y values
y2 = y
x3 = x0 # for AUCLST without trailing zero y values
y3 = y0
}
}
if (is.null(slopePoints)) {
tRes__ <- BestSlope_IQRnca(xslope__, yslope__, adm, TOL=TOL,FLAGivSlopeCmax)
# Return indices of points used for slope calculation
Indices_pointsSlope__ <- rev(seq(length(xslope__), by=-1, length.out = tRes__["LAMZNPT"]))
IXpointsSlope__ <- (1:length(x))[xorig__ %in% xslope__][Indices_pointsSlope__]
} else {
# Call the getSlope_IQRnca function
tRes__ <- getSlope_IQRnca(xorig__[slopePoints], log(yorig__[slopePoints]))
# Return indices of points used for slope calculation
IXpointsSlope__ <- slopePoints
}
if (length(tRes__) != 9) tRes__ = c(NA, NA, 0, NA, NA, NA, NA, NA, NA)
Res__[c("R2", "R2ADJ", "LAMZNPT", "LAMZ", "b0", "CORRXY", "LAMZLL", "LAMZUL", "CLSTP")] = tRes__
tabAUC = AUC_IQRnca(x3, y3, down)
Res__[c("AUCLST","AUMCLST")] = tabAUC[length(x3),]
Res__["AUCALL"] = AUC_IQRnca(x2, y2, down)[length(x2),1]
Res__["LAMZHL"] = log(2)/Res__["LAMZ"]
Res__["TMAX"] = x[which.max(y)][1]
Res__["CMAX"] = max(y)
Res__["TLST"] = x[iLastNonZero__]
Res__["CLST"] = y[iLastNonZero__]
Res__["AUCIFO"] = Res__["AUCLST"] + Res__["CLST"]/Res__["LAMZ"]
Res__["AUCIFP"] = Res__["AUCLST"] + Res__["CLSTP"]/Res__["LAMZ"]
Res__["AUCPEO"] = (1 - Res__["AUCLST"]/Res__["AUCIFO"])*100
Res__["AUCPEP"] = (1 - Res__["AUCLST"]/Res__["AUCIFP"])*100
Res__["AUMCIFO"] = Res__["AUMCLST"] + Res__["CLST"]*Res__["TLST"]/Res__["LAMZ"] + Res__["CLST"]/Res__["LAMZ"]/Res__["LAMZ"]
Res__["AUMCIFP"] = Res__["AUMCLST"] + Res__["CLSTP"]*Res__["TLST"]/Res__["LAMZ"] + Res__["CLSTP"]/Res__["LAMZ"]/Res__["LAMZ"]
Res__["AUMCPEO"] = (1 - Res__["AUMCLST"]/Res__["AUMCIFO"])*100
Res__["AUMCPEP"] = (1 - Res__["AUMCLST"]/Res__["AUMCIFP"])*100
if (!is.na(dose) & dose > 0) {
Res__["CMAXD"] = Res__["CMAX"]/dose
Res__["AUCIFOD"] = Res__["AUCIFO"]/dose
Res__["AUCIFPD"] = Res__["AUCIFP"]/dose
}
if (UT_IQRnca(adm) == "BOLUS") {
Res__["C0"] = C0 # Phoneix WinNonlin 6.4 User's Guide p27
Res__["AUCPBEO"] = tabAUC[2,1]/Res__["AUCIFO"]*100
Res__["AUCPBEP"] = tabAUC[2,1]/Res__["AUCIFP"]*100
} else {
if (sum(y0==0) > 0) Res__["TLAG"] = x0[max(which(y0==0))] # Trailing zero should not exist
else Res__["TLAG"] = 0
if (!is.na(x0[x0==0][1])) {
if (y0[x0==0] > 0) Res__["TLAG"] = 0 # This is WinNonlin logic
}
}
# Calculate AUCTAU using IntAUC_IQRnca for SS data
if (FLAGsteadystate) {
Res__["AUCTAU"] = IntAUC_IQRnca(x, y, 0, tau, Res__, down=down)
Res__["AUCTAUD"] = Res__["AUCTAU"]/dose
Res__["TAU"] = tau
}
if (UT_IQRnca(adm) == "EXTRAVASCULAR") {
if (!FLAGsteadystate) {
# For single dose data:
Res__["VZFO"] = dose/Res__["AUCIFO"]/Res__["LAMZ"]
Res__["VZFP"] = dose/Res__["AUCIFP"]/Res__["LAMZ"]
Res__["CLFO"] = dose/Res__["AUCIFO"]
Res__["CLFP"] = dose/Res__["AUCIFP"]
} else {
# For steady-state data only SS values are determined
Res__["CLFSS"] = dose/Res__["AUCTAU"]
Res__["VZFSS"] = dose/Res__["AUCTAU"]/Res__["LAMZ"]
}
Res__["MRTEVLST"] = Res__["AUMCLST"]/Res__["AUCLST"]
Res__["MRTEVIFO"] = Res__["AUMCIFO"]/Res__["AUCIFO"]
Res__["MRTEVIFP"] = Res__["AUMCIFP"]/Res__["AUCIFP"]
} else {
if (!FLAGsteadystate) {
# For single dose data:
Res__["VZO"] = dose/Res__["AUCIFO"]/Res__["LAMZ"]
Res__["VZP"] = dose/Res__["AUCIFP"]/Res__["LAMZ"]
Res__["CLO"] = dose/Res__["AUCIFO"]
Res__["CLP"] = dose/Res__["AUCIFP"]
Res__["VSSO"] = Res__["MRTIVIFO"]*Res__["CLO"]
Res__["VSSP"] = Res__["MRTIVIFP"]*Res__["CLP"]
} else {
# For steady-state data only SS values are determined
Res__["CLSS"] = dose/Res__["AUCTAU"]
Res__["VZSS"] = dose/Res__["AUCTAU"]/Res__["LAMZ"]
}
Res__["MRTIVLST"] = Res__["AUMCLST"]/Res__["AUCLST"] - dur/2
Res__["MRTIVIFO"] = Res__["AUMCIFO"]/Res__["AUCIFO"] - dur/2
Res__["MRTIVIFP"] = Res__["AUMCIFP"]/Res__["AUCIFP"] - dur/2
}
# Remove NA elments in Res__
Res__ <- Res__[!is.na(Res__)]
# Get units and potential conversion factors
Units = Unit_IQRnca(doseUnit=doseUnit, timeUnit=timeUnit, concUnit=concUnit, MW=MW)
# Apply unit conversion factors
for (i in 1:length(Res__)) Res__[i] = Res__[i] * Units[names(Res__[i]),2]
# Determine interval AUC
if (is.data.frame(iAUC)) { # eg iAUC = data.frame(Name=c("AUC[0-12h]","AUC[0-24h]"), Start=c(0,0), End=c(12,24))
niAUC__ = nrow(iAUC)
if (niAUC__ > 0) {
RetNames = union(RetNames1__, as.character(iAUC[,"Name"]))
for (i in 1:niAUC__) {
if (adm == "BOLUS") Res__[as.character(iAUC[i,"Name"])] = IntAUC_IQRnca(x2, y2, iAUC[i,"Start"], iAUC[i,"End"], Res__, down=down)
else Res__[as.character(iAUC[i,"Name"])] = IntAUC_IQRnca(x, y, iAUC[i,"Start"], iAUC[i,"End"], Res__, down=down)
Units <- rbind(Units, structure(Units[rownames(Units) == "AUCLST", ], row.names = as.character(iAUC[i,"Name"])))
}
}
} else {
niAUC__ = 0
}
if (is.null(IQRversion)) {
# To ensure compatibility with IQR Tools versions until 1.0.8
units__ <- sapply(names(Res__), function (n) Units[n,1])
} else {
# From IQR 1.0.9 its different
# Return units as named vector - contains all ... matching later
units__ <- Units$Unit
names(units__) <- rownames(Units)
}
attr(Res__, "units") <- units__
attr(Res__, "slopePoints") <- IXpointsSlope__
return(Res__)
}
########################################################################### ###
# getSlope_IQRnca function
########################################################################### ###
getSlope_IQRnca <- function(x, y)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.20
# Called by: BestSlope
# Calls: none except base
# INPUT
# x: time
# y: natural log of concentration
# RETURNS
Result__ = c(R2 = NA_real_, # R square
R2ADJ = NA_real_, # R square adjusted
LAMZNPT = 0, # Number of points for Lambda z
LAMZ = NA_real_, # Lambda z, terminal slope as a positive number
b0 = NA_real_, # intercept from OLS, i.e. simple linear regression
CORRXY = NA_real_, # Correlation of x, y
LAMZLL = NA_real_, # Lower time for lambda z
LAMZUL = NA_real_, # Upper time for lambda z
CLSTP = NA_real_) # Concentration last predicted in original scale
# Input Check
n = length(x)
if (n == 1 | n != length(y) | !is.numeric(x) | !is.numeric(y)) {
return(Result__) # return default
}
#
mx = mean(x)
my = mean(y)
Sxx = sum((x - mx)*(x - mx))
Sxy = sum((x - mx)*(y - my))
Syy = sum((y - my)*(y - my))
b1 = Sxy/Sxx
if (is.nan(b1) | b1 >= 0) {
return(Result__) # return default
}
# Expectedly
Result__["LAMZNPT"] = n
Result__["LAMZ"] = -b1
Result__["b0"] = my - b1*mx
Result__["R2"] = b1 * Sxy/Syy
Result__["R2ADJ"] = 1 - (1 - Result__["R2"])*(n - 1)/(n - 2)
Result__["CORRXY"] = sign(b1)*sqrt(Result__["R2"])
Result__["LAMZLL"] = x[1]
Result__["LAMZUL"] = x[n]
Result__["CLSTP"] = exp(Result__["b0"] + b1 * x[n])
return(Result__)
}
########################################################################### ###
# LogAUC_IQRnca function
########################################################################### ###
LogAUC_IQRnca <- function(x, y) # down="Log" means Linear-Up Log-Down
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: none except base
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# RETURNS
Result__ = c(AUC = NA_real_, # AUC by log down method
AUMC = NA_real_) # AUMC by log down method
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
auc = 0
aumc = 0
for (i in 2:n) {
if (y[i] < y[i-1] & y[i] > 0) {
k = (log(y[i-1]) - log(y[i]))/(x[i] - x[i-1]) # -k slope in y-log scale
auc = auc + (y[i-1] - y[i])/k
aumc = aumc + (x[i-1]*y[i-1] - x[i]*y[i])/k + (y[i-1] - y[i])/k/k
} else {
auc = auc + (x[i] - x[i-1])*(y[i] + y[i-1])/2
aumc = aumc + (x[i] - x[i-1])*(y[i]*x[i] + y[i-1]*x[i-1])/2
}
}
Result__["AUC"] = auc
Result__["AUMC"] = aumc
return(Result__)
}
########################################################################### ###
# LinAUC_IQRnca function
########################################################################### ###
LinAUC_IQRnca <- function(x, y) # down="Linear"
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: none except base
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# RETURNS
Result__ = c(AUC = NA_real_, # AUC by linear down method
AUMC = NA_real_) # AUMC by linear down method
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
Result__["AUC"] = sum((x[-1] - x[-n])*(y[-1] + y[-n]))/2
Result__["AUMC"] = sum((x[-1] - x[-n])*(x[-1]*y[-1] + x[-n]*y[-n]))/2
return(Result__)
}
###############################################################################
# Interpol_IQRnca function
###############################################################################
Interpol_IQRnca <- function(x, y, xnew__, Slope=0, b0=0, down="Linear")
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: UT
# INPUT
# x: time or similar vector
# y: concentration or similear vector
# xnew__: new x point to be interpolated
# Slope: terminal slope of log(y) ~ x
# b0: intercept of log(y) ~ x
# down: "Linear" or "Log"
# RETURNS
Result__ = list(x, y) # Default return value
# Input check
n = length(x)
if (n != length(y)) {
warning("Interpol: Length of x and y are different!")
newN = min(n, length(y))
x = x[1:newN]
y = y[1:newN]
}
if (!is.numeric(x) | !is.numeric(y) | !is.character(down)) {
return(Result__)
}
if (xnew__ %in% x) return(Result__)
if (sum(x < xnew__) > 0) {
LEFT__ = TRUE
x1 = x[max(which(x < xnew__))]
y1 = y[max(which(x < xnew__))]
} else LEFT__ = FALSE
if (sum(x > xnew__) > 0) {
RIGHT__ = TRUE
x2 = x[min(which(x > xnew__))]
y2 = y[min(which(x > xnew__))]
} else RIGHT__ = FALSE
if (LEFT__==TRUE & RIGHT__==TRUE) {
if (UT_IQRnca(down)=="LOG" & y2 < y1 & y2 > 0) {
ynew = exp(log(y1) + (log(y2) - log(y1))/(x2 - x1)*(xnew__ - x1))
} else {
ynew = y1 + (y2 - y1)/(x2 - x1)*(xnew__ - x1)
}
}
if (LEFT__==TRUE & RIGHT__==FALSE) ynew = exp(b0 - Slope*xnew__) # NOT ynew = exp(log(y1) - Slope*(xnew__ - x1))
if (LEFT__==FALSE & RIGHT__==TRUE) ynew = y2/x2 * xnew__
if (LEFT__==FALSE & RIGHT__==FALSE) return(Result__)
Result__ = list(sort(c(x, xnew__)), c(y, ynew)[order(c(x, xnew__))])
return(Result__)
}
###############################################################################
# IntAUC_IQRnca function
###############################################################################
IntAUC_IQRnca <- function(x, y, t1, t2, Res__, down="Linear")
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: singleNCA_IQRnca
# Calls: Interpol, LinAUC, LogAUC, UT
# INPUT
# x: time
# y: natural log of concentration
# t1: time of start. This could be interpolated, if this is not included in x.
# t2: time of end. This could be interpolated, if this is not included in x.
# Res__: result of singleNCA_IQRnca
# down: "Linear" or "Log"
# RETURN
Result__ = NA_real_ # Interval AUC
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
if (t1 > Res__["TLST"]) return(Result__)
tL__ = Res__["TLST"]
if (t2 > tL__ & is.na(Res__["LAMZ"])) return(Result__)
newSeries__ = Interpol_IQRnca(x, y, t1, Res__["LAMZ"], Res__["b0"], down=down)
newSeries__ = Interpol_IQRnca(newSeries__[[1]], newSeries__[[2]], t2, Res__["LAMZ"], Res__["b0"], down=down)
x = newSeries__[[1]]
y = newSeries__[[2]]
if (UT_IQRnca(down)=="LINEAR") {
if (t2 <= tL__) {
Result__ = LinAUC_IQRnca(x[x>=t1 & x<=t2], y[x>=t1 & x<=t2])[[1]]
} else {
Result__ = LinAUC_IQRnca(x[x>=t1 & x<=tL__], y[x>=t1 & x<=tL__])[[1]] + LogAUC_IQRnca(x[x>=tL__ & x<=t2], y[x>=tL__ & x<=t2])[[1]]
}
} else if (UT_IQRnca(down)=="LOG") {
Result__ = LogAUC_IQRnca(x[x>=t1 & x<=t2], y[x>=t1 & x<=t2])[[1]]
} else {
Result__ = NA_real_
}
return(Result__)
}
###############################################################################
# BestSlope_IQRnca function
###############################################################################
BestSlope_IQRnca <- function(x, y, adm="Extravascular", TOL=1e-4,FLAGivSlopeCmax=FALSE)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.19
# Called by : singleNCA_IQRnca
# Calls : Slope, UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# adm: method of drug administration "Bolus", "Infusion", or "Extravascular"
# TOL: Tolerance, see Phoneix WinNonlin 6.4 User's Guide p33
# RETURNS
Result__ = c(R2 = NA, # R square
R2ADJ = NA, # R square adjusted
LAMZNPT = 0, # Number of points for Lambda z
LAMZ = NA, # Lambda z, terminal slope as a positive number
b0 = NA, # intercept from OLS, i.e. simple linear regression
CORRXY = NA, # Correlation of x, y
LAMZLL = NA, # Lower time for lambda z
LAMZUL = NA, # Upper time for lambda z
CLSTP = NA) # Concentration last predicted in original scale
# Input Check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y) | length(y[y < 0]) > 0) {
Result__["LAMZNPT"] = 0
return(Result__)
}
if (length(unique(y)) == 1) { # Case of all the same values
Result__["LAMZNPT"] = 0
Result__["b0"] = unique(y)
return(Result__)
}
if (UT_IQRnca(adm) == "BOLUS" | (FLAGivSlopeCmax & UT_IQRnca(adm) == "INFUSION")) {
locStart__ = which.max(y) # From Tmax (for Bolus and (Infusion if FLAGivSlopeCmax=TRUE))
} else {
locStart__ = which.max(y) + 1 # From next after Tmax (for the others)
}
locLast__ = max(which(y > 0)) # Till non-zero concentration
if (locLast__ - locStart__ < 2) { # Too few to fit, if this is 2, R2ADJ becomes NaN.
Result__["LAMZNPT"] = 0
return(Result__)
}
tmpMat__ = matrix(nrow=(locLast__ - locStart__ - 1), ncol=length(Result__))
colnames(tmpMat__) = names(Result__)
for (i in locStart__:(locLast__ - 2)) {
tmpMat__[i - locStart__ + 1,] = getSlope_IQRnca(x[i:locLast__], log(y[i:locLast__]))
}
# Keep only slopes for which more than 2 points were used
tmpMat__ = tmpMat__[tmpMat__[,"LAMZNPT"] > 2,,drop=FALSE]
# Exclusion of slopes for which R2ADJ could not be determined
if (any(is.na(tmpMat__[,"R2ADJ"]))) {
warning("Adjusted Rsquared could not be determined for at least one of possible set of slope points. These sets cannot be considered.")
tmpMat__ <- tmpMat__[,!is.na(tmpMat__[,"R2ADJ"])]
}
if (nrow(tmpMat__) > 0) {
maxAdjRsq = max(tmpMat__[,"R2ADJ"])
OKs__ = ifelse(abs(maxAdjRsq - tmpMat__[,"R2ADJ"]) < TOL, TRUE, FALSE)
nMax__ = max(tmpMat__[OKs__,"LAMZNPT"])
Result__ = tmpMat__[OKs__ & tmpMat__[,"LAMZNPT"]==nMax__,]
} else {
Result__["LAMZNPT"] = 0
}
return(Result__)
}
###############################################################################
# AUC_IQRnca function
###############################################################################
AUC_IQRnca <- function(x, y, down)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by: singleNCA_IQRnca
# Calls: UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# down: LINEAR or LOG
# The switching from linear to log in case of down="LOG"
# will not be done on an individual sample pair basis but based on
# pre Tmax and post Tmax samples - as this is how WinNonlin is doing it.
n = length(x)
# RETURNS
Result__ = c(AUC = rep(NA_real_, n), # vector of cumulative AUC
AUMC = rep(NA_real_, n)) # vector of cumulative AUMC
# Input check
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
Res__ = matrix(nrow=n, ncol=2) # temporary result
Res__[1,] = c(0, 0)
for (i in 2:n) {
if (x[i-1] < min(x[which(y==max(y))])) { # New approach ... based on Tmax
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else if (UT_IQRnca(down) == "LINEAR") { # downward & linear
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else if (UT_IQRnca(down) == "LOG") { # downward & log
k = (log(y[i - 1]) - log(y[i]))/(x[i] - x[i-1]) # -k slope in y-log scale
if (k==0) {
# Special case of using linear as same value
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else {
# Use log
Res__[i,1] = (y[i-1] - y[i])/k
Res__[i,2] = (x[i-1]*y[i-1] - x[i]*y[i])/k + (y[i-1] - y[i])/k/k
}
} else {
Res__[i,1] = NA_real_
Res__[i,2] = NA_real_
}
}
Result__ = cbind(AUC=cumsum(Res__[,1]), AUMC=cumsum(Res__[,2]))
return(Result__)
}
###############################################################################
# UT_IQRnca function
###############################################################################
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
UT_IQRnca <- function(x) toupper(gsub("^\\s+|\\s+$", "", x))
###############################################################################
# UnitUrine_IQRnca function
###############################################################################
UnitUrine_IQRnca <- function(conU="ng/mL", volU="mL", amtU__="mg", MW=0)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by:
# Calls:
# INPUT
# conU: concentration unit
# volU: volume unit
# amtU__: amount unit, upper and lower case should match
# durU: duration unit, "h" fixed
# rateU: rate unit, amtU__/durU
# MW: molecular weight
# RETURNS conversion factor for the mulplication: conc * vol * factor = amount in the given unit
Result__ = 0 # default value to return
# Input check
if (!is.numeric(MW)) return (Result__)
if (MW < 0) return (Result__)
#
rGram__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rGram__) = c("g", "mg", "ug", "ng", "pg")
rVol__ = c(1, 1e1, 1e3, 1e6)
names(rVol__) = c("L", "dL", "mL", "uL")
rMol__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rMol__) = c("mol", "mmol", "umol", "nmol", "pmol")
amtU__ = tolower(amtU__)
if (toupper(conU) == toupper("mg/mL")) conU = "g/L"
if (toupper(conU) == toupper("ug/mL")) conU = "mg/L"
if (toupper(conU) == toupper("ng/mL")) conU = "ug/L"
if (toupper(conU) == toupper("pg/mL")) conU = "ng/L"
if (toupper(conU) == toupper("mmol/mL")) conU = "mol/L"
if (toupper(conU) == toupper("umul/mL")) conU = "mmol/L"
if (toupper(conU) == toupper("nmol/mL")) conU = "umol/L"
if (toupper(conU) == toupper("pmol/mL")) conU = "nmol/L"
amtU1__ = strsplit(conU, "/")[[1]][1]
if ((amtU1__ %in% names(rMol__) & amtU__ %in% names(rGram__)) | (amtU1__ %in% names(rGram__) & amtU__ %in% names(rMol__))) {
if (MW <= 0) {
warning("Molecular weight should be given!")
return (Result__)
}
}
if (amtU__ %in% names(rMol__) & amtU1__ %in% names(rGram__)) Result__ = rMol__[amtU__]/rGram__[amtU1__] / rVol__[volU] / MW
else if (amtU__ %in% names(rGram__) & amtU1__ %in% names(rMol__)) Result__ = rGram__[amtU__]/rMol__[amtU1__] / rVol__[volU] * MW
else if (amtU__ %in% names(rGram__) & amtU1__ %in% names(rGram__)) Result__ = rGram__[amtU__]/rGram__[amtU1__] / rVol__[volU]
else if (amtU__ %in% names(rMol__) & amtU1__ %in% names(rMol__)) Result__ = rMol__[amtU__]/rMol__[amtU1__] / rVol__[volU]
else {Result__ = 0}
return (Result__)
}
## CDISC TERMS PPTESTCD
# RCAMLST (recovery amount last)
# RCPCLST (recovery percent last) FracExcr Fraction excreted unchanged
# ERINT (excretion rate interval)
# LSTURINE: Last urine collection time from PCENDTC
# RENALCL (renal clearance) = RCAMLST/AUC[0-LSTURINE]
###############################################################################
# Unit_IQRnca function
###############################################################################
Unit_IQRnca <- function(code__="", timeUnit="h", concUnit="ng/mL", doseUnit="mg", MW=0)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by:
# Calls:
# INPUT
# code__: SDTM PPTESTCD
# timeUnit: time unit
# concUnit: concentration unit
# doseUnit: dose unit, this should not be amount per kg (BWT) or per m2 (BSA)
# MW: molecular weight
# RETURNS
Result__ = c(Unit = NA_character_, # unit of SDTM PPTESTCD like AUCLST, CMAX, CMAXD, ...
Factor = NA_real_) # conversion factor used internally
# Input check
if (length(strsplit(doseUnit, "/")[[1]]) != 1) return(Result__)
if (!is.numeric(MW)) return(Result__)
if (MW < 0) return(Result__)
#
rGram__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rGram__) = c("g", "mg", "ug", "ng", "pg")
rMol__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rMol__) = c("mol", "mmol", "umol", "nmol", "pmol")
doseUnit = tolower(doseUnit)
timeUnit = tolower(timeUnit)
if (toupper(concUnit) == toupper("mg/mL")) concUnit = "g/L"
if (toupper(concUnit) == toupper("ug/mL")) concUnit = "mg/L"
if (toupper(concUnit) == toupper("ng/mL")) concUnit = "ug/L"
if (toupper(concUnit) == toupper("pg/mL")) concUnit = "ng/L"
if (toupper(concUnit) == toupper("mmol/mL")) concUnit = "mol/L"
if (toupper(concUnit) == toupper("umol/mL")) concUnit = "mmol/L"
if (toupper(concUnit) == toupper("nmol/mL")) concUnit = "umol/L"
if (toupper(concUnit) == toupper("pmol/mL")) concUnit = "nmol/L"
tConc__ = strsplit(concUnit, "/")[[1]]
uAmt__ = tConc__[1]
uVol__ = tConc__[2]
if ((uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) | (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__))) {
if (is.na(MW)) warning("Molecular weight should be given for more informative results!")
if (MW <= 0) warning("Molecular weight should be given for more informative results!")
}
TestCD__ = c("b0", "CMAX", "CMAXD", "TMAX", "TLAG", "CLST", "CLSTP", "TLST", "LAMZHL", "LAMZ",
"LAMZLL", "LAMZUL", "LAMZNPT", "CORRXY", "R2", "R2ADJ", "C0", "AUCLST", "AUCALL",
"AUCIFO", "AUCIFOD", "AUCIFP", "AUCIFPD", "AUCPEO", "AUCPEP", "AUCPBEO", "AUCPBEP",
"AUMCLST", "AUMCIFO", "AUMCIFP", "AUMCPEO", "AUMCPEP",
"MRTIVLST", "MRTIVIFO", "MRTIVIFP", "MRTEVLST", "MRTEVIFO", "MRTEVIFP",
"VZO", "VZP", "VZFO", "VZFP", "CLO", "CLP", "CLFO", "CLFP", "VSSO", "VSSP","VZFSS","VZSS","CLFSS","CLSS","AUCTAU","AUCTAUD","TAU")
nTestCD = length(TestCD__)
Res__ = data.frame(Unit=rep("",nTestCD), Factor=rep(1,nTestCD), stringsAsFactors = FALSE)
rownames(Res__) = TestCD__
for (i in 1:nTestCD) {
Code__ = TestCD__[i]
if (Code__ %in% c("CMAX", "CLST", "CLSTP", "C0")) Res__[i, 1] = concUnit
if (Code__ == "CMAXD") Res__[i, 1] = paste0(concUnit,"/",doseUnit)
if (Code__ %in% c("TMAX", "TLAG", "TLST", "LAMZHL", "LAMZLL", "LAMZUL", "MRTIVLST",
"MRTIVIFO", "MRTIVIFP", "MRTEVLST", "MRTEVIFO", "MRTEVIFP","TAU")) Res__[i, 1] = timeUnit
if (Code__ == "LAMZ") Res__[i, 1] = paste0("1/",timeUnit)
if (Code__ %in% c("b0", "LAMZNPT", "CORRXY", "R2", "R2ADJ")) Res__[i, 1] = ""
if (Code__ %in% c("AUCLST", "AUCALL", "AUCIFO", "AUCIFP","AUCTAU")) Res__[i, 1] = paste0(timeUnit," * ",concUnit)
if (Code__ %in% c("AUCIFOD", "AUCIFPD","AUCTAUD")) Res__[i, 1] = paste0(timeUnit," * ",concUnit,"/",doseUnit)
if (Code__ %in% c("AUCPEO", "AUCPEP", "AUCPBEO", "AUCPBEP", "AUMCPEO", "AUMCPEP")) Res__[i, 1] = "%"
if (Code__ %in% c("AUMCLST", "AUMCIFO", "AUMCIFP")) Res__[i, 1] = paste0(timeUnit,"2 * ",concUnit)
if (Code__ %in% c("VZO", "VZP", "VZFO", "VZFP", "VSSO", "VSSP","VZFSS","VZSS")) {
if (uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) Res__[i, ] = c(uVol__, rMol__[uAmt__]/rGram__[doseUnit] / MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__)) Res__[i, ] = c(uVol__, rGram__[uAmt__]/rMol__[doseUnit] * MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rGram__)) Res__[i, ] = Res__[i, ] = c(uVol__, rGram__[uAmt__]/rGram__[doseUnit])
else Res__[i, ] = c(uVol__, rMol__[uAmt__]/rMol__[doseUnit])
}
if (Code__ %in% c("CLO", "CLP", "CLFO", "CLFP","CLFSS","CLSS")) {
if (uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rMol__[uAmt__]/rGram__[doseUnit] / MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__)) Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rGram__[uAmt__]/rMol__[doseUnit] * MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rGram__)) Res__[i, ] = Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rGram__[uAmt__]/rGram__[doseUnit])
else Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rMol__[uAmt__]/rMol__[doseUnit])
}
}
Res__[,2] = as.numeric(Res__[,2])
Res__[Res__[,2] == 0 | Res__[,2] == Inf, 2] = NA
if (code__ == "") Result__ = Res__ # return all codes
else return(Result__ = Res__[code__,]) # return only specific codes
return(Result__)
}
IntiQuan/NonCompartIQR documentation built on Nov. 8, 2019, 2:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions singleNCA_IQRnca getSlope_IQRnca LogAUC_IQRnca LinAUC_IQRnca Interpol_IQRnca IntAUC_IQRnca BestSlope_IQRnca AUC_IQRnca UT_IQRnca UnitUrine_IQRnca Unit_IQRnca
# Original Author: Kyun-Seop Bae k@acr.kr
# Adapted by: Henning Schmidt henning.schmidt@intiquan.com
# Based on NonCompart package on CRAN
########################################################################### ###
# singleNCA_IQRnca function
########################################################################### ###
singleNCA_IQRnca <- function(
x,
y,
dose=0,
adm="Extravascular",
dur=0,
doseUnit="mg",
timeUnit="h",
concUnit="ug/L",
NCAinfo=0, # Dummy for licensing check
iAUC="",
down="Log",
MW=0,
slopePoints = NULL,
NOSLOPE = NULL,
TOL = 1e-4,
FLAGsteadystate = FALSE,
tau = NULL,
IQRversion = NULL,
FLAGivSlopeCmax = FALSE) {
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.8.14
# Called by: tblNCA
# Calls: BestSlope, IntAUC, Unit, UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# dose: dose, actually given dose, not per kg or per square meter dose!
# adm: administration method, "Extravascular", "Bolus", or "Infusion"
# dur: duration of infusion
# doseUnit: dose unit (not per kg or per m2)
# timeUnit: time unit
# concUnit: concentration unit
# iAUC: data frame for interval AUC. See the example for the detail
# down: trapezoidal downward caculation. "Linear" or "Log"
# MW: molecular weight
# slopePoints: vector with indices of observations to use for slope calculation
# NOSLOPE: logical vector TRUE=>do not consider point for bestslope, FALSE: do consider
# IQRversion: version number for handling compatibility outside
# FLAGivSlopeCmax: if set to TRUE then Cmax obs will be allowed to be considered
# for bestslope calc. useful for 1 cpt dynamics and availability of very few points for NCA
# cannot be implemented by default as then not matching with other examples with richer data and
# comparison to WinNonLin
# RETURNS
# NCA result
# For steady-state data: Require presence of pre-dose sample at time=0 if infusion or extravascular administration
if (FLAGsteadystate) {
if (x[1] != 0) {
stop("Analysis of steady-state data requires availability of pre-dose sample\nin the data for infusion and extravascular administration")
}
if (is.na(y[1])) {
stop("Analysis of steady-state data requires availability of pre-dose sample (non-NA)\nin the data for infusion and extravascular administration")
}
# Check if last observation time point larger or equal to tau - otherwise error
if (tau > max(x)) {
stop("For steady-state profile analysis the last sample time needs to be larger than tau")
}
}
# Other checks
if (!(is.numeric(x) & is.numeric(y) & is.numeric(dose) & is.numeric(dur) & is.character(adm) & is.character(down))) {
stop("Check input types!")
}
if (UT_IQRnca(adm) == "INFUSION" & !all(dur > 0)) stop("Infusion mode should have dur larger than 0!")
###################################### ###
# Special Handling for slope calculation
###################################### ###
# Save original x and y values for determination of the right slope points after
# slope calculation (done via x1 to xorig__ matching). In addition xorig__ and yorig__
# are also used when manually points are selected for slope calculation
xorig__ <- x
yorig__ <- y
# Get xslope__, yslope__ to be used for bestslope algorithm
xslope__ <- x[!NOSLOPE] # Remove the points to be avoided in bestslope algorithm
yslope__ <- y[!NOSLOPE] # Remove the points to be avoided in bestslope algorithm
xslope__ = xslope__[yslope__ != 0] # Remove all points with zeros in y (including mid) for LAMZ
yslope__ = yslope__[yslope__ != 0] # Remove all points with zeros in y
NApoints__ <- is.na(xslope__) | is.na(yslope__)
xslope__ = xslope__[!NApoints__] # remove NA points in x
yslope__ = yslope__[!NApoints__] # remove NA points in y
###################################### ###
# Standard handling
###################################### ###
# Remove NA points
NApoints__ <- is.na(x) | is.na(y)
x <- x[!NApoints__] # remove NA points in x
y <- y[!NApoints__] # remove NA points in y
n <- length(x)
RetNames1__ = c("b0", "CMAX", "CMAXD", "TMAX", "TLAG", "CLST", "CLSTP", "TLST", "LAMZHL", "LAMZ",
"LAMZLL", "LAMZUL", "LAMZNPT", "CORRXY", "R2", "R2ADJ", "AUCLST", "AUCALL",
"AUCIFO", "AUCIFOD", "AUCIFP", "AUCIFPD", "AUCPEO", "AUCPEP",
"AUMCLST", "AUMCIFO", "AUMCIFP", "AUMCPEO", "AUMCPEP")
# If steady-state profile then add AUCTAU
if (FLAGsteadystate) {
RetNames1__ = c(RetNames1__, "AUCTAU","AUCTAUD","TAU")
}
if (UT_IQRnca(adm) == "BOLUS") {
RetNames1__ = union(RetNames1__, c("C0", "AUCPBEO", "AUCPBEP"))
}
if (UT_IQRnca(adm) == "EXTRAVASCULAR") {
# Some of them are dependent of SS and SD ... removed later what is not needed
RetNames1__ = union(RetNames1__, c("VZFO", "VZFP", "CLFO", "CLFP", "CLFSS", "VZFSS", "MRTEVLST", "MRTEVIFO", "MRTEVIFP"))
} else {
# Some of them are dependent of SS and SD ... removed later what is not needed
RetNames1__ = union(RetNames1__, c("VZO", "VZP", "CLO", "CLP", "CLSS", "VZSS", "MRTIVLST", "MRTIVIFO", "MRTIVIFP", "VSSO", "VSSP"))
}
Res__ = rep(NA_real_, length(RetNames1__))
names(Res__) = RetNames1__
if (n == 0 | n != length(y) | length(y[y < 0]) > 0) {
Res__["LAMZNPT"] = 0
return(Res__)
}
uY__ = unique(y)
if (length(uY__) == 1) { # Case of all the same values
Res__["CMAX"] = uY__
if (dose > 0) Res__["CMAXD"] = uY__ / dose
Res__["TMAX"] = x[y==uY__][1] # First Tmax
if (which(y==uY__)[1] > 1) {
Res__["TLAG"] = x[which(y==uY__) - 1]
} else {
Res__["TLAG"] = 0
}
Res__["CLST"] = uY__
Res__["TLST"] = x[y==uY__][1]
Res__["LAMZNPT"] = 0
Res__["b0"] = uY__
return(Res__)
}
iLastNonZero__ = max(which(y > 0)) # Index of last non-zero y
x0 = x[1:iLastNonZero__] # Till Non-zero concentration. i.e. removing trailing zeros
y0 = y[1:iLastNonZero__] # Till Non-zero concentration. i.e. removing trailing zeros
if (UT_IQRnca(adm) == "BOLUS") {
# Bolus - add an initial 0 at time 0 if not included in vector
# Works for steady-state and single dose data
if (is.na(x[x==0][1])) {
# Determine concentration at time 0
if (y[1] > y[2] & y[2] > 0) {
# Use extrapolation
C0 = exp(-x[1]*(log(y[2]) - log(y[1]))/(x[2] - x[1]) + log(y[1]))
} else {
# Use first non-zero value
C0 = y[x==min(x[y > 0])]
}
x2 = c(0, x)
y2 = c(C0, y)
x3 = c(0, x0)
y3 = c(C0, y0)
} else {
# 0 in x as first element. Check that y[1] > 0 otherwise error
if (y[1] <= 0) stop("For bolus administration the concentration at time 0 is not allowed to be 0")
x2 = x
y2 = y
x3 = x0
y3 = y0
# Set initial concentration at time 0 to observed one
C0 = y[1]
}
} else {
# Infusion & extravascular - add an initial 0 at time 0 if not included in vector
# Works for steady-state and single dose data
if (is.na(x[x==0][1])) {
# This case cannot happen for steady-state data as the 0 time point is required to be present and
# populated from the data.
x2 = c(0, x) # for AUCALL including trailing zero y values
y2 = c(0, y)
x3 = c(0, x0) # for AUCLST without trailing zero y values
y3 = c(0, y0)
} else {
# Time 0 included in x ... no need to add
# Works for single dose and steady-state dose!
x2 = x # for AUCALL including trailing zero y values
y2 = y
x3 = x0 # for AUCLST without trailing zero y values
y3 = y0
}
}
if (is.null(slopePoints)) {
tRes__ <- BestSlope_IQRnca(xslope__, yslope__, adm, TOL=TOL,FLAGivSlopeCmax)
# Return indices of points used for slope calculation
Indices_pointsSlope__ <- rev(seq(length(xslope__), by=-1, length.out = tRes__["LAMZNPT"]))
IXpointsSlope__ <- (1:length(x))[xorig__ %in% xslope__][Indices_pointsSlope__]
} else {
# Call the getSlope_IQRnca function
tRes__ <- getSlope_IQRnca(xorig__[slopePoints], log(yorig__[slopePoints]))
# Return indices of points used for slope calculation
IXpointsSlope__ <- slopePoints
}
if (length(tRes__) != 9) tRes__ = c(NA, NA, 0, NA, NA, NA, NA, NA, NA)
Res__[c("R2", "R2ADJ", "LAMZNPT", "LAMZ", "b0", "CORRXY", "LAMZLL", "LAMZUL", "CLSTP")] = tRes__
tabAUC = AUC_IQRnca(x3, y3, down)
Res__[c("AUCLST","AUMCLST")] = tabAUC[length(x3),]
Res__["AUCALL"] = AUC_IQRnca(x2, y2, down)[length(x2),1]
Res__["LAMZHL"] = log(2)/Res__["LAMZ"]
Res__["TMAX"] = x[which.max(y)][1]
Res__["CMAX"] = max(y)
Res__["TLST"] = x[iLastNonZero__]
Res__["CLST"] = y[iLastNonZero__]
Res__["AUCIFO"] = Res__["AUCLST"] + Res__["CLST"]/Res__["LAMZ"]
Res__["AUCIFP"] = Res__["AUCLST"] + Res__["CLSTP"]/Res__["LAMZ"]
Res__["AUCPEO"] = (1 - Res__["AUCLST"]/Res__["AUCIFO"])*100
Res__["AUCPEP"] = (1 - Res__["AUCLST"]/Res__["AUCIFP"])*100
Res__["AUMCIFO"] = Res__["AUMCLST"] + Res__["CLST"]*Res__["TLST"]/Res__["LAMZ"] + Res__["CLST"]/Res__["LAMZ"]/Res__["LAMZ"]
Res__["AUMCIFP"] = Res__["AUMCLST"] + Res__["CLSTP"]*Res__["TLST"]/Res__["LAMZ"] + Res__["CLSTP"]/Res__["LAMZ"]/Res__["LAMZ"]
Res__["AUMCPEO"] = (1 - Res__["AUMCLST"]/Res__["AUMCIFO"])*100
Res__["AUMCPEP"] = (1 - Res__["AUMCLST"]/Res__["AUMCIFP"])*100
if (!is.na(dose) & dose > 0) {
Res__["CMAXD"] = Res__["CMAX"]/dose
Res__["AUCIFOD"] = Res__["AUCIFO"]/dose
Res__["AUCIFPD"] = Res__["AUCIFP"]/dose
}
if (UT_IQRnca(adm) == "BOLUS") {
Res__["C0"] = C0 # Phoneix WinNonlin 6.4 User's Guide p27
Res__["AUCPBEO"] = tabAUC[2,1]/Res__["AUCIFO"]*100
Res__["AUCPBEP"] = tabAUC[2,1]/Res__["AUCIFP"]*100
} else {
if (sum(y0==0) > 0) Res__["TLAG"] = x0[max(which(y0==0))] # Trailing zero should not exist
else Res__["TLAG"] = 0
if (!is.na(x0[x0==0][1])) {
if (y0[x0==0] > 0) Res__["TLAG"] = 0 # This is WinNonlin logic
}
}
# Calculate AUCTAU using IntAUC_IQRnca for SS data
if (FLAGsteadystate) {
Res__["AUCTAU"] = IntAUC_IQRnca(x, y, 0, tau, Res__, down=down)
Res__["AUCTAUD"] = Res__["AUCTAU"]/dose
Res__["TAU"] = tau
}
if (UT_IQRnca(adm) == "EXTRAVASCULAR") {
if (!FLAGsteadystate) {
# For single dose data:
Res__["VZFO"] = dose/Res__["AUCIFO"]/Res__["LAMZ"]
Res__["VZFP"] = dose/Res__["AUCIFP"]/Res__["LAMZ"]
Res__["CLFO"] = dose/Res__["AUCIFO"]
Res__["CLFP"] = dose/Res__["AUCIFP"]
} else {
# For steady-state data only SS values are determined
Res__["CLFSS"] = dose/Res__["AUCTAU"]
Res__["VZFSS"] = dose/Res__["AUCTAU"]/Res__["LAMZ"]
}
Res__["MRTEVLST"] = Res__["AUMCLST"]/Res__["AUCLST"]
Res__["MRTEVIFO"] = Res__["AUMCIFO"]/Res__["AUCIFO"]
Res__["MRTEVIFP"] = Res__["AUMCIFP"]/Res__["AUCIFP"]
} else {
if (!FLAGsteadystate) {
# For single dose data:
Res__["VZO"] = dose/Res__["AUCIFO"]/Res__["LAMZ"]
Res__["VZP"] = dose/Res__["AUCIFP"]/Res__["LAMZ"]
Res__["CLO"] = dose/Res__["AUCIFO"]
Res__["CLP"] = dose/Res__["AUCIFP"]
Res__["VSSO"] = Res__["MRTIVIFO"]*Res__["CLO"]
Res__["VSSP"] = Res__["MRTIVIFP"]*Res__["CLP"]
} else {
# For steady-state data only SS values are determined
Res__["CLSS"] = dose/Res__["AUCTAU"]
Res__["VZSS"] = dose/Res__["AUCTAU"]/Res__["LAMZ"]
}
Res__["MRTIVLST"] = Res__["AUMCLST"]/Res__["AUCLST"] - dur/2
Res__["MRTIVIFO"] = Res__["AUMCIFO"]/Res__["AUCIFO"] - dur/2
Res__["MRTIVIFP"] = Res__["AUMCIFP"]/Res__["AUCIFP"] - dur/2
}
# Remove NA elments in Res__
Res__ <- Res__[!is.na(Res__)]
# Get units and potential conversion factors
Units = Unit_IQRnca(doseUnit=doseUnit, timeUnit=timeUnit, concUnit=concUnit, MW=MW)
# Apply unit conversion factors
for (i in 1:length(Res__)) Res__[i] = Res__[i] * Units[names(Res__[i]),2]
# Determine interval AUC
if (is.data.frame(iAUC)) { # eg iAUC = data.frame(Name=c("AUC[0-12h]","AUC[0-24h]"), Start=c(0,0), End=c(12,24))
niAUC__ = nrow(iAUC)
if (niAUC__ > 0) {
RetNames = union(RetNames1__, as.character(iAUC[,"Name"]))
for (i in 1:niAUC__) {
if (adm == "BOLUS") Res__[as.character(iAUC[i,"Name"])] = IntAUC_IQRnca(x2, y2, iAUC[i,"Start"], iAUC[i,"End"], Res__, down=down)
else Res__[as.character(iAUC[i,"Name"])] = IntAUC_IQRnca(x, y, iAUC[i,"Start"], iAUC[i,"End"], Res__, down=down)
Units <- rbind(Units, structure(Units[rownames(Units) == "AUCLST", ], row.names = as.character(iAUC[i,"Name"])))
}
}
} else {
niAUC__ = 0
}
if (is.null(IQRversion)) {
# To ensure compatibility with IQR Tools versions until 1.0.8
units__ <- sapply(names(Res__), function (n) Units[n,1])
} else {
# From IQR 1.0.9 its different
# Return units as named vector - contains all ... matching later
units__ <- Units$Unit
names(units__) <- rownames(Units)
}
attr(Res__, "units") <- units__
attr(Res__, "slopePoints") <- IXpointsSlope__
return(Res__)
}
########################################################################### ###
# getSlope_IQRnca function
########################################################################### ###
getSlope_IQRnca <- function(x, y)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.20
# Called by: BestSlope
# Calls: none except base
# INPUT
# x: time
# y: natural log of concentration
# RETURNS
Result__ = c(R2 = NA_real_, # R square
R2ADJ = NA_real_, # R square adjusted
LAMZNPT = 0, # Number of points for Lambda z
LAMZ = NA_real_, # Lambda z, terminal slope as a positive number
b0 = NA_real_, # intercept from OLS, i.e. simple linear regression
CORRXY = NA_real_, # Correlation of x, y
LAMZLL = NA_real_, # Lower time for lambda z
LAMZUL = NA_real_, # Upper time for lambda z
CLSTP = NA_real_) # Concentration last predicted in original scale
# Input Check
n = length(x)
if (n == 1 | n != length(y) | !is.numeric(x) | !is.numeric(y)) {
return(Result__) # return default
}
#
mx = mean(x)
my = mean(y)
Sxx = sum((x - mx)*(x - mx))
Sxy = sum((x - mx)*(y - my))
Syy = sum((y - my)*(y - my))
b1 = Sxy/Sxx
if (is.nan(b1) | b1 >= 0) {
return(Result__) # return default
}
# Expectedly
Result__["LAMZNPT"] = n
Result__["LAMZ"] = -b1
Result__["b0"] = my - b1*mx
Result__["R2"] = b1 * Sxy/Syy
Result__["R2ADJ"] = 1 - (1 - Result__["R2"])*(n - 1)/(n - 2)
Result__["CORRXY"] = sign(b1)*sqrt(Result__["R2"])
Result__["LAMZLL"] = x[1]
Result__["LAMZUL"] = x[n]
Result__["CLSTP"] = exp(Result__["b0"] + b1 * x[n])
return(Result__)
}
########################################################################### ###
# LogAUC_IQRnca function
########################################################################### ###
LogAUC_IQRnca <- function(x, y) # down="Log" means Linear-Up Log-Down
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: none except base
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# RETURNS
Result__ = c(AUC = NA_real_, # AUC by log down method
AUMC = NA_real_) # AUMC by log down method
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
auc = 0
aumc = 0
for (i in 2:n) {
if (y[i] < y[i-1] & y[i] > 0) {
k = (log(y[i-1]) - log(y[i]))/(x[i] - x[i-1]) # -k slope in y-log scale
auc = auc + (y[i-1] - y[i])/k
aumc = aumc + (x[i-1]*y[i-1] - x[i]*y[i])/k + (y[i-1] - y[i])/k/k
} else {
auc = auc + (x[i] - x[i-1])*(y[i] + y[i-1])/2
aumc = aumc + (x[i] - x[i-1])*(y[i]*x[i] + y[i-1]*x[i-1])/2
}
}
Result__["AUC"] = auc
Result__["AUMC"] = aumc
return(Result__)
}
########################################################################### ###
# LinAUC_IQRnca function
########################################################################### ###
LinAUC_IQRnca <- function(x, y) # down="Linear"
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: none except base
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# RETURNS
Result__ = c(AUC = NA_real_, # AUC by linear down method
AUMC = NA_real_) # AUMC by linear down method
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
Result__["AUC"] = sum((x[-1] - x[-n])*(y[-1] + y[-n]))/2
Result__["AUMC"] = sum((x[-1] - x[-n])*(x[-1]*y[-1] + x[-n]*y[-n]))/2
return(Result__)
}
###############################################################################
# Interpol_IQRnca function
###############################################################################
Interpol_IQRnca <- function(x, y, xnew__, Slope=0, b0=0, down="Linear")
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: IntAUC
# Calls: UT
# INPUT
# x: time or similar vector
# y: concentration or similear vector
# xnew__: new x point to be interpolated
# Slope: terminal slope of log(y) ~ x
# b0: intercept of log(y) ~ x
# down: "Linear" or "Log"
# RETURNS
Result__ = list(x, y) # Default return value
# Input check
n = length(x)
if (n != length(y)) {
warning("Interpol: Length of x and y are different!")
newN = min(n, length(y))
x = x[1:newN]
y = y[1:newN]
}
if (!is.numeric(x) | !is.numeric(y) | !is.character(down)) {
return(Result__)
}
if (xnew__ %in% x) return(Result__)
if (sum(x < xnew__) > 0) {
LEFT__ = TRUE
x1 = x[max(which(x < xnew__))]
y1 = y[max(which(x < xnew__))]
} else LEFT__ = FALSE
if (sum(x > xnew__) > 0) {
RIGHT__ = TRUE
x2 = x[min(which(x > xnew__))]
y2 = y[min(which(x > xnew__))]
} else RIGHT__ = FALSE
if (LEFT__==TRUE & RIGHT__==TRUE) {
if (UT_IQRnca(down)=="LOG" & y2 < y1 & y2 > 0) {
ynew = exp(log(y1) + (log(y2) - log(y1))/(x2 - x1)*(xnew__ - x1))
} else {
ynew = y1 + (y2 - y1)/(x2 - x1)*(xnew__ - x1)
}
}
if (LEFT__==TRUE & RIGHT__==FALSE) ynew = exp(b0 - Slope*xnew__) # NOT ynew = exp(log(y1) - Slope*(xnew__ - x1))
if (LEFT__==FALSE & RIGHT__==TRUE) ynew = y2/x2 * xnew__
if (LEFT__==FALSE & RIGHT__==FALSE) return(Result__)
Result__ = list(sort(c(x, xnew__)), c(y, ynew)[order(c(x, xnew__))])
return(Result__)
}
###############################################################################
# IntAUC_IQRnca function
###############################################################################
IntAUC_IQRnca <- function(x, y, t1, t2, Res__, down="Linear")
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.24
# Called by: singleNCA_IQRnca
# Calls: Interpol, LinAUC, LogAUC, UT
# INPUT
# x: time
# y: natural log of concentration
# t1: time of start. This could be interpolated, if this is not included in x.
# t2: time of end. This could be interpolated, if this is not included in x.
# Res__: result of singleNCA_IQRnca
# down: "Linear" or "Log"
# RETURN
Result__ = NA_real_ # Interval AUC
# Input check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
if (t1 > Res__["TLST"]) return(Result__)
tL__ = Res__["TLST"]
if (t2 > tL__ & is.na(Res__["LAMZ"])) return(Result__)
newSeries__ = Interpol_IQRnca(x, y, t1, Res__["LAMZ"], Res__["b0"], down=down)
newSeries__ = Interpol_IQRnca(newSeries__[[1]], newSeries__[[2]], t2, Res__["LAMZ"], Res__["b0"], down=down)
x = newSeries__[[1]]
y = newSeries__[[2]]
if (UT_IQRnca(down)=="LINEAR") {
if (t2 <= tL__) {
Result__ = LinAUC_IQRnca(x[x>=t1 & x<=t2], y[x>=t1 & x<=t2])[[1]]
} else {
Result__ = LinAUC_IQRnca(x[x>=t1 & x<=tL__], y[x>=t1 & x<=tL__])[[1]] + LogAUC_IQRnca(x[x>=tL__ & x<=t2], y[x>=tL__ & x<=t2])[[1]]
}
} else if (UT_IQRnca(down)=="LOG") {
Result__ = LogAUC_IQRnca(x[x>=t1 & x<=t2], y[x>=t1 & x<=t2])[[1]]
} else {
Result__ = NA_real_
}
return(Result__)
}
###############################################################################
# BestSlope_IQRnca function
###############################################################################
BestSlope_IQRnca <- function(x, y, adm="Extravascular", TOL=1e-4,FLAGivSlopeCmax=FALSE)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification: 2017.7.19
# Called by : singleNCA_IQRnca
# Calls : Slope, UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# adm: method of drug administration "Bolus", "Infusion", or "Extravascular"
# TOL: Tolerance, see Phoneix WinNonlin 6.4 User's Guide p33
# RETURNS
Result__ = c(R2 = NA, # R square
R2ADJ = NA, # R square adjusted
LAMZNPT = 0, # Number of points for Lambda z
LAMZ = NA, # Lambda z, terminal slope as a positive number
b0 = NA, # intercept from OLS, i.e. simple linear regression
CORRXY = NA, # Correlation of x, y
LAMZLL = NA, # Lower time for lambda z
LAMZUL = NA, # Upper time for lambda z
CLSTP = NA) # Concentration last predicted in original scale
# Input Check
n = length(x)
if (n != length(y) | !is.numeric(x) | !is.numeric(y) | length(y[y < 0]) > 0) {
Result__["LAMZNPT"] = 0
return(Result__)
}
if (length(unique(y)) == 1) { # Case of all the same values
Result__["LAMZNPT"] = 0
Result__["b0"] = unique(y)
return(Result__)
}
if (UT_IQRnca(adm) == "BOLUS" | (FLAGivSlopeCmax & UT_IQRnca(adm) == "INFUSION")) {
locStart__ = which.max(y) # From Tmax (for Bolus and (Infusion if FLAGivSlopeCmax=TRUE))
} else {
locStart__ = which.max(y) + 1 # From next after Tmax (for the others)
}
locLast__ = max(which(y > 0)) # Till non-zero concentration
if (locLast__ - locStart__ < 2) { # Too few to fit, if this is 2, R2ADJ becomes NaN.
Result__["LAMZNPT"] = 0
return(Result__)
}
tmpMat__ = matrix(nrow=(locLast__ - locStart__ - 1), ncol=length(Result__))
colnames(tmpMat__) = names(Result__)
for (i in locStart__:(locLast__ - 2)) {
tmpMat__[i - locStart__ + 1,] = getSlope_IQRnca(x[i:locLast__], log(y[i:locLast__]))
}
# Keep only slopes for which more than 2 points were used
tmpMat__ = tmpMat__[tmpMat__[,"LAMZNPT"] > 2,,drop=FALSE]
# Exclusion of slopes for which R2ADJ could not be determined
if (any(is.na(tmpMat__[,"R2ADJ"]))) {
warning("Adjusted Rsquared could not be determined for at least one of possible set of slope points. These sets cannot be considered.")
tmpMat__ <- tmpMat__[,!is.na(tmpMat__[,"R2ADJ"])]
}
if (nrow(tmpMat__) > 0) {
maxAdjRsq = max(tmpMat__[,"R2ADJ"])
OKs__ = ifelse(abs(maxAdjRsq - tmpMat__[,"R2ADJ"]) < TOL, TRUE, FALSE)
nMax__ = max(tmpMat__[OKs__,"LAMZNPT"])
Result__ = tmpMat__[OKs__ & tmpMat__[,"LAMZNPT"]==nMax__,]
} else {
Result__["LAMZNPT"] = 0
}
return(Result__)
}
###############################################################################
# AUC_IQRnca function
###############################################################################
AUC_IQRnca <- function(x, y, down)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by: singleNCA_IQRnca
# Calls: UT
# INPUT
# x: time or similar vector
# y: concentration or similar vector
# down: LINEAR or LOG
# The switching from linear to log in case of down="LOG"
# will not be done on an individual sample pair basis but based on
# pre Tmax and post Tmax samples - as this is how WinNonlin is doing it.
n = length(x)
# RETURNS
Result__ = c(AUC = rep(NA_real_, n), # vector of cumulative AUC
AUMC = rep(NA_real_, n)) # vector of cumulative AUMC
# Input check
if (n != length(y) | !is.numeric(x) | !is.numeric(y)) return(Result__)
Res__ = matrix(nrow=n, ncol=2) # temporary result
Res__[1,] = c(0, 0)
for (i in 2:n) {
if (x[i-1] < min(x[which(y==max(y))])) { # New approach ... based on Tmax
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else if (UT_IQRnca(down) == "LINEAR") { # downward & linear
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else if (UT_IQRnca(down) == "LOG") { # downward & log
k = (log(y[i - 1]) - log(y[i]))/(x[i] - x[i-1]) # -k slope in y-log scale
if (k==0) {
# Special case of using linear as same value
Res__[i,1] = (x[i] - x[i-1])*(y[i] + y[i-1])/2
Res__[i,2] = (x[i] - x[i-1])*(x[i]*y[i] + x[i-1]*y[i-1])/2
} else {
# Use log
Res__[i,1] = (y[i-1] - y[i])/k
Res__[i,2] = (x[i-1]*y[i-1] - x[i]*y[i])/k + (y[i-1] - y[i])/k/k
}
} else {
Res__[i,1] = NA_real_
Res__[i,2] = NA_real_
}
}
Result__ = cbind(AUC=cumsum(Res__[,1]), AUMC=cumsum(Res__[,2]))
return(Result__)
}
###############################################################################
# UT_IQRnca function
###############################################################################
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
UT_IQRnca <- function(x) toupper(gsub("^\\s+|\\s+$", "", x))
###############################################################################
# UnitUrine_IQRnca function
###############################################################################
UnitUrine_IQRnca <- function(conU="ng/mL", volU="mL", amtU__="mg", MW=0)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by:
# Calls:
# INPUT
# conU: concentration unit
# volU: volume unit
# amtU__: amount unit, upper and lower case should match
# durU: duration unit, "h" fixed
# rateU: rate unit, amtU__/durU
# MW: molecular weight
# RETURNS conversion factor for the mulplication: conc * vol * factor = amount in the given unit
Result__ = 0 # default value to return
# Input check
if (!is.numeric(MW)) return (Result__)
if (MW < 0) return (Result__)
#
rGram__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rGram__) = c("g", "mg", "ug", "ng", "pg")
rVol__ = c(1, 1e1, 1e3, 1e6)
names(rVol__) = c("L", "dL", "mL", "uL")
rMol__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rMol__) = c("mol", "mmol", "umol", "nmol", "pmol")
amtU__ = tolower(amtU__)
if (toupper(conU) == toupper("mg/mL")) conU = "g/L"
if (toupper(conU) == toupper("ug/mL")) conU = "mg/L"
if (toupper(conU) == toupper("ng/mL")) conU = "ug/L"
if (toupper(conU) == toupper("pg/mL")) conU = "ng/L"
if (toupper(conU) == toupper("mmol/mL")) conU = "mol/L"
if (toupper(conU) == toupper("umul/mL")) conU = "mmol/L"
if (toupper(conU) == toupper("nmol/mL")) conU = "umol/L"
if (toupper(conU) == toupper("pmol/mL")) conU = "nmol/L"
amtU1__ = strsplit(conU, "/")[[1]][1]
if ((amtU1__ %in% names(rMol__) & amtU__ %in% names(rGram__)) | (amtU1__ %in% names(rGram__) & amtU__ %in% names(rMol__))) {
if (MW <= 0) {
warning("Molecular weight should be given!")
return (Result__)
}
}
if (amtU__ %in% names(rMol__) & amtU1__ %in% names(rGram__)) Result__ = rMol__[amtU__]/rGram__[amtU1__] / rVol__[volU] / MW
else if (amtU__ %in% names(rGram__) & amtU1__ %in% names(rMol__)) Result__ = rGram__[amtU__]/rMol__[amtU1__] / rVol__[volU] * MW
else if (amtU__ %in% names(rGram__) & amtU1__ %in% names(rGram__)) Result__ = rGram__[amtU__]/rGram__[amtU1__] / rVol__[volU]
else if (amtU__ %in% names(rMol__) & amtU1__ %in% names(rMol__)) Result__ = rMol__[amtU__]/rMol__[amtU1__] / rVol__[volU]
else {Result__ = 0}
return (Result__)
}
## CDISC TERMS PPTESTCD
# RCAMLST (recovery amount last)
# RCPCLST (recovery percent last) FracExcr Fraction excreted unchanged
# ERINT (excretion rate interval)
# LSTURINE: Last urine collection time from PCENDTC
# RENALCL (renal clearance) = RCAMLST/AUC[0-LSTURINE]
###############################################################################
# Unit_IQRnca function
###############################################################################
Unit_IQRnca <- function(code__="", timeUnit="h", concUnit="ng/mL", doseUnit="mg", MW=0)
{
# Author: Kyun-Seop Bae k@acr.kr
# Last modification; 2017-07-24
# Called by:
# Calls:
# INPUT
# code__: SDTM PPTESTCD
# timeUnit: time unit
# concUnit: concentration unit
# doseUnit: dose unit, this should not be amount per kg (BWT) or per m2 (BSA)
# MW: molecular weight
# RETURNS
Result__ = c(Unit = NA_character_, # unit of SDTM PPTESTCD like AUCLST, CMAX, CMAXD, ...
Factor = NA_real_) # conversion factor used internally
# Input check
if (length(strsplit(doseUnit, "/")[[1]]) != 1) return(Result__)
if (!is.numeric(MW)) return(Result__)
if (MW < 0) return(Result__)
#
rGram__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rGram__) = c("g", "mg", "ug", "ng", "pg")
rMol__ = c(1, 1e3, 1e6, 1e9, 1e12)
names(rMol__) = c("mol", "mmol", "umol", "nmol", "pmol")
doseUnit = tolower(doseUnit)
timeUnit = tolower(timeUnit)
if (toupper(concUnit) == toupper("mg/mL")) concUnit = "g/L"
if (toupper(concUnit) == toupper("ug/mL")) concUnit = "mg/L"
if (toupper(concUnit) == toupper("ng/mL")) concUnit = "ug/L"
if (toupper(concUnit) == toupper("pg/mL")) concUnit = "ng/L"
if (toupper(concUnit) == toupper("mmol/mL")) concUnit = "mol/L"
if (toupper(concUnit) == toupper("umol/mL")) concUnit = "mmol/L"
if (toupper(concUnit) == toupper("nmol/mL")) concUnit = "umol/L"
if (toupper(concUnit) == toupper("pmol/mL")) concUnit = "nmol/L"
tConc__ = strsplit(concUnit, "/")[[1]]
uAmt__ = tConc__[1]
uVol__ = tConc__[2]
if ((uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) | (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__))) {
if (is.na(MW)) warning("Molecular weight should be given for more informative results!")
if (MW <= 0) warning("Molecular weight should be given for more informative results!")
}
TestCD__ = c("b0", "CMAX", "CMAXD", "TMAX", "TLAG", "CLST", "CLSTP", "TLST", "LAMZHL", "LAMZ",
"LAMZLL", "LAMZUL", "LAMZNPT", "CORRXY", "R2", "R2ADJ", "C0", "AUCLST", "AUCALL",
"AUCIFO", "AUCIFOD", "AUCIFP", "AUCIFPD", "AUCPEO", "AUCPEP", "AUCPBEO", "AUCPBEP",
"AUMCLST", "AUMCIFO", "AUMCIFP", "AUMCPEO", "AUMCPEP",
"MRTIVLST", "MRTIVIFO", "MRTIVIFP", "MRTEVLST", "MRTEVIFO", "MRTEVIFP",
"VZO", "VZP", "VZFO", "VZFP", "CLO", "CLP", "CLFO", "CLFP", "VSSO", "VSSP","VZFSS","VZSS","CLFSS","CLSS","AUCTAU","AUCTAUD","TAU")
nTestCD = length(TestCD__)
Res__ = data.frame(Unit=rep("",nTestCD), Factor=rep(1,nTestCD), stringsAsFactors = FALSE)
rownames(Res__) = TestCD__
for (i in 1:nTestCD) {
Code__ = TestCD__[i]
if (Code__ %in% c("CMAX", "CLST", "CLSTP", "C0")) Res__[i, 1] = concUnit
if (Code__ == "CMAXD") Res__[i, 1] = paste0(concUnit,"/",doseUnit)
if (Code__ %in% c("TMAX", "TLAG", "TLST", "LAMZHL", "LAMZLL", "LAMZUL", "MRTIVLST",
"MRTIVIFO", "MRTIVIFP", "MRTEVLST", "MRTEVIFO", "MRTEVIFP","TAU")) Res__[i, 1] = timeUnit
if (Code__ == "LAMZ") Res__[i, 1] = paste0("1/",timeUnit)
if (Code__ %in% c("b0", "LAMZNPT", "CORRXY", "R2", "R2ADJ")) Res__[i, 1] = ""
if (Code__ %in% c("AUCLST", "AUCALL", "AUCIFO", "AUCIFP","AUCTAU")) Res__[i, 1] = paste0(timeUnit," * ",concUnit)
if (Code__ %in% c("AUCIFOD", "AUCIFPD","AUCTAUD")) Res__[i, 1] = paste0(timeUnit," * ",concUnit,"/",doseUnit)
if (Code__ %in% c("AUCPEO", "AUCPEP", "AUCPBEO", "AUCPBEP", "AUMCPEO", "AUMCPEP")) Res__[i, 1] = "%"
if (Code__ %in% c("AUMCLST", "AUMCIFO", "AUMCIFP")) Res__[i, 1] = paste0(timeUnit,"2 * ",concUnit)
if (Code__ %in% c("VZO", "VZP", "VZFO", "VZFP", "VSSO", "VSSP","VZFSS","VZSS")) {
if (uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) Res__[i, ] = c(uVol__, rMol__[uAmt__]/rGram__[doseUnit] / MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__)) Res__[i, ] = c(uVol__, rGram__[uAmt__]/rMol__[doseUnit] * MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rGram__)) Res__[i, ] = Res__[i, ] = c(uVol__, rGram__[uAmt__]/rGram__[doseUnit])
else Res__[i, ] = c(uVol__, rMol__[uAmt__]/rMol__[doseUnit])
}
if (Code__ %in% c("CLO", "CLP", "CLFO", "CLFP","CLFSS","CLSS")) {
if (uAmt__ %in% names(rMol__) & doseUnit %in% names(rGram__)) Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rMol__[uAmt__]/rGram__[doseUnit] / MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rMol__)) Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rGram__[uAmt__]/rMol__[doseUnit] * MW)
else if (uAmt__ %in% names(rGram__) & doseUnit %in% names(rGram__)) Res__[i, ] = Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rGram__[uAmt__]/rGram__[doseUnit])
else Res__[i, ] = c(paste0(uVol__,"/",timeUnit), rMol__[uAmt__]/rMol__[doseUnit])
}
}
Res__[,2] = as.numeric(Res__[,2])
Res__[Res__[,2] == 0 | Res__[,2] == Inf, 2] = NA
if (code__ == "") Result__ = Res__ # return all codes
else return(Result__ = Res__[code__,]) # return only specific codes
return(Result__)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.