Achievwin
★★  

US,
2023-07-05 20:38
(437 d 01:44 ago)

Posting: # 23655
Views: 3,303
 

 any SAS or R-macro for 3 period two sequence RSABEIs [RSABE / ABEL]

Post reply
Hi

Is there any SAS or R-Macro or program for performing RSABR in a 3-period 2-sequnce balanced replicate design? or how a doubt given below can be clarified?

A 3-period, 2-sequence (TRR or RTT) is proposed assuming it is a balanced design! how to clarify the reference-scaled average bioequivalence approach will be applied to the proposed sequences.

Thinking Rationale: was that the design was preferred in the current statistical approaches guidance current SAS macros ignore any missing sequence and provide RSABE conclusions with existing macros and secondly any unknown carry-over effect is best addressed in the current sequence by balance effect.

Any input is appreciated.
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2023-07-06 10:46
(436 d 11:35 ago)

@ Achievwin
Posting: # 23657
Views: 2,951
 

 Use the one for TRTR | RTRT

Post reply
Hi Achievwin,

❝ Is there any SAS or R-Macro or program for performing RSABR in a 3-period 2-sequnce balanced replicate design?

Use the SAS code for the 4-period 2-sequence replicate design given in the Progesterone and ANDA guidances. Since you have also two sequences, only trivial modifications to select the treatments in the first part of the code. My guess for TRR | RTT (I don’t speak SAS):

data test1;
  set test;
  if (seq=1 and per=1);
  lat1t=lauct;
run;
data test2;
  set test;
  if (seq=2 and per=2) or (seq=2 and per=3);
  lat2t=lauct;
run;
data ref1;
  set ref;
  if (seq=1 and per=2) or (seq=1 and per=3);
  lat1r=lauct;
run;
data ref2;
  set ref;
  if (seq=2 and per=1);
  lat2r=lauct;
run;

Also doable in R (I think Detlew posted sumfink back in the day; too lazy to search).

❝ or how a doubt given below can be clarified?

❝ A 3-period, 2-sequence (TRR or RTT) is proposed assuming it is a balanced design! how to clarify the reference-scaled average bioequivalence approach will be applied to the proposed sequences.

I guess you are referring to the draft guidance of December 2022. You can use also TRT | RTR (see page 7). Accepting 3-period 2-sequence full replicate designs is a step forward because when using the SAS code for the lousy TRR | RTR | RRT (ABE if swR < 0.294), the optimizer may fail to converge (see this article).

❝ Thinking Rationale: was that the design was preferred in the current statistical approaches guidance current SAS macros ignore any missing sequence and provide RSABE conclusions with existing macros and secondly any unknown carry-over effect is best addressed in the current sequence by balance effect.

Not sure whether I understand you correctly. What is stated in the guidances about carry-over is nothing more than :blahblah:. Unequal carry-over – which would bias the treatment effect – cannot be ‘corrected’ statistically and only avoided by sufficiently long washouts (see this article).
—
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
 
mittyri
★★  

Russia,
2023-07-06 16:07
(436 d 06:15 ago)

@ Helmut
Posting: # 23658
Views: 2,888
 

 RSABE R coded by Detlew

Post reply
Hi Achievwin and Helmut,

❝ Also doable in R (I think Detlew posted sumfink back in the day; too lazy to search).


here it is. The best ever. FDA guys were happy.
—
Kind regards,
Mittyri
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2023-07-07 16:41
(435 d 05:41 ago)

@ mittyri
Posting: # 23660
Views: 2,827
 

 ‘dlat’ for 3-period full replicate?

Post reply
Hi mittyri,

THX for discovering this gem!

I have a problem (both in Phoenix/WinNonlin, R, and – likely – in SAS as well). To calculate swR the FDA wants a linear model (in SAS-lingo proc glm) of repeated reference data with subject as a factor. However, in a 3-period full replicate they are all in one sequence and the model collapses. Now what?

If you want to play with data sets
  1. https://bebac.at/downloads/ds01.csv (TRTR|RTRT; the EMA’s reference data set I)
  2. https://bebac.at/downloads/ds03.csv (TRT|RTR; period 4 of #1 removed)
—
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
 
mittyri
★★  

Russia,
2023-07-07 18:04
(435 d 04:18 ago)

@ Helmut
Posting: # 23661
Views: 2,871
 

 ‘dlat’ for 3-period full replicate?

Post reply
Hi Helmut!

❝ I have a problem (both in Phoenix/WinNonlin, R, and – likely – in SAS as well). To calculate swR the FDA wants a linear model (in SAS-lingo proc glm) of repeated reference data with subject as a factor. However, in a 3-period full replicate they are all in one sequence and the model collapses. Now what?

you discovered another point for a lecture ;-)
So the problem is here
proc glm data=scavbe;
  class seq;
  model dlat=seq;
  ods output overallanova=dglm1;
  ods output NObs=dglm3;
  title1 'scaled average BE';
run;

I would look at the original formula

\(s_\text{wR}^2=\Large\frac{\sum\limits_{i=1}^{m}\sum_\limits{j=1}^{n_i}\left(D_{ij}-\overline{D}_{i\,\cdot}\right)^2}{2\left(n-m\right)}\)

So it is impossible to get D for `TRT` sequence. Could we just omit one sequence here?
—
Kind regards,
Mittyri
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2023-07-07 22:06
(435 d 00:16 ago)

@ mittyri
Posting: # 23662
Views: 3,010
 

 Likely wrong workaround?

Post reply
Hi mittyri,

❝ I would look at the original formula

❝      \(s_\text{wR}^2=\Large\frac{\sum\limits_{i=1}^{m}\sum_\limits{j=1}^{n_i}\left(D_{ij}-\overline{D}_{i\,\cdot}\right)^2}{2\left(n-m\right)}\)

❝ So it is impossible to get D for `TRT` sequence. Could we just omit one sequence here?


A first attempt. [image].script at the end.
  • The 4-period full replicate
    swR       = 0.45168 (>= 0.294, RSABE acceptable)
    CVwR      =  47.57%
    PE        = 115.46% (pass)
    90% CI    = 106.39% - 125.31% (informative only)
    critbound = -0.09454 (<= 0, pass)
    aggregate = pass
    If I would apply the workaround though not necessary:
    swR       = 0.45874
    CVwR      =  48.40%
    In Phoenix/WinNonin (cross-validated with SAS):
    swR       = 0.44645 (>= 0.294, RSABE acceptable)
    CVwR      =  46.96%
    PE        = 115.46% (pass)
    90% CI    = 106.39% - 125.31% (informative only)
    critbound = -0.09208 (<= 0, pass)
    aggregate = pass
    Why? BTW, my code for ABE is seemingly completely wrong.

  • The 3-period full replicate
    swR       = 0.53372 (>= 0.294, RSABE acceptable)
    CVwR      =  57.41%
    PE        = 124.52% (pass)
    90% CI    = 113.72% - 136.34% (informative only)
    critbound = -0.096683 (<= 0, pass)
    aggregate = pass
    Note: Based on workaround; questionable!
Any ideas/suggestions?

# for Detlew’s original code see
# https://forum.bebac.at/forum_entry.php?id=15490
# libraries for the mixed-effects model of ABE
suppressMessages(library(lmerTest)) # requires also 'lme4' and 'Matrix'
alpha      <- 0.05 # level of the tests
ds         <- "per4full" # EMA data set I (TRTR|RTRT)
ds         <- "per3full" # period 3 of I removed (TRT|RTR)
if (ds == "per4full") {
  f <- "https://bebac.at/downloads/ds01.csv"
} else {
  f <- "https://bebac.at/downloads/ds03.csv"
}
data        <- read.csv(file = f, stringsAsFactors = FALSE)
names(data) <- c("subj", "per", "seq", "treat", "pk")
# sort by subject, treatment and period to get replicate number
# this works in contrast to the SAS code also for
# other sequences than TRTR|RTRT
data        <- data[order(data$subj, data$treat, data$per), ]
data$repl   <- 1
for (i in c(2:nrow(data))) {
  if (data$subj[i] != data$subj[i-1] | data$treat[i] != data$treat[i-1]) {
    data$repl[i] <- 1
  } else {
    data$repl[i] <- data$repl[i-1] + 1
  }
}
# make a column with treatment and repl, i.e T1,T2, R1,R2
data$code2 <- paste0(data$treat, data$repl)
data       <- data[order(data$subj, data$per), ]
# reshape to wide with colums subj, seq, pk.R1, pk.R2, pk.T1, pk.T2
# and log-transform pk
data_r     <- reshape(data = data, direction = "wide", v.names = "pk",
                      idvar = c("subj", "seq"), timevar = "code2",
                      drop = c("per", "treat", "repl"))
data_r[, 3:6]      <- log(data_r[, 3:6])
names(data_r)[3:6] <- sub("pk.", "logpk.", names(data_r))[3:6]
# now T-R analysis, ilat in the SAS code
# modified for TRT|RTR and TTR|RRT
data_r$TR <- NA
if (ds == "per4full") {
  data_r$TR <- (data_r$logpk.T1 + data_r$logpk.T2 -
                data_r$logpk.R1 - data_r$logpk.R2) / 2
} else {
  for (i in 1:nrow(data_r)) { # clumsy
    if (data_r$seq[i] == "RTR" | data_r$seq[i] == "RRT")
      data_r$TR[i] <- data_r$logpk.T1[i] -
                      (data_r$logpk.R1[i] + data_r$logpk.R2[i]) / 2
    if (data_r$seq[i] == "TRT" | data_r$seq[i] == "TTR")
      data_r$TR[i] <- (data_r$logpk.T1[i] + data_r$logpk.T2[i]) / 2 -
                       data_r$logpk.R1[i]
  }
}
# now get rid of subjects with missing periods (TR == NA)
data_r    <- data_r[complete.cases(data_r$TR), ]
# ilat analysis, ANOVA with seq as effect
data_r$seq <- as.factor(data_r$seq)
# with standard contr.treatment we get not the wanted intercept!
# with that the intercept is intercept+seq1
oc        <- options("contrasts") # save options
options(contrasts = c("contr.sum", "contr.poly"))
m1        <- lm(TR ~ seq, data = data_r)
# intercept is estimate of µt-µR
est       <- coef(m1)[[1]]
pe1       <- exp(est)
ci1       <- as.numeric(confint(m1, "(Intercept)", level = 1 - 2 * alpha))
dfTR      <- m1$df.residual
# stderr, "the unknown x" for unbiased estimate of (µT-µR)^2
stderr    <- summary(m1)$coefficients["(Intercept)","Std. Error"]
# linearized scABE criterion component x
x         <- est^2 - stderr^2
boundx    <- max(abs(ci1))^2
ci1       <- exp(ci1)
data_r$RR <- data_r$logpk.R1 - data_r$logpk.R2
# now get rid of subjects where R is not repeated (RR == NA)
data_r    <- data_r[complete.cases(data_r$RR), ]
if (ds == "per4full") {
  # now the equivalent of SAS code for R-R, dlat analysis
  m2      <- lm(RR ~ seq, data = data_r)
  dfRR    <- m2$df.residual
  s2wR    <- summary(m2)$sigma^2 / 2
} else {
  # workaround because lm() with only one sequence not possible
  # make a column with intra-subject variances
  for (i in 1:nrow(data_r)) {
    data_r$varR[i] <- var(c(data_r$logpk.R1[i],
                            data_r$logpk.R2[i]), na.rm = TRUE)
  }
  dfRR    <- nrow(data_r) - 2
  s2wR    <- mean(data_r$varR)
}
CVwR      <- sqrt(exp(s2wR) - 1)
# regulatory setting for HVD(P)s
theta     <- ((log(1.25)) / 0.25)^2
# linearized scABE criterion component y
y         <- -theta * s2wR
boundy    <- y * dfRR / qchisq(1 - alpha, dfRR)
swR       <- sqrt(s2wR)
# linearized scABE criterion, 95% upper bound
critbound <- (x +y ) + sqrt(((boundx - x)^2) + ((boundy - y)^2))
# ABE, subjects random: evaluation by lmer() of package lme4
# Satterthwaite df by package lmerTest
fac       <- c("subj", "per", "seq", "treat")
data[fac] <- lapply(data[fac], factor)
data$treat <- relevel(data$treat, "T")
m3        <- lmer(log(pk) ~ seq + per + treat + (1 | subj),
                            na.action = na.omit, data = data)
sumry     <- summary(m3, ddf = "Satterthwaite")
# direct access of "treatT" does not always work (sometimes "treat1")
r         <- 1 + length(unique(data$seq)) - 1 +
             length(unique(data$per)) - 1 +
             length(unique(data$treat)) - 1
pe2       <- exp(sumry$coef[r, "Estimate"])
ci2       <- exp(confint(m3, level = 1 - 2 * alpha, method = "Wald")[r + 2, ])
CVw       <- sqrt(exp(m3@sigma^2) - 1)
options(oc) # restore options
RSABE     <- ABE <- "fail" # the nulls
if (ds == "per3full") {
  note <- "\nNote: Based on workaround; questionable!\n"
} else {
  note <- "\n"
}
# assess RSABE or ABE dependent on swR
if (swR >= 0.294) { # assess RSABE (critbound <=0 and pe within 0.8–1.25)
  if (critbound <= 0 & (pe1 >= 0.8 & pe1 <= 1.25)) RSABE <- "pass"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9), ", RSABE acceptable)")
  } else {
    txt1 <- "(< 0.294, apply ABE instead)"
  }
  if (pe1 >= 0.8 & pe1 <= 1.25) txt2 <- "(pass)" else txt2 <- "(fail)"
  if (critbound <= 0) {
    txt3 <- paste0("(", substr("\\<= 0", 2, 9), ", pass)")
  } else {
    txt3 <- "(> 0, fail)"
  }
  cat("\nswR       =", sprintf("%.5g", swR), txt1,
      "\nCVwR      =", sprintf("%6.2f%%", 100 * CVwR),
      "\nPE        =", sprintf("%6.2f%%", 100 * pe1), txt2,
      "\n90% CI    =", sprintf("%6.2f%% \u002D %.2f%%", 100 * ci1[1], 100 * ci1[2]),
      "(informative only)",
      "\ncritbound =", sprintf("%+.5g", critbound), txt3,
      "\naggregate =", RSABE, note)
} else {            # assess ABE (ci within 0.8–1.25)
  if (round(ci2[1], 4) >= 0.8 & round(ci2[2], 4) <= 1.25) ABE <- "pass"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9), ", apply RSABE instead)")
  } else {
    txt1 <- "(< 0.294, ABE applied)"
  }
  cat("\nswR    =", sprintf("%.5g", swR), txt1,
      "\nCVw    =", sprintf("%6.2f%%", 100 * CVw),
      "\nPE     =", sprintf("%.2f%%", 100 * pe2),
      "\n90% CI =", sprintf("%.2f%% \u002D %.2f%%", 100 * ci2[1], 100 * ci2[2]),
      "\nABE    =", ABE, "\n")
}

—
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
 
mittyri
★★  

Russia,
2023-07-08 00:52
(434 d 21:30 ago)

@ Helmut
Posting: # 23663
Views: 2,823
 

 a bit corrected workaround

Post reply
Hi Helmut!

❝ if (ds == "per4full") {

❝   # now the equivalent of SAS code for R-R, dlat analysis

❝   m2 <- lm(RR ~ seq, data = data_r)

❝   dfRR <- m2$df.residual

❝   s2wR <- summary(m2)$sigma^2 / 2

❝ } else {

❝   # workaround because lm() with only one sequence not possible

❝   # make a column with intra-subject variances

❝   for (i in 1:nrow(data_r)) {

❝     data_r$varR[i] <- var(c(data_r$logpk.R1[i],

❝                             data_r$logpk.R2[i]), na.rm = TRUE)

❝   }

❝   dfRR <- nrow(data_r) - 2

❝   s2wR <- mean(data_r$varR)

❝ }


I suggest to substitute the workaround to dinosaur method:
  sumsq <- 0
  for (Seq in unique(data_r$seq)) {
    CurrentData <- data_r$RR[data_r$seq == Seq]
    sumsq <- sumsq + sum((CurrentData-mean(CurrentData))^2)
  }
  NSeqs <- length(unique(data_r$seq))
  s2wR    <- sumsq/2/(nrow(data_r) - NSeqs) 
  dfRR    <- nrow(data_r) - NSeqs

at least it gives exactly the same results as classical one.
Yes, ABE is still under fire, the best approach you've found is here

Sorry, right now don't have a clue regarding differences with Winnonlin in RSABE part
—
Kind regards,
Mittyri
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2023-07-08 14:15
(434 d 08:07 ago)

@ mittyri
Posting: # 23664
Views: 2,718
 

 THX! Another attempt…

Post reply
Hi mittyri,

❝ I suggest to substitute the workaround to dinosaur method:

❝ […]

❝ at least it gives exactly the same results as classical one.


THX, that’s what’s given in the formula of the guidance. Avoing its overhead, more than 10times faster than lm(). :-D

❝ Yes, ABE is still under fire, the best approach you've found is here


Of course, this model is not exactly the one of the FDA. Seemingly (always?) a bit more conservative. From a regulatory perspective this should be acceptable. But what if the FDA’s passes and the other fails? Sponsors will not like that.

New [image]-script at the end. Allows also testing with local files. Below always RSABE followed by ABE. Results are similar but not the same.
  • The 4-period full replicate
    swR           :   0.45168 (>= 0.294: RSABE acceptable)
    CVwR          :  47.57%
    CVwT          :  35.49%   (informative only)
    implied limits:  66.82% - 149.65%
    PE            : 115.46%   (within 80% - 125%: pass)
    90% CI        : 106.39% - 125.31% (informative only)
    critbound     : -0.09454  (nonpositive: pass)
    aggregate     : pass
    Phoenix/WinNonlin
    swR           :   0.44645
    CVwR          :  46.96%
    PE            : 115.46%
    90% CI        : 106.39% - 125.31%
    critbound     : -0.092076

    swR    :   0.45168 (>= 0.294: should apply RSABE instead)
    CVwR   :  48.09%   (informative only)
    CVwT   :  46.72%   (informative only)
    PE     : 115.80%
    90% CI : 107.13% - 125.17%
    ABE    : fail      (one CL outside 80.00% - 125.00%)
    Phoenix/WinNonlin
    CVwR   :  48.09%
    CVwT   :  35.29%
    PE     : 115.66%
    90% CI : 107.10% - 124.89%     (pass!)

  • The 3-period full replicate
    swR           :   0.54127 (>= 0.294: RSABE acceptable)
    CVwR          :  58.34%
    CVwT          :  30.67%   (informative only)
    implied limits:  61.69% - 162.11%
    PE            : 124.52%   (within 80% - 125%: pass)
    90% CI        : 113.72% - 136.34% (informative only)
    critbound     :  -0.1022  (nonpositive: pass)
    aggregate     : pass
    Phoenix/WinNonlin not possible

    swR    :   0.54127 (>= 0.294: should apply RSABE instead)
    CVwR   :  49.63%   (informative only)
    CVwT   :  46.67%   (informative only)
    PE     : 124.38%
    90% CI : 113.10% - 136.78%
    ABE    : fail      (one CL outside 80.00% - 125.00%)
    Phoenix/WinNonlin
    swR    :   0.46667
    CVwR   :  49.33%
    CVwT   :  31.50%
    PE     : 124.28%
    90% CI : 113.17% - 136.49%

❝ Sorry, right now don't have a clue regarding differences with Winnonlin in RSABE part


So do I… We cross-validated it, right?

It’s funny that the FDA’s draft guidance of December 2022 states:

For example, a three-period design [TRT|RTR], could be used. A fully replicated design can estimate the subject-by-formulation interaction variance components.

True, but not with the FDA’s code.

PS @Achievwin: You keep us busy. ;-)

get.data <- function(file) { # import data and attach libraries
  require(nlme)    # for the mixed-effects model of ABE
  require(emmeans) # for its PE and CI
  # must be a flat file with column separator ',' and decimal separator '.'
  # arbitrary column names with the order below; untransformed PK-response
  data        <- read.csv(file = file,
                          stringsAsFactors = FALSE) # for R <4.0.0
  names(data) <- c("subj", "per", "seq", "treat", "pk")
  nseq        <- length(unique(data$seq))
  nper        <- length(unique(data$per))
  if (nper == 2 & nseq == 2) stop ("Simple crossover not supported!")
  if (nper == 3 & nseq == 3) stop ("Partial replicate design not supported!")
  if (nper < 3 & !nseq == 2) stop ("Design not supported!")
  if (nper == 4) ds <- "repl4"
  if (nper == 3) ds <- "repl3"
  res <- list(data = data, ds = ds)
  return(res)
}
###########################################################
# The FDA’s RSABE and ABE in full replicate designs with  #
# two sequences (a desperate attempt)                     #
# ------------------------------------------------------- #
# Detlew’s original code                                  #
# https://forum.bebac.at/forum_entry.php?id=15490         #
# Michael’s direct implementation of dlat                 #
# https://forum.bebac.at/forum_entry.php?id=23663         #
# ABE with the mixed-effects model with unequal variances #
# https://forum.bebac.at/forum_entry.php?id=20503         #
###########################################################
# select one of the examples or a local file
file      <- "https://bebac.at/downloads/ds01.csv" # EMA data set I (TRTR|RTRT)
file      <- "https://bebac.at/downloads/ds03.csv" # period 3 of I removed (TRT|RTR)
res       <- get.data(file)
data      <- res$data
ds        <- res$ds
# sort by subject, treatment and period to get replicate number
# this works in contrast to the SAS code also for
# other sequences than TRTR|RTRT
data      <- data[order(data$subj, data$treat, data$per), ]
data$repl <- 1
for (i in c(2:nrow(data))) {
  if (data$subj[i] != data$subj[i-1] | data$treat[i] != data$treat[i-1]) {
    data$repl[i] <- 1
  } else {
    data$repl[i] <- data$repl[i-1] + 1
  }
}
# make a column with treatment and repl, i.e T1,T2, R1,R2
data$code2 <- paste0(data$treat, data$repl)
data       <- data[order(data$subj, data$per), ]
# reshape to wide with colums subj, seq, pk.R1, pk.R2, pk.T1, pk.T2
# and log-transform pk
data_r     <- reshape(data = data, direction = "wide", v.names = "pk",
                      idvar = c("subj", "seq"), timevar = "code2",
                      drop = c("per", "treat", "repl"))
data_r[, 3:6]      <- log(data_r[, 3:6])
names(data_r)[3:6] <- sub("pk.", "logpk.", names(data_r))[3:6]
# T-R analysis, ilat in the SAS code
# modified for TRT|RTR and TTR|RRT
data_r$TR <- NA
if (ds == "repl4") {          # 4-period full replicate
  data_r$TR <- (data_r$logpk.T1 + data_r$logpk.T2 -
                data_r$logpk.R1 - data_r$logpk.R2) / 2
} else {                      # 3-period full replicate
  for (i in 1:nrow(data_r)) {
    if (data_r$seq[i] == "RTR" | data_r$seq[i] == "RRT")
      data_r$TR[i] <- data_r$logpk.T1[i] -
                      (data_r$logpk.R1[i] + data_r$logpk.R2[i]) / 2
    if (data_r$seq[i] == "TRT" | data_r$seq[i] == "TTR")
      data_r$TR[i] <- (data_r$logpk.T1[i] + data_r$logpk.T2[i]) / 2 -
                       data_r$logpk.R1[i]   }
}
# get rid of subjects with missing periods (TR == NA)
data_r    <- data_r[complete.cases(data_r$TR), ]
# ilat analysis, ANOVA with seq as effect
data_r$seq <- as.factor(data_r$seq)
# with standard contr.treatment we get not the wanted intercept!
# with that the intercept is intercept+seq1
oc        <- options("contrasts") # save options
options(contrasts = c("contr.sum", "contr.poly"))
m1        <- lm(TR ~ seq, data = data_r)
# intercept is estimate of µt-µR
alpha     <- 0.05 # level of the tests
est       <- coef(m1)[[1]]
pe1       <- exp(est)
ci1       <- as.numeric(confint(m1, "(Intercept)", level = 1 - 2 * alpha))
dfTR      <- m1$df.residual
# stderr, "the unknown x" for unbiased estimate of (µT-µR)^2
stderr    <- summary(m1)$coefficients["(Intercept)", "Std. Error"]
# linearized scABE criterion component x
x         <- est^2 - stderr^2
boundx    <- max(abs(ci1))^2
ci1       <- exp(ci1)
data_r$RR <- data_r$logpk.R1 - data_r$logpk.R2
data_r$TT <- data_r$logpk.T1 - data_r$logpk.T2
# get rid of subjects where R or T is not repeated (RR == NA or TT = NA)
data_rr   <- data_r[complete.cases(data_r$RR), c(1:2, 8)]
data_tt   <- data_r[complete.cases(data_r$TT), c(1:2, 9)]
# dlat analysis acc. to the formula in the guidance; THX to Michael!
# lm() with only one sequence in 3-period full replicates not possible
# gives the same result as lm() for 4-period full replicate
# more than 10times faster than lm(RR ~ seq, data = data_rr)
seqs      <- unique(as.character(data_rr$seq))
n         <- nrow(data_rr)
m         <- length(seqs)
denom     <- 2 * (n - m)
numer     <- 0
for (i in seqs) {
  D       <- data_rr$RR[data_rr$seq == i]
  numer   <- numer + sum((D - mean(D, na.rm = TRUE))^2)
}
dfRR      <- n - m
s2wR      <- numer / denom
# the same for T-T (informative!)
seqs      <- unique(as.character(data_tt$seq))
n         <- nrow(data_tt)
m         <- length(seqs)
denom     <- 2 * (n - m)
numer     <- 0
for (i in seqs) {
  D       <- data_tt$TT[data_tt$seq == i]
  numer   <- numer + sum((D - mean(D, na.rm = TRUE))^2)
}
dfTT      <- n - m
s2wT      <- numer / denom
CVwR1     <- sqrt(exp(s2wR) - 1)
CVwT1     <- sqrt(exp(s2wT) - 1)
# regulatory setting for HVD(P)s
theta     <- ((log(1.25)) / 0.25)^2
# linearized scABE criterion component y
y         <- -theta * s2wR
boundy    <- y * dfRR / qchisq(1 - alpha, dfRR)
swR       <- sqrt(s2wR)
# linearized scABE criterion, 95% upper bound
critbound <- (x +y ) + sqrt(((boundx - x)^2) + ((boundy - y)^2))
# ABE: evaluation by lme() of 'nlme' and confint() of 'emmeans'
# Satterthwaite df
fac       <- c("subj", "per", "seq", "treat")
data[fac] <- lapply(data[fac], factor)
m2        <- lme(log(pk) ~ seq + per + treat,
                           random = ~ treat | subj,
                           data = data, weights = varIdent(~ treat),
                           method = "REML", na.action = na.exclude,
                           control = list(opt = "optim"))
res       <- confint(emmeans(m2, pairwise ~ treat, mode = "satterth"),
                     level = 1 - 2 * alpha)
# note: order reversed because none of relevel(), reverse = TRUE,
#       or ref = "R" works as desired (gives always contrast R - T)
pe2       <- exp(-res$contrasts$estimate)
ci2       <- exp(-c(res$contrasts$upper.CL, res$contrasts$lower.CL))
CVwR2     <- sqrt(exp(res$emmeans$SE[1] * 2) - 1)
CVwT2     <- sqrt(exp(res$emmeans$SE[2] * 2) - 1)
options(oc) # restore options
# assess RSABE or ABE dependent on swR
if (swR >= 0.294) { # assess RSABE (critbound <=0 and pe within 0.8–1.25)
  IL      <- setNames(exp(c(-1, +1) * log(1.25) / 0.25 *swR),
                      c("lower", "upper"))
  RSABE   <- "fail" # the null
  if (critbound <= 0 & (pe1 >= 0.8 & pe1 <= 1.25)) RSABE <- "pass"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9), ": RSABE acceptable)")
  } else {
    txt1 <- "(< 0.294, apply ABE instead)"
  }
  if (pe1 >= 0.8 & pe1 <= 1.25) {
    txt2 <- "  (within 80% \u002D 125%: pass)"
  } else {
    txt2 <- "  (outside 80% \u002D 125%: fail)"
  }
  if (critbound <= 0) {
    txt3 <- " (nonpositive: pass)"
  } else {
    txt3 <- " (positive: fail)"
  }
  cat("\nswR           :", sprintf("%9.5g", swR), txt1,
      "\nCVwR          :", sprintf("%6.2f%%", 100 * CVwR1),
      "\nCVwT          :", sprintf("%6.2f%%", 100 * CVwT1),
      "  (informative only)",
      "\nimplied limits:", sprintf("%6.2f%% \u002D %.2f%%",
                                   100 * IL[["lower"]], 100 * IL[["upper"]]),
      "\nPE            :", sprintf("%6.2f%%", 100 * pe1), txt2,
      "\n90% CI        :", sprintf("%6.2f%% \u002D %.2f%%",
                                   100 * ci1[1], 100 * ci1[2]),
      "(informative only)",
      "\ncritbound     :", sprintf("%+8.5g", critbound), txt3,
      "\naggregate     :", RSABE, "\n")
} else {            # assess ABE (ci within 0.8–1.25)
  fails <- 0
  if (round(ci2[1], 4) < 0.8) fails <- fails + 1
  if (round(ci2[2], 4) > 1.25) fails <- fails + 1
  if (fails == 1) ABE <- "fail      (one CL outside 80.00% \u002D 125.00%)"
  if (fails == 2) ABE <- "fail      (CI not within 80.00% \u002D 125.00%)"
  if (round(ci2[1], 4) > 1.25 | round(ci2[1], 4) < 0.8)
    ABE <- "inequivalent (CI entirely outside 80.00% \u002D 125.00%)"
  if (round(ci2[1], 4) >= 0.8 & round(ci2[2], 4) <= 1.25)
    ABE <- "pass     (CI within 80.00% \u002D 125.00%)"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9),
                   ": should apply RSABE instead)")
  } else {
    txt1 <- "(< 0.294, ABE applied)"
  }
  cat("\nswR    :", sprintf("%9.5g", swR), txt1,
      "\nCVwR   :", sprintf("%6.2f%%", 100 * CVwR2), "  (informative only)",
      "\nCVwT   :", sprintf("%6.2f%%", 100 * CVwT2), "  (informative only)",
      "\nPE     :", sprintf("%6.2f%%", 100 * pe2),
      "\n90% CI :", sprintf("%6.2f%% \u002D %.2f%%",
                            100 * ci2[1], 100 * ci2[2]),
      "\nABE    :", ABE, "\n")
}

—
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
 
mittyri
★★  

Russia,
2023-07-09 23:45
(432 d 22:37 ago)

@ Helmut
Posting: # 23665
Views: 2,696
 

 complete.cases for dlat and ilat?

Post reply
Hi Helmut!

Got a flashback. Take a look:

 data_r    <- data_r[complete.cases(data_r$TR), ]
 <...>
 data_r$RR <- data_r$logpk.R1 - data_r$logpk.R2
 data_r$TT <- data_r$logpk.T1 - data_r$logpk.T2
 data_rr   <- data_r[complete.cases(data_r$RR), c(1:2, 8)]  
 data_tt   <- data_r[complete.cases(data_r$TT), c(1:2, 9)]

data_r should not be used for dlat AND ilat (if the goal is to mimic Phoenix template where dlat is calculated for all data where RR is available, not only full periods data)
—
Kind regards,
Mittyri
 
Achievwin
★★  

US,
2023-07-10 22:43
(431 d 23:39 ago)

@ mittyri
Posting: # 23669
Views: 2,630
 

 a bit corrected workaround

Post reply
PS @Achievwin: You keep us busy. ;-)

Thats why you called Expert Gurus ;-):cool::-D
 
UA Flag
Activity
 Admin contact
23,224 posts in 4,878 threads, 1,652 registered users;
38 visitors (1 registered, 37 guests [including 9 identified bots]).
Forum time: 22:22 CEST (Europe/Vienna)

We absolutely must leave room for doubt
or there is no progress and no learning.
There is no learning without having to pose a question.
And a question requires doubt.
People search for certainty.
But there is no certainty.    Richard Feynman

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