billdenney / pknca

An R package is designed to perform all noncompartmental analysis (NCA) calculations for pharmacokinetic (PK) data.
http://billdenney.github.io/pknca/
GNU Affero General Public License v3.0
67 stars 24 forks source link

Feature Request: pk.calc.auctau #203

Closed asybing closed 7 months ago

asybing commented 1 year ago

Could a function be written that calculates phoenix's AUCtau (apologies if it has been written under a different function name)?

In most cases AUClast = AUCtau except when patients do not all have the same tlast. For example, for an interval of 24 hours, most patients will have concentrations at 24 hours but some will only have concentrations 16 hours (Not a BLQ at 24 hours, but the data is just not available yet). It would be helpful to add a function in PKNCA to calculate auctau as such:

pk.calc.auctau <- function(auclast,tau,clast.obs,clast.pred,tlast,lambda.z) {

  ctau.pred = clast.pred*exp(-lambda.z*(tau-tlast))

  # linear calculation
  auctaulin = auclast + 0.5*(tau-tlast)*(clast.obs+ctau.pred)

  # Log down calculation
  auctaulog = auclast + (tau-tlast)*(clast.obs-ctau.pred)/log(clast.obs/ctau.pred) # Phoenix default
}

Andrew

billdenney commented 1 year ago

Hi @asybing , I think that you're looking for the pk.calc.aucint function. Please let me know if that's what you're thinking of.

billdenney commented 1 year ago

If you're still having an issue, please reopen the issue.

asybing commented 1 year ago

Hi Bill, sorry for the late response I didn't get a notification.

aucint.all works in cases where end=tlast, but not in cases where end>tlast. I have the PKNCA auc.method="linear" to match phoenix's default, however since ctau needs to be extrapolated, phoenix uses the "log-linear" trapezoidal rule to calculate, whereas it seems PKNCA is still using a "linear" partial AUC calculation. Using the pk.calc.auctau calculation I've written above, I match Phoenix's AUCtau exactly when end>tlast.

Could you share the extrapolation partial AUC calculation aucint.all uses in cases where end>tlast?

billdenney commented 12 months ago

I believe that the parameter you're looking for is aucint.inf.obs. See the example below for a comparison of some of the aucint options:

library(PKNCA)
#> 
#> Attaching package: 'PKNCA'
#> The following object is masked from 'package:stats':
#> 
#>     filter

d_conc <-
  data.frame(
    conc = c(0, 1, 2, 0.75, 0.5, 0.2),
    time = c(0, 1, 2, 3, 4, 5)
  )

o_conc <- PKNCAconc(data = d_conc, conc~time)
d_intervals <-
  data.frame(
    start = 0, end = 6,
    aucint.inf.obs = TRUE,
    aucint.all = TRUE,
    aucint.last = TRUE
  )
o_data <- PKNCAdata(o_conc, intervals = d_intervals)
as.data.frame(
  pk.nca(o_data)
)
#> No dose information provided, calculations requiring dose will return NA.
#> # A tibble: 14 × 5
#>    start   end PPTESTCD            PPORRES exclude
#>    <dbl> <dbl> <chr>                 <dbl> <chr>  
#>  1     0     6 aucint.last           4.22  <NA>   
#>  2     0     6 aucint.all            4.32  <NA>   
#>  3     0     6 tmax                  2     <NA>   
#>  4     0     6 tlast                 5     <NA>   
#>  5     0     6 clast.obs             0.2   <NA>   
#>  6     0     6 lambda.z              0.661 <NA>   
#>  7     0     6 r.squared             0.953 <NA>   
#>  8     0     6 adj.r.squared         0.905 <NA>   
#>  9     0     6 lambda.z.time.first   3     <NA>   
#> 10     0     6 lambda.z.n.points     3     <NA>   
#> 11     0     6 clast.pred            0.218 <NA>   
#> 12     0     6 half.life             1.05  <NA>   
#> 13     0     6 span.ratio            1.91  <NA>   
#> 14     0     6 aucint.inf.obs        4.36  <NA>

Created on 2023-09-09 with reprex v2.0.2

asybing commented 12 months ago

Hi Bill, It looks like the aucint.inf.obs using a linear extrapolation. Is there a way to change it to a log extrapolation (while keeping PKNCA auc.method="linear")? This is what Phoenix AUCtau does with AUC calculation set to linear, extrapolation is done using log.

Andrew

billdenney commented 12 months ago

Let's chat about it.

I just took a look in the extrapolation code, and it gives a concentration of 0.1032796 at 6 hr which is extrapolation from the 5-hour concentration with a lambda.z of 0.6608779. I'm not sure how that would be giving a linear extrapolation (aufint.all does linear extrapolation).

asybing commented 12 months ago

the Cpredtau calculation Cpredlast*exp(-lambdaz*(tau-tpredlast)) is the same in Phoenix, but the difference is the linear vs log AUC calculation below

linear: 0.5*(tau-tlast)*(clast.obs+ctau.pred)

log: (tau-tlast)*(clast.obs-ctau.pred)/log(clast.obs/ctau.pred)

Andrew

billdenney commented 12 months ago

I just double-checked, and it is using the log-down calculation. What I did below is calculate through 5 hr (no extrapolation) and 6 hr (with extrapolation), then I calculated the values manually, and they line up with what you indicated so that the AUC between 5 and 6 hr is 0.1463515.

library(PKNCA)
#> 
#> Attaching package: 'PKNCA'
#> The following object is masked from 'package:stats':
#> 
#>     filter
d_conc <-
  data.frame(
    conc = c(0, 1, 2, 0.75, 0.5, 0.2),
    time = c(0, 1, 2, 3, 4, 5)
  )

o_conc <- PKNCAconc(data = d_conc, conc~time)
d_intervals <-
  data.frame(
    start = 0, end = c(5, 6),
    aucint.inf.obs = TRUE#,
    # aucint.all = TRUE,
    # aucint.last = TRUE
  )
o_data <- PKNCAdata(o_conc, intervals = d_intervals)
as.data.frame(as.data.frame(
  pk.nca(o_data)
))
#> No dose information provided, calculations requiring dose will return NA.
#>    start end            PPTESTCD   PPORRES exclude
#> 1      0   5                tmax 2.0000000    <NA>
#> 2      0   5               tlast 5.0000000    <NA>
#> 3      0   5           clast.obs 0.2000000    <NA>
#> 4      0   5            lambda.z 0.6608779    <NA>
#> 5      0   5           r.squared 0.9525736    <NA>
#> 6      0   5       adj.r.squared 0.9051472    <NA>
#> 7      0   5 lambda.z.time.first 3.0000000    <NA>
#> 8      0   5   lambda.z.n.points 3.0000000    <NA>
#> 9      0   5          clast.pred 0.2177734    <NA>
#> 10     0   5           half.life 1.0488279    <NA>
#> 11     0   5          span.ratio 1.9068906    <NA>
#> 12     0   5      aucint.inf.obs 4.2184147    <NA>
#> 13     0   6                tmax 2.0000000    <NA>
#> 14     0   6               tlast 5.0000000    <NA>
#> 15     0   6           clast.obs 0.2000000    <NA>
#> 16     0   6            lambda.z 0.6608779    <NA>
#> 17     0   6           r.squared 0.9525736    <NA>
#> 18     0   6       adj.r.squared 0.9051472    <NA>
#> 19     0   6 lambda.z.time.first 3.0000000    <NA>
#> 20     0   6   lambda.z.n.points 3.0000000    <NA>
#> 21     0   6          clast.pred 0.2177734    <NA>
#> 22     0   6           half.life 1.0488279    <NA>
#> 23     0   6          span.ratio 1.9068906    <NA>
#> 24     0   6      aucint.inf.obs 4.3647661    <NA>

#(tau-tlast)*(clast.obs-ctau.pred)/log(clast.obs/ctau.pred)
(6-5)*(0.2-0.1032796)/log(0.2/0.1032796)
#> [1] 0.1463515
4.2184147 + 0.1463515
#> [1] 4.364766

Created on 2023-09-09 with reprex v2.0.2

billdenney commented 12 months ago

@asybing, Could you try the simple dataset above to see if this is an issue? Specifically, could you try calculating AUCTAU from 0 to 5 hours and from 0 to 6 hours and indicate the results with at least 5 digits of precision after the decimal place? It would also be useful to have lambda.z from that to ensure that there is nothing different happening there, too.

From our discussion and my testing, I think that I'm already doing what you're asking, but as there appears to still be an issue, I'm not sure I can track it down with more granularity without some detailed specifics on what is happening with the values you expect.

asybing commented 12 months ago

@billdenney here is the code below. When I switch PKNCA.options(auc.method="linear"), the extrapolation also switches to linear. To match Phoenix the AUC calculation should be linear, but the extrapolation should be log.

R4.2.3|2020-09-01> library(PKNCA); PKNCA.options(auc.method="linear")
R4.2.3|2020-09-01> 
R4.2.3|2020-09-01> d_conc <-
+   data.frame(
+     ID = c(1,1,1,1,1,1),
+     conc = c(0, 1, 2, 0.75, 0.5, 0.2),
+     time = c(0, 1, 2, 3, 4, 5)
+   )
R4.2.3|2020-09-01> 
R4.2.3|2020-09-01> o_conc <- PKNCAconc(conc~time|ID,data = d_conc)
R4.2.3|2020-09-01> d_intervals <-
+   data.frame(
+     start = 0, end = c(5, 6),
+     aucint.inf.obs = TRUE
+   )
R4.2.3|2020-09-01> o_data <- pk.nca(PKNCAdata(o_conc, intervals = d_intervals))
No dose information provided, calculations requiring dose will return NA.
R4.2.3|2020-09-01> ncaresult <- as.data.frame(o_data);ncaresult
   start end ID            PPTESTCD   PPORRES exclude
1      0   5  1                tmax 2.0000000    <NA>
2      0   5  1               tlast 5.0000000    <NA>
3      0   5  1           clast.obs 0.2000000    <NA>
4      0   5  1            lambda.z 0.6608779    <NA>
5      0   5  1           r.squared 0.9525736    <NA>
6      0   5  1       adj.r.squared 0.9051472    <NA>
7      0   5  1 lambda.z.time.first 3.0000000    <NA>
8      0   5  1   lambda.z.n.points 3.0000000    <NA>
9      0   5  1          clast.pred 0.2177734    <NA>
10     0   5  1           half.life 1.0488279    <NA>
11     0   5  1          span.ratio 1.9068906    <NA>
12     0   5  1      aucint.inf.obs 4.3500000    <NA>
13     0   6  1                tmax 2.0000000    <NA>
14     0   6  1               tlast 5.0000000    <NA>
15     0   6  1           clast.obs 0.2000000    <NA>
16     0   6  1            lambda.z 0.6608779    <NA>
17     0   6  1           r.squared 0.9525736    <NA>
18     0   6  1       adj.r.squared 0.9051472    <NA>
19     0   6  1 lambda.z.time.first 3.0000000    <NA>
20     0   6  1   lambda.z.n.points 3.0000000    <NA>
21     0   6  1          clast.pred 0.2177734    <NA>
22     0   6  1           half.life 1.0488279    <NA>
23     0   6  1          span.ratio 1.9068906    <NA>
24     0   6  1      aucint.inf.obs 4.5016398    <NA>
R4.2.3|2020-09-01> 
R4.2.3|2020-09-01> 
R4.2.3|2020-09-01> aucinf.obs5 <- ncaresult$PPORRES[12]; aucinf.obs6 <- ncaresult$PPORRES[24]
R4.2.3|2020-09-01> aucinf.obs6-aucinf.obs5
[1] 0.1516398
R4.2.3|2020-09-01> 
R4.2.3|2020-09-01> lambdaz <- ncaresult$PPORRES[4]
R4.2.3|2020-09-01> ctau.pred <- .2*exp(-lambdaz*(6-5))
R4.2.3|2020-09-01> 0.5*(6-5)*(0.2+ctau.pred)
[1] 0.1516398
R4.2.3|2020-09-01> (6-5)*(0.2-ctau.pred)/log(0.2/ctau.pred)
[1] 0.1463515
billdenney commented 12 months ago

AHA! I think that I finally understand the issue. (I think that I was the dense one here.) Let me say it all back to you to see if I've really got it now:

When

  1. calculating an AUC with a partial extrapolation (aucint.inf.obs in PKNCA language, AUCTAU in WinNonlin language) and
  2. the end of the interval is after the last concentration above the LLOQ and
  3. the linear trapezoidal rule is being used for integration

the extrapolation after tlast should use log-down integration instead of linear trapezoidal integration.

Did I get all of that right?

If I got it right (finally), I think that I agree and the fix should be pretty easy.

asybing commented 12 months ago

yes that’s exactly right

glad to hear it’s an easy fix!

billdenney commented 12 months ago

For future reference, I'm going to simplify the overall issue to the following: If extrapolating using half-life, use log-down extrapolation. And, this should occur regardless of the AUC calculation method.

asybing commented 12 months ago

great, thank you bill! when do you think the change will be available?

asybing commented 7 months ago

thanks bill. happy new year!!Andrew--Andrew B. @. | (240) 505-5990On Jan 12, 2024, at 3:46 PM, Bill Denney @.> wrote: Closed #203 as completed via #263.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

asybing commented 3 months ago

Hi @billdenney, hope you’ve been well! Are you planning an update soon to incorporate this change?

billdenney commented 2 months ago

Version 0.11.0 is now on CRAN! Thanks for your patience.