version 0.8.0

Author: Kyun-Seop Bae

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: wnl
-Version: 0.7.3
+Version: 0.8.0
 Title: Minimization Tool for Pharmacokinetic-Pharmacodynamic Data
         Analysis
 Description: This is a set of minimization tools (maximum likelihood estimation and least square fitting) to solve examples in the Johan Gabrielsson and Dan Weiner's book "Pharmacokinetic and Pharmacodynamic Data Analysis - Concepts and Applications" 5th ed. (ISBN:9198299107). Examples include linear and nonlinear compartmental model, turn-over model, single or multiple dosing bolus/infusion/oral models, allometry, toxicokinetics, reversible metabolism, in-vitro/in-vivo extrapolation, enterohepatic circulation, metabolite modeling, Emax model, inhibitory model, tolerance model, oscillating response model, enantiomer interaction model, effect compartment model, drug-drug interaction model, receptor occupancy model, and rebound phenomena model. 
@@ -12,5 +12,5 @@ NeedsCompilation: no
 LazyLoad: yes
 Repository: CRAN
 URL: https://cran.r-project.org/package=wnl
-Packaged: 2023-04-17 21:02:37 UTC; Kyun-SeopBae
-Date/Publication: 2023-04-25 03:30:03 UTC
+Packaged: 2024-02-22 05:56:36 UTC; Kyun-SeopBae
+Date/Publication: 2024-02-22 06:20:02 UTC

MD5

@@ -1,25 +1,25 @@
-dcd531c376e9100eee542e50ae21b588 *DESCRIPTION
+b4ec2498de8b0b66680530cb18261254 *DESCRIPTION
 390de0872c5e034f3bb65922b8b78051 *NAMESPACE
-83d6bb93ca7f2ebe336d3d1c70e6d5de *R/BasicUtil.R
+afc57c175ab7a762077e71ac3effeed1 *R/BasicUtil.R
 ba1255e411249f4db6f2caa66ee30cfe *R/Comp1.R
 6ee187dd007ec9a4d1193101fb331c53 *R/EnvObj.R
 cf5ef0c31a8433b8f966716f203471b9 *R/ExpandDH.R
-82c430639b2aec87a2d8257e1ed87b19 *R/Objs.R
-f7a61f78602501325418063cc00df4d4 *R/Secondary.R
+a3c994e66e093b3202e7ab69da35fc93 *R/Objs.R
+8f65358ebfbc3a2450251fc0f686d17b *R/Secondary.R
 5d9a05e422d4e7fd3413f288aaea93a8 *R/SolComp2.R
 7791bb80ef7ff7935d783d3550d08bb7 *R/SolComp3.R
 b801ef09ac2e1d9c931f1d08fe4a6b90 *R/cmpChi.R
 752433b1b66c620b33c3def811efdc8f *R/dx.R
 9bf4095904c0bf59982cb603c2b69e5c *R/hSkew.R
 ca032758dcb79150beac3ef4564b7cb8 *R/nComp.R
-6c3ab16618f11b24fa73adf2d1083d1e *R/nlr.R
+b087399546d4cedaf961d69ceb6238d8 *R/nlr.R
 969cb4bf5cab158d9f4b9fbd260ae4e1 *R/pComp.R
 29c541a8d30080132cd50a98484542db *R/pProf.R
 3a85a468ae0dc2bbb899524ea01feccf *R/sysdata.rda
 1c9e1b28728c0f1fc96017efc560b0d6 *R/wnl5.R
-f9e97cf825c853b1c7161b06fb90cc50 *inst/NEWS.Rd
+96768ba9460fc589c2a4a40084325289 *inst/NEWS.Rd
 75463c4d3bd42432719b597f93dd859d *inst/doc/Compartment-vignette.pdf
-13988a98e0e282530ebd7bdcf151ab0b *inst/doc/wnl-manual.pdf
+00ccc7e8b63db769da280561f0329a64 *inst/doc/wnl-manual.pdf
 57010db082832f52da4856e4b3521524 *man/BasicUtil.Rd
 f9ab397c3e89d4acba0f7c5f06e91923 *man/Comp1.Rd
 6c87b1f721849279be48c7474dee83b5 *man/DAT.Rd
@@ -34,7 +34,7 @@ e1e7543809d439770aa7d1cb6e897193 *man/dx.Rd
 ee0ea2d2db69f309039e2d70a791ac0a *man/e.Rd
 d60a1682e9d37ff616877f1ce4738665 *man/hSkew.Rd
 62aedc1f108c168426e5402b454b82e1 *man/nComp.Rd
-4990a914b5d7730be858fd9462f9e750 *man/nlr.Rd
+ae6d73065d75619d1a1a83086569cdf2 *man/nlr.Rd
 3a1b70982f3cf42a7bd0a912b81ca182 *man/pComp.Rd
 2c87a9fd98a811d39fd3431176c7b46e *man/pProf.Rd
 8640ba18696633870f706cd63f6d440c *man/wnl-package.Rd

R/BasicUtil.R

@@ -164,16 +164,14 @@ nHessian = function(fx, x)
   return(H)
 }
 
-g2inv = function(A, Augmented=FALSE, eps=1e-8)
+g2inv = function(A, eps=1e-8)
 {
   idx = abs(diag(A)) > eps
   p = sum(idx, na.rm=T)
-  p0 = ifelse(Augmented, p - 1, p)
-  if (p == 0 | p0 < 1) { A[, ] = 0 ; attr(A, "rank") = 0 ; return(A) }
+  if (p == 0) { M[, ] = 0 ; attr(M, "rank") = 0 ; return(M) }
   B = A[idx, idx, drop=F]
-
   r = 0
-  for (k in 1:p0) {
+  for (k in 1:p) {
     d = B[k, k]
     if (abs(d) < eps) { B[k, ] = 0 ; B[, k] = 0 ; next }
     B[k, ] = B[k, ]/d
@@ -188,10 +186,11 @@ g2inv = function(A, Augmented=FALSE, eps=1e-8)
     B[k, k] = 1/d
   }
 
-  A[!idx, !idx] = 0
-  A[idx, idx] = B
-  attr(A, "rank") = r
-  return(A)
+  M = matrix(0, nrow=NCOL(A), ncol=NROW(A))
+  M[1:r, 1:r] = B[1:r, 1:r]
+  attr(M, "rank") = r
+
+  return(M)
 }
 
 Hougaard = function(J, H, ssq)

R/Objs.R

@@ -1,33 +1,46 @@
+#ObjEst = function(vPara)
+#{
+#  b0 = exp(vPara - e$alpha)
+#  vPara = vector(length=e$nPara0)
+#  vPara[e$toEst] = b0/(b0 + 1)*(e$UB - e$LB) + e$LB
+#  if (e$fix[1] != 0) vPara[e$fix] = e$IE0[e$fix]
+#  return(e$Obj(vPara))
+#}
+
 ObjEst = function(vPara)
 {
   b0 = exp(vPara - e$alpha)
   vPara = b0/(b0 + 1)*(e$UB - e$LB) + e$LB
   return(e$Obj(vPara))
 }
 
-ObjDef = function(vPara) # Default Obj for Covariance Step (original obj)
+ObjDef = function(vPara0) # Default Obj for Covariance Step (original obj)
 {
-  Fi = e$Fx(vPara[1:e$nTheta])
+  vPara = vector(length=e$nPara0)
+  vPara[e$toEst] = vPara0
+  if (e$fix[1] != 0) vPara[e$fix] = e$IE0[e$fix]
+
+  Fi = e$Fx(vPara[1:e$nTheta0])
   Ri = e$Y - Fi
   e$Fi = Fi
 
   if (e$Error == "A") {
-    Ci = rep(vPara[e$SGindex], e$nRec)
+    Ci = rep(vPara[e$SGindex0], e$nRec)
   } else if (e$Error == "POIS") {
-    Ci = vPara[e$SGindex]*Fi     # Caution on Fi is zero
+    Ci = vPara[e$SGindex0]*Fi     # Caution on Fi is zero
     Ci[Fi == 0] = 1 # Weight cannot be calculated with zero values
     Ri[Fi == 0] = 0
   } else if (e$Error == "P") {
-    Ci = vPara[e$SGindex]*Fi*Fi  # Caution on Fi is zero
+    Ci = vPara[e$SGindex0]*Fi*Fi  # Caution on Fi is zero
     Ci[Fi == 0] = 1 # Weight cannot be calculated with zero values
     Ri[Fi == 0] = 0
   } else if (e$Error == "C") {
-    Ci = rep(vPara[e$SGindex[1]], e$nRec) + vPara[e$SGindex[2]]*Fi*Fi
+    Ci = rep(vPara[e$SGindex0[1]], e$nRec) + vPara[e$SGindex0[2]]*Fi*Fi
 #  } else if (e$Error == "CFA") { # Additive fixed propertional
 #    Ci = (e$FASD + sqrt(vPara[e$SGindex[1]])*Fi)^2
   } else if (e$Error == "S") {
-    Si = e$Sx(vPara[1:e$nTheta])
-    Ci = vPara[e$SGindex]*Si
+    Si = e$Sx(vPara[1:e$nTheta0])
+    Ci = vPara[e$SGindex0]*Si
     Ci[Si == 0] = 1 # Weight cannot be calculated with zero values
     Ri[Si == 0] = 0
   }

R/Secondary.R

@@ -4,9 +4,9 @@ Secondary = function(Formula, PE, COV)
   gr0 = deriv(Formula, names(PE), function.arg=names(PE), func=TRUE)
   gr1 = do.call("gr0", as.list(PE))
   PE2 = gr1[1]
-  gr2 = attr(gr1,"gradient") 
+  gr2 = attr(gr1,"gradient")
   SE2 = sqrt(gr2 %*% COV %*% t(gr2))
-  CV2 = SE2/PE2*100  
+  CV2 = SE2/PE2*100
   Result = c(PE2, SE2, CV2)
   names(Result) = c("PE", "SE", "RSE")
   return(Result)

R/nlr.R

@@ -1,4 +1,4 @@
-nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms, Method="L-BFGS-B", Sx, conf.level=0.95, k)
+nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms, Method="L-BFGS-B", Sx, conf.level=0.95, k, fix=0)
 {
 #  e = new.env(parent=globalenv()) # environment should exist before call this function
   t1 = Sys.time()
@@ -10,14 +10,13 @@ nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames,
     e$nRec = nrow(e$Y)
   } else {
     if (!("DV" %in% colnames(Data))) stop("Data should have 'DV' column.")
-    if (mean(abs(Data[, "DV"])) < 1e-6 | mean(abs(Data[,"DV"])) > 1e6) warning("DV is too large or too small. Rescale it.")
+    if (mean(abs(Data[, "DV"])) < 1e-6 | mean(abs(Data[, "DV"])) > 1e6) warning("DV is too large or too small. Rescale it.")
     e$Y  = Data[, "DV"] # Observation values, Data should have "DV" column.
     e$nRec = length(e$Y)
     if (sum(is.na(Data[, "DV"])) > 0) stop("DV column should not have NAs.")
   }
 
   if (length(pNames) != length(IE)) stop("pNames and IE should match.")
-
   e$pNames = pNames # parameter names in the order of Fx arguments
   e$IE = IE # initial estimate of Fx arguments
   e$nTheta = length(IE)
@@ -86,10 +85,37 @@ nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames,
 #   e$Ci = cbind(rep(SG1*SG1, e$nRec), rep(SG2*SG2, e$nRec))
 #   e$SumLogCi = sum(log(e$Ci))
   }
+
+  e$IE0 = e$IE
+  e$nPara0 = e$nPara
+  e$nTheta0 = e$nTheta
+  e$SGindex0 = e$SGindex
+  e$fix0 = fix
+  e$fix = e$fix0
+  e$toEst = which(!(1:e$nPara %in% e$fix))
+  e$pNames= e$pNames[e$toEst]
+  e$IE = e$IE[e$toEst]
+  e$LB = e$LB[e$toEst]
+  e$UB = e$UB[e$toEst]
+  e$nTheta = sum(1:e$nTheta %in% e$toEst)
+  e$nPara = length(e$toEst)
+  e$nEps = sum((1 + e$nTheta0):e$nPara0 %in% e$toEst)
+  if (e$nEps > 0) {
+    e$SGindex = (e$nTheta + 1):(e$nTheta + e$nEps) # index of error variance(s)
+  } else {
+    e$SGindex = 0
+  }
+  e$alpha = e$alpha[e$toEst]
+
   e$r = optim(rep(0.1, e$nPara), ObjEst, method=Method)
   e$PE = exp(e$r$par - e$alpha)/(exp(e$r$par - e$alpha) + 1)*(e$UB - e$LB) + e$LB
 
+#  e$PE0 = vector(length=e$nPara0)
+#  e$PE0[e$toEst] = e$PE
+#  e$PE0[e$fix] = e$IE0[e$fix]
+#  e$InvCov = hessian(e$Obj, e$PE0)/2  # FinalEst from EstStep()
   e$InvCov = hessian(e$Obj, e$PE)/2  # FinalEst from EstStep()
+  e$InvCov = e$InvCov
   e$Cov = try(solve(e$InvCov), silent=T)
   if (!is.matrix(e$Cov)) {
     e$Cov = g2inv(e$InvCov)
@@ -109,6 +135,19 @@ nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames,
   e$Correl = cov2cor(e$Cov)
   e$EigenVal = eigen(e$Correl)$values
 
+## Hougaard Skewness
+  if (e$Error == "A") {
+    e$PE0 = vector(length=e$nPara0)
+    e$PE0[e$toEst] = e$PE
+    e$PE0[e$fix] = e$IE0[e$fix]
+    
+    e$J = nGradient(e$Fx, e$PE0[1:e$nTheta0])
+    e$H = nHessian(e$Fx, e$PE0[1:e$nTheta0])
+    e$HouSkew = Hougaard(e$J, e$H, e$PE0[e$SGindex0])
+    names(e$HouSkew) = pNames
+  }
+
+## SE
   if (e$SGindex[1] > 0) {
     e$PE = c(e$PE, sqrt(e$PE[e$SGindex]))
     e$SE = c(e$SE, e$SE[e$SGindex]/2/sqrt(e$PE[e$SGindex])) # Delta method, See Wackerly p484, gr = (x^0.5)' = 0.5x^(-0.5), gr^2 = 1/(4*x)
@@ -158,13 +197,6 @@ nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames,
     e$nRec = sum(e$Y != -1)
   }
 
-## Hougaard Skewness
-  if (e$Error == "A") {
-    e$J = nGradient(e$Fx, e$PE[1:e$nTheta])
-    e$H = nHessian(e$Fx, e$PE[1:e$nTheta])
-    e$HouSkew = Hougaard(e$J, e$H, e$PE[e$SGindex])
-  }
-
 ## Likelihood Profile
   cAdd = e$nRec*log(2*pi)
   nRes = 51
@@ -255,14 +287,35 @@ nlr = function(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames,
 
   options(warn=0) # end of calling uniroot
 
+  if (e$SGindex[1] > 0) {
+    e$LI = cbind(e$LI, sqrt(e$LI[, e$SGindex]))
+    colnames(e$LI) = colnames(e$Est)[1:(e$nPara + length(e$SGindex))]
+  }
+
+  if (!missing(SecNames)) {
+    nSec = length(SecNames)
+    tRes = matrix(nrow=2, ncol=nSec)
+    colnames(tRes) = SecNames
+    rownames(tRes) = c("LL", "UL")
+    for (i in 1:nSec) {
+      fx2 = deriv(SecForms[[i]], e$pNames, function.arg=e$pNames, func=TRUE)
+      tv1 = do.call("fx2", as.list(e$LI[1, e$pNames]))
+      tv2 = do.call("fx2", as.list(e$LI[2, e$pNames]))
+      tRes[, i] = sort(c(tv1, tv2))
+    }
+    e$LI = cbind(e$LI, tRes)
+  }
+  colnames(e$LI) = colnames(e$Est)
+  attr(e$LI, "k") = exp(logk)
+
   e$mOFV = matrix(nrow=nRes, ncol=e$nPara)
   e$mPar = matrix(nrow=nRes, ncol=e$nPara)
   colnames(e$mOFV) = e$pNames[1:e$nPara]
   colnames(e$mPar) = e$pNames[1:e$nPara]
   for (j in 1:e$nPara) {
     e$mPar[, j] = seq(e$mParB[1, j], e$mParB[2, j], length.out=nRes)
     for (i in 1:nRes) {
-      tPar = e$PE
+      tPar = e$PE[1:e$nPara0]
       tPar[j] = e$mPar[i, j]
       e$mOFV[i, j] = cAdd + e$Obj(tPar) # -2LL
     }

inst/NEWS.Rd

@@ -1,10 +1,14 @@
 \name{NEWS}
 \title{News for Package \pkg{wnl}}
-
+\section{Version 0.8.0 (2024-02-22)}{
+  \itemize{
+    \item{g2inv utility function is revised to handle non-square matrix.}
+  }
+}
 
 \section{Version 0.7.3 (2023-04-25)}{
   \itemize{
-    \item{If default k value in nlr() is adjusted.}
+    \item{Default k value in nlr() is adjusted.}
   }
 }
 

inst/doc/wnl-manual.pdf

@@ -6,7 +6,7 @@ It performs nonlinear regression usually for pharmacokinetic and pharmacodynamic
 }
 \usage{
 nlr(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms, 
-    Method="L-BFGS-B", Sx, conf.level=0.95, k)
+    Method="L-BFGS-B", Sx, conf.level=0.95, k, fix=0)
 }
 \arguments{
   \item{Fx}{Function for structural model. It should return a vector of the same length to observations.}
@@ -23,6 +23,7 @@ nlr(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms,
   \item{Sx}{Scale function. This is usually the inverse of weight. It should return the same length(nrow) of Y. When Error="S", Scale function should be provided as \code{Sx}.}
   \item{conf.level}{Confidence level for confidence interval}
   \item{k}{1/k likelihood interval(LI) will be provided. Currently recommended value is exp(qf(1 - alpha, 1, nRec-nPara)/2) + 1.}
+  \item{fix}{indices of parameters to fix}
 }
 \details{
 This uses scaled transformed parameters and environment \code{e} internally.
@@ -75,5 +76,15 @@ This uses scaled transformed parameters and environment \code{e} internally.
     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)
 }
 \keyword{Maximum Likelihood Estimation}