Achievwin
★

US,
2023-07-05 20:38
(343 d 15:57 ago)

Posting: # 23655
Views: 2,430

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

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

Vienna, Austria,
2023-07-06 10:46
(343 d 01:48 ago)

@ Achievwin
Posting: # 23657
Views: 2,154

## Use the one for TRTR | RTRT

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 . 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 🖖🏼 Довге життя Україна!
Helmut Schütz

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

Russia,
2023-07-06 16:07
(342 d 20:28 ago)

@ Helmut
Posting: # 23658
Views: 2,110

## RSABE R coded by Detlew

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

Vienna, Austria,
2023-07-07 16:41
(341 d 19:54 ago)

@ mittyri
Posting: # 23660
Views: 2,040

## ‘dlat’ for 3-period full replicate?

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

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

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

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

@ Helmut
Posting: # 23661
Views: 2,091

## ‘dlat’ for 3-period full replicate?

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

Vienna, Austria,
2023-07-07 22:06
(341 d 14:29 ago)

@ mittyri
Posting: # 23662
Views: 2,056

## Likely wrong workaround?

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. .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 🖖🏼 Довге життя Україна! Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes mittyri ★★ Russia, 2023-07-08 00:52 (341 d 11:43 ago) @ Helmut Posting: # 23663 Views: 2,026 ## a bit corrected workaround 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 ★★★ Vienna, Austria, 2023-07-08 14:15 (340 d 22:20 ago) @ mittyri Posting: # 23664 Views: 1,942 ## THX! Another attempt… 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() ❝ 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 -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 🖖🏼 Довге життя Україна!
Helmut Schütz

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

Russia,
2023-07-09 23:45
(339 d 12:50 ago)

@ Helmut
Posting: # 23665
Views: 1,879

## complete.cases for dlat and ilat?

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
(338 d 13:52 ago)

@ mittyri
Posting: # 23669
Views: 1,856

## a bit corrected workaround

PS @Achievwin: You keep us busy.

Thats why you called Expert Gurus