## Likely wrong workaround? [RSABE / ABEL]

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