gracehung
☆    

Taiwan,
2013-05-20 09:35
(3965 d 14:20 ago)

Posting: # 10603
Views: 19,756
 

 handling of missing data [NCA / SHAM]

I have one question about the handling of the missing data, and need your advice.
My question is:
As in BE study, One sample’s actual sampling time was miss recorded. However as according to protocol, this sample was analyzed by our analytical laboratory.
  1. Considering from PK personnel’s aspect, as it was not specified in protocol, how would you handling this missed time point at statistical analysis?
  2. Will this data be included in statistical analysis? If yes, what is the handling procedure?
  3. If the data is excluded, what kind of method would you use in calculate your AUC?
    1. Would you adapt interpolation method?
    2. If yes, would you adapt logarithmic interpolation or linear interpolation?
  4. Is there any guideline or reference document we could follow?
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-20 16:22
(3965 d 07:34 ago)

@ gracehung
Posting: # 10604
Views: 17,690
 

 Lin-up/log-down trapezoidal avoids most trouble

Hi Gracehung,

please follow the Forum’s Policy in the future. THX.

❝ As in BE study, One sample’s actual sampling time was miss recorded. However as according to protocol, this sample was analyzed by our analytical laboratory.


Perform an audit. Why was an existing sample recorded as missing?

❝ 1. it was not specified in protocol…


Don’t worry, it’s too late. :-D

❝ … how would you handling this missed time point at statistical analysis?


In the future describe how you would handle missing samples in an SOP or – much better IMHO – in the protocol. In the latter case any procedure is more transparent to an assessor (i.e., approved by the IEC and the local authority).

❝ 2. Will this data be included in statistical analysis? If yes, what is the handling procedure?


Depends on the intellectual horsepower you are willing/able to invest. Missings are generally not problematic if unambiguously located in the absorption1 or distribution/elimination2 phase. If the missing might be Cmax3 data imputation is problematic. Regulators don’t like PK modeling, so I have tried regressive splines with not very convincing results (references there). :-(

❝ 3. If the data is excluded, what kind of method would you use in calculate your AUC?

❝ a. Would you adapt interpolation method?

❝ b. If yes, would you adapt logarithmic interpolation or linear interpolation?


I always suggest to use the lin-up/log-down trapezoidal method. I think that Pharsight will make it even the default in the next version of Phoenix/WinNonlin. No need to impute data. Fine unless the missing might be Cmax… See also this presentation about different integration algos (especially slides 16–17 on what happens with missing values).

❝ 4. Is there any guideline or reference document we could follow?


Not that I recall.


Let’s denote denote the missing as Ci. Then if
  1. Ci+1 > Ci-1 ~ absorption
  2. Ci+1 < Ci-1 ~ distribution/elimination
  3. Ci+1 = Ci-1 ~ Cmax or pure chance? Any algo will come up with the linear interpolation Ci = Ci-1 + |(ti ti-1)/(ti+1 ti-1)|(Ci+1 Ci-1) – which likely underestimates the true value.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-22 10:25
(3963 d 13:30 ago)

@ Helmut
Posting: # 10614
Views: 17,330
 

 Lin-up/log-down trapezoidal avoids most trouble

Dear Helmut,

Great. However is there any special consideration when using lin-up/log-down trapezoidal rule to calculate AUCs if there is more than one absorption peak? Thanks in advance.

❝ ...

❝ I always suggest to use the lin-up/log-down trapezoidal method. I think that Pharsight will make it even the default in the next version of Phoenix/WinNonlin. No need to impute data. Fine unless the missing might be Cmax… See also this presentation about different integration algos (especially slides 16–17 on what happens with missing values).


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-22 22:25
(3963 d 01:31 ago)

@ yjlee168
Posting: # 10620
Views: 17,657
 

 Multiple peaks: fallback to linear trapezoidal

Hi Yung-jin!

❝ However is there any special consideration when using lin-up/log-down trapezoidal rule to calculate AUCs if there is more than one absorption peak?


Not that I am aware of. If values increase again the algo falls back to the linear trapezoidal. Try this clumsy code (a bimodal release formulation with two lag-times, two absorption rate constants, and a common elimination):

#set.seed(29081957) # keep this line to compare results only
require(truncnorm)
n       <- 48        # no of profiles
t       <- c(0,0.25,0.5,0.75,1,1.25,1.5,2,2.5,3,3.5,
             4,4.25,4.5,4.75,5,5.25,5.5,6,6.5,7,7.75,9,12,16,24)
D       <- 40
D.f     <- c(0.5, 0.5) # dose fractions of components
# Index ".d" denotes default values, ".c" CV
F.d     <- 0.5      # fraction absorbed (BA)
F.c     <- 0.1      # ±0.1, uniform
V.d     <- 1        # volume of distribution
V.c     <- 0.1      # CV 10%, lognormal
k01.1.d <- 1.4      # absorption rate constant (1st component)
k01.2.d <- 0.7      # absorption rate constant (2nd component)
k01.c   <- 0.2      # CV 20%, lognormal
k10.d   <- log(2)/3 # elimination rate constant (common to both components)
k10.c   <- 0.2      # CV 20%, lognormal
tlag1.d <- 0.5      # 1st lag-time
tlag1.l <- 0        # lower truncation 1
tlag1.u <- 1.5      # upper truncation 1
tlag2.d <- 4        # 2nd lag-time
tlag2.l <- 3        # lower truncation 2
tlag2.u <- 6        # upper truncation 2
tlag.c  <- 0.5      # CV 50%, truncated normal
tmax.th <- max((log(k01.1.d/k10.d)/(k01.1.d-k10.d)+tlag1.d),
               (log(k01.2.d/k10.d)/(k01.2.d-k10.d)+tlag2.d-tlag1.d))
Cmax.th <- F.d*D*D.f[1]/V.d*k01.1.d/(k01.1.d-k10.d)*
             (exp(-k10.d*(tmax.th-tlag1.d))-exp(-k01.1.d*(tmax.th-tlag1.d)))+
           F.d*D*D.f[2]/V.d*k01.2.d/(k01.2.d-k10.d)*
             (exp(-k10.d*(tmax.th-tlag2.d))-exp(-k01.2.d*(tmax.th-tlag2.d)))
x       <- seq(min(t), max(t), length.out=((max(t)-min(t))*4+1))
C1      <- F.d*D*D.f[1]/V.d*k01.1.d/(k01.1.d-k10.d)*
             (exp(-k10.d*(x-tlag1.d))-exp(-k01.1.d*(x-tlag1.d)))
C1[C1<0]<- 0        # set negative values of 1st component to zero
C2      <- F.d*D*D.f[2]/V.d*k01.2.d/(k01.2.d-k10.d)*
            (exp(-k10.d*(x-tlag2.d))-exp(-k01.2.d*(x-tlag2.d)))
C2[C2<0]<- 0        # set negative values of 2nd component to zero
C       <- C1+C2    # add components
AErr    <- 0.10     # analytical error CV 10%, lognormal
LLOQ    <- 0.01*Cmax.th # 1% of theoretical Cmax
split.screen(c(2,1))
for (s in 1:2) {
  screen(s)
  if (s == 1){
    plot(x, C, type="l", lwd=3, col="red", ylim=c(0, Cmax.th*1.75),
      main="linear trapezoidal", cex.main=1.1,
      xlab="time", ylab="concentration")
  } else {
    plot(x, C, type="l", lwd=3, col="red", ylim=c(0, Cmax.th*1.75),
      main="lin-up/log-down trapezoidal", cex.main=1.1,
      xlab="time", ylab="concentration")
  }
  rug(t, side=3)
  abline(v=0, lty=3)
  abline(h=0, lty=3)
}
AUCt.th <- NULL; AUCt.lin <- NULL; AUCt.linlog <- NULL
for (i in 1:n) {
  F      <- runif(n=1, min=F.d-F.c, max=F.d+F.c)
  V      <- rlnorm(n=1, meanlog=log(V.d)-0.5*log(V.c^2+1),
              sdlog=sqrt(log(V.c^2+1)))
  k01.1  <- rlnorm(n=1, meanlog=log(k01.1.d)-0.5*log(k01.c^2+1),
              sdlog=sqrt(log(k01.c^2+1)))
  k01.2  <- rlnorm(n=1, meanlog=log(k01.2.d)-0.5*log(k01.c^2+1),
              sdlog=sqrt(log(k01.c^2+1)))
  k10    <- rlnorm(n=1, meanlog=log(k10.d)-0.5*log(k10.c^2+1),
              sdlog=sqrt(log(k10.c^2+1)))
  tlag1  <- rtruncnorm(n=1, a=tlag1.l, b=tlag1.u, mean=tlag1.d, sd=tlag.c)
  tlag2  <- rtruncnorm(n=1, a=tlag2.l, b=tlag2.u, mean=tlag2.d, sd=tlag.c)
  C0.1   <- F*D*D.f[1]/V  # fictive iv concentration at t=0
  B1.1   <- C0.1*k01.1/(k01.1-k10)*exp(k10*tlag1)
  B2.1   <- C0.1*k01.1/(k01.1-k10)*exp(k01.1*tlag1)
  C1     <- B1.1*exp(-k10*t)-B2.1*exp(-k01.1*t)
  C1[C1<0] <- 0
  C0.2   <- F*D*D.f[2]/V  # fictive iv concentration at t=0
  B1.2   <- C0.2*k01.2/(k01.2-k10)*exp(k10*tlag2)
  B2.2   <- C0.2*k01.2/(k01.2-k10)*exp(k01.2*tlag2)
  C2     <- B1.2*exp(-k10*t)-B2.2*exp(-k01.2*t)
  C2[C2<0] <- 0
  C <- C1+C2
  for (j in 1:length(t)) { # analytical error (multiplicative 10% C)
    if (j==1)
      AErr1 <- rnorm(n=1, mean=0, sd=abs(C[j]*AErr))
    else
      AErr1 <- c(AErr1, rnorm(n=1, mean=0, sd=abs(C[j]*AErr)))
  }
  C <- C+AErr1 # add analytical error
  C[C<LLOQ] <- NA  # assign NAs to Cs below LLOQ
  tfirst <- (t[!is.na(C)])[1]              # 1st C >= LLOQ
  t0     <- t[length(head((t[is.na(C)])))] # last LLOQ before tfirst
  tlast  <- rev(t[!is.na(C)])[1]           # last C >= LLOQ
  screen(1, FALSE)
  plot(t, C, type="l", axes=FALSE,
    ylim=c(0, Cmax.th*1.75), xlab="", ylab="") # linear spaghetti plots
  screen(2, FALSE)
  for (j in 1:(length(t)-1)) {
    if (is.na(C[j]) == FALSE & is.na(C[j+1]) == FALSE) {
      if (C[j+1] >= C[j]) {
        lines(x=c(t[j],t[j+1]), y=c(C[j],C[j+1])) # increasing segments
      } else {
        tk <- NULL; Ck <- NULL # interpolated segments
        for (k in seq(t[j], t[j+1], length.out=10)) {
          tk <- c(tk, k)
          Ck <- c(Ck, exp(log(C[j])+abs((k-t[j])/(t[j+1]-t[j]))*(log(C[j+1])-log(C[j]))))
        }
        lines(tk, Ck) # decreasing segments
      }
    }
  }
  Cwo.err1 <- function(x) {B1.1*exp(-k10*x)-B2.1*exp(-k01.1*x)}
  Cwo.err2 <- function(x) {B1.2*exp(-k10*x)-B2.2*exp(-k01.2*x)}
  AUCt.1 <- integrate(Cwo.err1, lower=tlag1, upper=max(t))
  AUCt.2 <- integrate(Cwo.err2, lower=tlag2, upper=max(t))
  AUCt <- AUCt.1$value+AUCt.2$value # AUCt without analytical error
  AUCt.th <- c(AUCt.th, AUCt)
  pAUC.lin <- 0 # linear trapezidal rule
  for (j in seq(which(t == t0), which(t == tlast)-1)){
    if (t[j] == t0) { # 1st triangle
      pAUC.lin <- pAUC.lin + 0.5*(t[j+1]-t[j])*C[j+1]
    } else {
      pAUC.lin <- pAUC.lin + 0.5*(t[j+1]-t[j])*(C[j+1]+C[j])
    }
  }
  AUCt.lin <- c(AUCt.lin, pAUC.lin)
  pAUC.linlog <- 0 # lin-up/log-down trapezidal rule
  for (j in seq(which(t == t0), which(t == tlast)-1)){
    if (t[j] == t0) { # 1st triangle
      pAUC.linlog <- pAUC.linlog + 0.5*(t[j+1]-t[j])*C[j+1]
      } else {
      if (C[j+1] >= C[j]) { # linear trap.
        pAUC.linlog <- pAUC.linlog + 0.5*(t[j+1]-t[j])*(C[j+1]+C[j])
      } else { # loglinear trap.
        pAUC.linlog <- pAUC.linlog + (t[j+1]-t[j])*(C[j+1]-C[j])/(log(C[j+1]/C[j]))
      }
    }
  }
  AUCt.linlog <- c(AUCt.linlog, pAUC.linlog)
}
close.screen(all = TRUE)
par(ask=T)
RE.AUCt.lin <- 100*(AUCt.lin-AUCt.th)/AUCt.th
RE.AUCt.linlog <- 100*(AUCt.linlog-AUCt.th)/AUCt.th
boxplot(RE.AUCt.lin, RE.AUCt.linlog, ylab="% RE",
  ylim=c(-max(RE.AUCt.lin, RE.AUCt.linlog)*1.1, max(RE.AUCt.lin, RE.AUCt.linlog)*1.1),
  main="Bias of algos", cex.main=1.1,
  names=c("linear", "lin-up/log-down"))
abline(h=0, lty=3)
cat("\nCVtotal; linear trap.:", round(sd(AUCt.lin)/mean(AUCt.lin),4),
  "– lin-up/log-down trap.:", round(sd(AUCt.linlog)/mean(AUCt.linlog),4), "\n\n")
wilcox.test(RE.AUCt.lin, mu=0, conf.int=T)
wilcox.test(RE.AUCt.linlog, mu=0, conf.int=T)
wilcox.test(RE.AUCt.lin, RE.AUCt.linlog, conf.int=T)
par(ask=F)


No big deal, but the lin-up/log-down should perform better if there are missings.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-22 23:57
(3962 d 23:59 ago)

@ Helmut
Posting: # 10622
Views: 17,343
 

 Multiple peaks: fallback to linear trapezoidal

Dear Helmut,

Thank you for your messages and codes. I tried to run your codes without luck, but I will try it later. I will change the way of AUC calculation with bear later. I remembered that I read something about method comparisons for trapezoidal rule in one textbook, but just cannot remember which book.

❝ Not that I am aware of. If values increase again the algo falls back to the linear trapezoidal. Try this clumsy code (a bimodal release formulation with two lag-times, two absorption rate constants and a common elimination):


...

  B1.1   <- C0.1*k01.1/(k01.1-K)*exp(k10*tlag1)

  B2.1   <- C0.1*k01.1/(k01.1-K)*exp(k01.1*tlag1)


Here I got Error: object 'K' not found and some error messages after this.

One more question, in your above presentation (slide#17[edited]; missing data occurred in R at 12h), it was a lin-up/log-down plot. As you said that no need to do data imputation if using lin-up/log-down. Do you do any line smoothing with your data first? Otherwise, the line should not be connected in that way. It should be connected as the line that the arrow points. Do I miss anything here?

[image]

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-23 01:25
(3962 d 22:31 ago)

@ yjlee168
Posting: # 10623
Views: 17,477
 

 No interpolation

Hi Yung-jin,

❝ I will change the way of AUC calculation with bear later.


Relax. Would make only sense if the user has a choice. Some people still get crazy if one tries anything else than the linear trapezoidal.
Don’t use my code without improving it! In its current form it will not work with NAs and/or isolated BQLs within the profile.

❝ I remembered that I read something about method comparisons for trapezoidal rule in one textbook, but just cannot remember which book.


Yes. I think I have posted even references somewhere. To lazy to search now.

❝ ❝ ...

❝ ❝   B1.1   <- C0.1*k01.1/(k01.1-K)*exp(k10*tlag1)

❝ ❝   B2.1   <- C0.1*k01.1/(k01.1-K)*exp(k01.1*tlag1)

❝ Here I got Error: object 'K' not found and some error messages after this.


Sorry, I made a copy/paste error. Should work now (use k10 instead of K).

❝ One more question, in your above presentation (slide#17[edited]; missing data occurred in R at 12h), it was a lin-up/log-down plot.


Nope. 16/17 are linear plots. In slide 16 I connected all points by straight lines (what the linear trap. would do), and in slide 17 by straight lines if ascending and by an exponential if descending.

❝ As you said that no need to do data imputation if using lin-up/log-down. Do you do any line smoothing with your data first?


No.

❝ Otherwise, the line should not be connected in that way. It should be connected as the line that the arrow points.


Why?

❝ Do I miss anything here?


Yes, you do. ;-) If the 12 h sample is missing the linear trapezoidal works as if we would impute a value by linear interpolation. The lin-up/log-down works as if we would impute a value by logarithmic interpolation. But the trick is that nothing has to be imputed in the lin-up/log-down at all!
  t    C
  8  33.17
 12  miss.
 16   7.86


Interpolated values for imputation:
linear: C12*=33.17+|(12-8)/(16-8)|(7.86-33.17)=20.52
linlog: C12*=exp(log(33.17)+|(12-8)/(16-8)|log(7.86/33.17))=16.15


Calculation based on imputed values:
linear/lin imp.: pAUC8-16=0.5(12- 8)(33.17+20.52)+
                          0.5(16-12)(20.52+ 7.86)=164.1
linear/log imp.: pAUC8-16=0.5(12- 8)(33.17+20.52)+
                          0.5(16-12)(20.52+ 7.86)=146.7
linlog/log imp.: pAUC8-16=(12- 8)(16.15-33.17)/log(16.15/33.17)+
                          (16-12)( 7.86-16.15)/log( 7.86/16.15)=140.6


Direct calculation without imputation:
linear: pAUC8-16=0.5(16-8)(33.17+7.86)=164.1
linlog: pAUC8-16=(16-8)(7.86-33.17)/log(7.86/33.17)=140.6


BTW, the theoretical pAUC8-16 based on the model is 140.62…

To make a long story short: Only if one has to use the linear trapezoidal (for whatever wacky reasons) I would impute an estimate based on the log/linear interpolation (second variant above). Otherwise the bias is substantial. Q.E.D.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-23 13:43
(3962 d 10:12 ago)

@ Helmut
Posting: # 10626
Views: 17,010
 

 No interpolation

Dear Helmut,

Thanks for your messages, and the codes works well now.

❝ ...

❝ Relax. Would make only sense if the user has a choice. Some people still get crazy if one tries anything else than the linear trapezoidal.


Yes, that's my concern too. I should provide users options for this step in NCA. It should be easy to implement such a function with bear at the moment.

❝ Don’t use my code without improving it! In its current form it will not work with NAs and/or isolated BQLs within the profile.


Sure and thanks for your reminding.

❝ ...

❝ ❝ Otherwise, the line should not be connected in that way. It should be connected as the line that the arrow points.


❝ Why?


Sorry, I thought it was using lin-up/log-down algo and plots in slide#17.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-23 17:15
(3962 d 06:41 ago)

@ yjlee168
Posting: # 10630
Views: 17,325
 

 Spaghetti & other pasta

Dear Yung-jin,

❝ ❝ ❝ Otherwise, the line should not be connected in that way. It should be connected as the line that the arrow points.

❝ ❝ Why?

❝ Sorry, I thought it was using lin-up/log-down algo and plots in slide#17.


We are all guilty in not representing in spaghettis plots the method we use for calcu­lat­ing the AUC. For simplicity we connect point with straight lines, regardless whether we use a linear or semilogarithmic scale (all left panels below).
If we use the linear trapezoidal, the semilog plot is wrong. What we are actually cal­cu­lating would be reflected in the lower right.

[image]

Weird, isn’t it? But then hopefully everybody would realize that the algo overestimates AUCs in the distribution / elimination phase.

On the other hand if we use the lin-up/log-down trapezoidal the linear plot is wrong; should use the upper right one instead.

[image]

Have you ever seen something like this before? I didn’t. Would be a tough job to come up with such a setup in commercial software. :-D

PS: Very rarely you find publications where in mean plots the points are not connected by lines at all.
PPS: What do you think about adding rug(t, side=3) to plots in bear?

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-23 22:47
(3962 d 01:08 ago)

@ Helmut
Posting: # 10632
Views: 17,105
 

 Spaghetti & other pasta

Dear Helmut,

Wow! That really quite surprises me. I know that all errors may result from "the approximation of trapezoid" (or assumption of trapezoid area between two points) when calculating AUC. Yes, right now we just use line connection between points to plot linear/semilog graphs. From your plots, it looks like lin-up/log-down method is more close to actual AUC. In this case, it does no matter what there is any missing data or not. Lin-up/log-down is more accurate than linear, isn't?

❝ If we use the linear trapezoidal, the semilog plot is wrong. What we are actually calculating would be reflected in the lower right.


[image]


❝ Weird, isn’t it? But then hopefully everybody would realize that the algo overestimates AUCs in the distribution / elimination phase.


Weird indeed but very convincible. Great demo. How do you plot the low-right graph? Is it possible to overlap both the lower-left and the lower-right graph as one plot? That make differences more clear.

❝ Have you ever seen something like this before? I didn’t. Would be a tough job to come up with such a setup in commercial software. :-D


Me either. Why do you think that it would be a tough job for commercial software?

❝ PPS: What do you think about adding rug(t, side=3) to plots in bear?


Should be do-able in bear. But I don't quite understand why you want to add rug() plot? We don't have many data points in each subject (usually less than 50's in each period) for a BE/BA study.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-23 23:29
(3962 d 00:27 ago)

@ yjlee168
Posting: # 10633
Views: 17,219
 

 Spaghetti & other pasta

Dear Yung-jin!

❝ […] From your plots, it looks like lin-up/log-down method is more close to actual AUC. In this case, it does no matter what there is any missing data or not. Lin-up/log-down is more accurate than linear, isn't?


I think so. Whatever PK a drug follows, distribution/elimination likely is first order (well, leaving :party: and Michaelis-Menten aside).

❝ Weird indeed but very convincible. […] How do you plot the low-right graph? Is it possible to overlap both the lower-left and the lower-right graph as one plot?


Sure. Just interpolate linear within each decreasing segment. Alternatively to lines 12–17 you could set up a function and use curve(). Was too lazy.
t <- c(0,0.5,0.75,1,1.5,2,2.5,3,4,6,8,16,24)
C <- c(0,69.59,84.64,92.61,96.73,93.45,87.51,
       80.85,68.02,47.54,33.17,7.86,1.86)
plot(t, C, type="p", ylim=c(1, max(C)*1.1), log="y", las=1,
  main="semilogarithmic", cex.main=1.1,
  xlab="time", ylab="concentration")
lines(t, C, col="blue")
for (j in 1:(length(t)-1)) {
  if (C[j+1] >= C[j]) {
    lines(x=c(t[j], t[j+1]), y=c(C[j], C[j+1]), col="blue")
  } else {
    tk <- NULL; Ck <- NULL
    for (k in seq(t[j], t[j+1], length.out=10)) {
      tk <- c(tk, k)
      Ck <- c(Ck, C[j]+abs((k-t[j])/(t[j+1]-t[j]))*(C[j+1]-C[j]))
    }
    lines(tk, Ck, col="red")
  }
}
rug(t, side=1)
legend("topright",
  c("trapezoidal rules", "linear", "lin-up/log-down"),
  lty=c(NA, 1, 1), col=c(NA, "red", "blue"))


[image]

If you want to use the code in a linear plot, interpolate logarithmically, i.e., change line 15 to
Ck <- c(Ck, exp(log(C[j])+abs((k-t[j])/(t[j+1]-t[j]))*(log(C[j+1])-log(C[j]))))

❝ Me either. Why do you think that it would be a tough job for commercial software?


I have no experiences with SAS, but in Phoenix/WinNonlin it is tricky indeed.

❝ […] But I don't quite understand why you want to add rug() plot? We don't have many data points in each subject (usually less than 50's in each period) for a BE/BA study.


It was just an idea to display the sampling schedule. Forget it. ;-)

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-24 00:46
(3961 d 23:10 ago)

@ Helmut
Posting: # 10634
Views: 16,991
 

 Spaghetti & other pasta

Dear Helmut,

Thanks again. Now I think that I have fully understood the differences between the linear and the lin-up/log-down method for AUC calculations. Actually, lin-up/log-down method is still a kind of data smoothing approach. It is applied to calculate more accurate AUC than the linear one. Some people should not get crazy when using lin-up/log-down if they can see the differences here. Just because the line is too smooth to be true? That may scare people a little bit. That's why I asked if you used some kinds of data smoothing methods in previous message. The key point is that involved calculation step still uses the raw data to calculate the interpolated data. As we can see, the errors between linear and lin-up/log-down become significantly apparent, especially when sampling time between two data points is getting greater. Therefore, when there is missing data occurring, the lin-up/log-down method shows more accurate than the linear.

❝ [...]


❝ I have no experiences with SAS, but in Phoenix/WinNonlin it is tricky indeed.


should be no problem at all. They are professional after all.

❝ [...]


❝ It was just an idea to display the sampling schedule. Forget it. ;-)


I see. However, it could be misunderstood or confused by people to interpret it (the rug) as a 'log scale' with x axis.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-24 03:24
(3961 d 20:32 ago)

@ yjlee168
Posting: # 10635
Views: 16,985
 

 Spaghetti & other pasta

Dear Yung-jin!

❝ Now I think that I have fully understood the differences between the linear and the lin-up/log-down method for AUC calculations.


Hhm.

❝ Actually, lin-up/log-down method is still a kind of data smoothing approach.


No, it isn’t.

❝ The key point is that involved calculation step still uses the raw data to calculate the interpolated data.


There is no interpolation done at all. Have a look at my example in the previous post. I gave the interpolation/imputation formulas only as a hint what might happen if we need an intermediate data point. If we want to calculate the area in the interval [t2, t2] no interpolation, smoothing, whatsoever is done. Both methods are straightforward:

linear: AUCt1→t2=0.5(t2–t1)(C2+C1)
linlog: AUCt1→t2=(t2–t1)(C2–C1)/log(C2/C1)


❝ As we can see, the errors between linear and lin-up/log-down become significantly apparent, especially when sampling time between two data points is getting greater.


Right.

❝ Therefore, when there is missing data occurring, the lin-up/log-down method shows more accurate than the linear.


Xactly. See also Nathan’s PK/PD Blog.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-24 11:17
(3961 d 12:39 ago)

@ Helmut
Posting: # 10638
Views: 17,073
 

 Spaghetti & other pasta

Dear Helmut,

❝ ❝ Actually, lin-up/log-down method is still a kind of data smoothing approach.


❝ No, it isn’t.


You are right. It isn't. Sorry about this.

❝ [...] I gave the interpolation/imputation formulas only as a hint what might happen if we need an intermediate data point. If we want to calculate the area in the interval [t2, t2] no interpolation, smoothing, whatsoever is done. Both methods are straightforward:


I see. Thank you for your detailed explanation once again.

linear: AUCt1→t2=0.5(t2–t1)(C2+C1)

linlog: AUCt1→t2=(t2–t1)(C2–C1)/log(C2/C1)


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Ken Peh
★    

Malaysia,
2013-05-30 21:55
(3955 d 02:01 ago)

@ Helmut
Posting: # 10688
Views: 16,869
 

 No interpolation

Dear Helmut,

❝ Interpolated values for imputation:

linear: C12*=33.17+|(12-8)/(16-8)|(7.86-33.17)=20.52


ok.

linlog: C12*=exp(log(33.17)+|(12-8)/(16-8)|log(7.86/33.17))=16.15


How did you get the readings of 16.15? I tried to follow your equation but I could not get the answer.:confused:

❝ Calculation based on imputed values:

linear/lin imp.: pAUC8-16=0.5(12- 8)(33.17+20.52)+

                          0.5(16-12)(20.52+ 7.86)=164.1

linear/log imp.: pAUC8-16=0.5(12- 8)(33.17+20.52)+

                          0.5(16-12)(20.52+ 7.86)=146.7


linear/lin imp and linear/log imp have the same data set. Why the output is different ?:confused:

linlog/log imp.: pAUC8-16=(12- 8)(16.15-33.17)/log(16.15/33.17)+

                          (16-12)( 7.86-16.15)/log( 7.86/16.15)=140.6


❝ Direct calculation without imputation:

linear: pAUC8-16=0.5(16-8)(33.17+7.86)=164.1


ok.

linlog: pAUC8-16=(16-8)(7.86-33.17)/log(7.86/33.17)=140.6


Similarly, I could not get the answer of 140.61. How did you calculate ?:confused:

Please kindly elaborate.

Thank you.

Ken
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-31 00:38
(3954 d 23:18 ago)

@ Ken Peh
Posting: # 10690
Views: 17,051
 

 Different algos!

Hi Ken!

❝ ❝ linlog: C12*=exp(log(33.17)+|(12-8)/(16-8)|log(7.86/33.17))=16.15

❝ How did you get the readings of 16.15? I tried to follow your equation but I could not get the answer.


R
round(exp(log(33.17)+abs((12-8)/(16-8))*log(7.86/33.17)), 2)
[1] 16.15
[image]
SAS
x=round(exp(log(33.17)+abs((12-8)/(16-8))*log(7.86/33.17)), 2);
16.15
[image]
Excel / OO Calc / Gnumeric
=ROUND(EXP(LN(33.17)+ABS((12-8)/(16-8))*LN(7.86/33.17)), 2)
16.15
[image]

❝ linear/lin imp and linear/log imp have the same data set. Why the output is different ?


Since I’m using different algorithms. Once more: Forget linear interpolation in the descending part of the profile.

[image]

Have you ever asked yourself why the arithmetic and geometric means of the same data are different?
Compare (7.86+33.17)*0.5
with     (7.86*33.17)^0.5
Are these values vaguely reminiscent of something you have seen before? :lookaround:

❝ ❝ linlog: pAUC8-16=(16-8)(7.86-33.17)/log(7.86/33.17)=140.6

❝ Similarly, I could not get the answer of 140.61.


round((16-8)*(7.86-33.17)/log(7.86/33.17), 1)
[1] 140.6
[image]
=ROUND((16-8)*(7.86-33.17)/LN(7.86/33.17), 1)
140.6
[image]

❝ How did you calculate ?


See above. And you?

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Ken Peh
★    

Malaysia,
2013-06-03 19:52
(3951 d 04:04 ago)

@ Helmut
Posting: # 10717
Views: 16,664
 

 Different algos!

Dear Helmut,

Thank you very much.

❝ ❝ How did you calculate ?


I know where my mistake is. It is the positive and negative signs. Brain is jammed up. Need more practice.

May I know which software in R you refer to that can be used to calculate the data ?

Regards,
Ken
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-06-03 20:11
(3951 d 03:45 ago)

@ Ken Peh
Posting: # 10718
Views: 16,864
 

 Calculate what?

Hi Ken!

❝ May I know which software in R you refer to that can be used to calculate the data ?


… to calculate what? For imputation use my code above (simply give the time­points in the t-vector and concentrations in C). For the lin-up/log-down trape­zoidal rule in the calculation of AUC see the end of my code in this post.
I guess it’s easier to wait for the next version of bear.* ;-)


    See there for the new version.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-25 17:09
(3960 d 06:47 ago)

@ yjlee168
Posting: # 10656
Views: 17,274
 

 Lin-up/log-down trapezoidal example

Dear Yung-jin,

❝ […] is there any special consideration when using lin-up/log-down trapezoidal rule to calculate AUCs if there is more than one absorption peak?


I checked one study of bimodal release formulations. Imbalanced, n=15. I evaluated AUCt based on the linear trapezoidal and the lin-up/log-down method. Then I introduced a failure of the centrifuge resulting in a complete loss of 12 h samples in the second period. Let’s see:

                               PE       90% CI     %CVintra  %CVinter
complete, linear              99.89  95.23–104.78    7.38     33.75
complete, lin-up/log-down     99.78  95.10–104.69    7.42     33.40
incomplete, linear           100.08  95.64–104.72    7.00     33.84
incomplete, lin-up/log-down   99.93  95.37–104.72    7.23     33.39

With the linear trapezoidal results change whereas with the lin-up/log-down results are affected to a lesser degree. I didn’t check but expect to get identical results in a balanced case.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2013-05-25 21:45
(3960 d 02:11 ago)

@ Helmut
Posting: # 10658
Views: 16,902
 

 Lin-up/log-down trapezoidal example

Dear Helmut,

I think the differences at elimination phase between linear and log-down may results from the degree of concavity ([C1-C2]) within a know time interval (t2-t1). I also did a simple test with two different scenarios:

** if t2-t1 = 12 h; C1 = 4.1, C2 =1.1 (the curve is more concave between t1 and t2)

  AUC(t1-t2) linear = 31.2; log-down = 27.362 (log-down should be more accurate in this situation)
 
** if t2-t1 = 12 h; C1 = 4.1, C2 =4.05 (the curve is less concave between t1 and t2; AUC(t1-t2) should be more close to real trapezoidal shape)

  AUC(t1-t2) linear = 48.9; log-down = 48.9 (before truncated was 48.89939) (the linear is supposed to be better than log-down; however log-down still keeps the accuracy as good as linear and almost as same as the linear.)
 

Therefore, you're right about log-down method. The log-down can be set as default method for AUC calculation safely, accurately and reliably.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Ohlbe
★★★

France,
2013-05-20 23:45
(3965 d 00:11 ago)

@ gracehung
Posting: # 10609
Views: 17,523
 

 handling of missing data

Dear Gracehung,

❝ As in BE study, One sample’s actual sampling time was miss recorded. However as according to protocol, this sample was analyzed by our analytical laboratory.


I understand you question differently from Helmut. My understanding is that the sample was collected, but the phlebotomist failed to record the actual blood sampling time. In such a case I would just file a protocol deviation and re-train the staff, full stop. Even if there were a sampling time deviation, the influence on the AUC 90 % CI is so limited that it would only change the outcome of the trial in a very borderline situation.

Regards
Ohlbe
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-21 15:53
(3964 d 08:03 ago)

@ Ohlbe
Posting: # 10610
Views: 17,192
 

 Sorry

Hi Ohlbe,

❝ I understand you question differently from Helmut.


Though I understood it correctly (see the start of my post), I am notorious for answering unasked questions as well. :-D

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
gracehung
☆    

Taiwan,
2013-05-23 03:06
(3962 d 20:50 ago)

@ Helmut
Posting: # 10624
Views: 17,239
 

 Sorry

Dear all:

Thanks for your reply!
The dilemma is we have the concentration data but don't have the actual sampling time. If we treat it as missing data, don't use it in bioequivalence evaluation, would assessor think we are hiding data? Should we include the data in PK analysis and use schedule time to be the "best guess"? One literature suggest to use log linear interpolation (at elimination phase) to minimize the bias (also as your reply). But someone said the assessor will think you are manipulating data. :confused: which one would be the best choice?

Grace
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-05-23 03:35
(3962 d 20:21 ago)

@ gracehung
Posting: # 10625
Views: 17,205
 

 Uncertain time point

Hi Grace!

❝ If we treat it as missing data, don't use it in bioequivalence evaluation, would assessor think we are hiding data?


Duno. But if you treat the value as missing – depending on the algo – you might get a sub­stantially biased AUC (see my conversation with Yung-jin above).

❝ Should we include the data in PK analysis and use schedule time to be the "best guess"?


Definitely (as Ohlbe suggested). I once had to assess deficiencies where in two studies the scheduled time points (instead of the actual ones) were used. Impact on the AUC comparison was in the second decimal place. And you have only one (!) value in an entire study…
BTW, scheduled time points were my standard method in hundreds (!) of studies. So what?

❝ One literature suggest to use log linear interpolation (at elimination phase) to minimize the bias (also as your reply).


Yes, but only if the value is missing. You know the value – only its time point is uncertain.

❝ But someone said the assessor will think you are manipulating data.


I would not use an estimate. Use the measured one at the scheduled time point. Report the deviation from the protocol, don’t panic, & sleep well.

Another idea – only if you really have nothing better to do: Most protocols require a justi­fi­cation of a time deviation to be recorded in CRF only if a certain time interval is exceeded (f.i. ±2’ in the very early phase, ±5’ later on, and ±15’ for t ≥24 h). Let’s assume only the actual sampling was not recorded, but still was within the allowed interval (otherwise we have two errors instead of one). Set up two data sets: One assuming the sample was taken at the maxi­mum allowed negative deviation (too early) and the other one with the max. positive deviation (too late). Go ahead with the BE and assess the impact on the outcome.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,636 registered users;
90 visitors (0 registered, 90 guests [including 6 identified bots]).
Forum time: 22:56 CET (Europe/Vienna)

Nothing shows a lack of mathematical education more
than an overly precise calculation.    Carl Friedrich Gauß

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5