Bioequivalence and Bioavailability Forum

Hi mittyri,

❝ The Tmax evaluation for parallel design is not implemented in Phoenix Winnonlin,

Correct. IIRC I put it years ago on Pharsight’s wish-list. Not sure.

❝ Are there some ways around?

In PHX/WNL – none. I suggest .

HL <- function(x, na.action=na.omit) { # Hodges-Lehmann estimator of x
        if (!is.vector(x) || !is.numeric(x))
          stop("'x' must be a numeric vector")
        x <- na.action(x)
        y <- outer(x, x, '+')
        median(y[row(y) >= col(y)])/2
}
T  <- c(1.5, 3.5, 3.0, 3.0, 3.0,
        2.0, 2.5, 3.0, 3.0, 3.5,
        2.5, 2.0, 2.5, 2.0, 4.5)
R  <- c(3.0, 4.0, 2.0, 2.5, 2.5,
        2.5, 3.5, 2.5, 3.0, 2.5,
        2.0, 4.0, 1.5, 3.0, 4.0,
        2.5)
n1 <- length(T)
n2 <- length(R)
N  <- n1+n2
tmax <- data.frame(t=c(T, R),
          treat=factor(c(rep("T", n1), rep("R", n2))))
wilcox.test(tmax$t ~ relevel(tmax$treat, ref="R"), data=tmax,
  alternative="two.sided", paired=FALSE, exact=TRUE,
  conf.int=TRUE, conf.level=0.90)


Throws warnings because of ties and automatically switches to the normal approximation with continuity cor­rection. Furthermore, the exact method in wilcox.test() is limited to <50 values. Note that for this data set medians are dif­ferent, whereas the Hodges-Lehmann estimates are equal. The test is based on the latter. Therefore, may­be better …

txt1 <- paste0("Medians of T and R                 : ",
          median(T), ", ", median(R))
txt1 <- paste0(txt1, "\nHodges-Lehmann estimates of T and R: ",
          HL(T), ", ", HL(R))
if (length(unique(sort(rank(tmax$t)))) == N) {
  if (length(tmax$t) < 50) {
    txt1 <- c(txt1, "\nNo ties in ranked data and <50 values; exact method used.\n")
    exact <- TRUE
  } else {
    txt1 <- c(txt1, paste("\nNo ties in ranked data, but", N,
      "values (\u226550); approximation used.\n"))
    exact <- FALSE
  }
} else {
  txt1 <- c(txt1, "\n<50 values, but ties in ranked data; normal approximation used.\n")
  exact <- FALSE
}
wt1 <- wilcox.test(tmax$t ~ relevel(tmax$treat, ref="R"), data=tmax,
  alternative="two.sided", paired=FALSE, exact=exact,
  conf.int=TRUE, conf.level=0.90)
cat(txt1); print(wt1)

… which gives:

Medians of T and R                 : 3, 2.5
Hodges-Lehmann estimates of T and R: 2.75, 2.75
<50 values, but ties in ranked data; normal approximation used.

        Wilcoxon rank sum test with continuity correction

data:  tmax$t by relevel(tmax$treat, ref = "R")
W = 122, p-value = 0.9516
alternative hypothesis: true location shift is not equal to 0
90 percent confidence interval:
 -0.4999373  0.5000596
sample estimates:
difference in location
          6.736596e-05

Consider package coin which is able to handle ties …

library(coin)
txt2 <- paste0("Hodges-Lehmann estimates of T and R: ",
          HL(T), ", ", HL(R), "\n")
wt2  <- wilcox_test(tmax$t ~ relevel(tmax$treat, ref="R"), data=tmax,
  alternative="two.sided", paired=FALSE, distribution="exact",
  conf.int=TRUE, conf.level=0.90)
cat(txt2); print(wt2)

… and gives exact results:

Hodges-Lehmann estimates of T and R: 2.75, 2.75

        Exact Wilcoxon-Mann-Whitney Test

data:  tmax$t by relevel(tmax$treat, ref = "R") (R, T)
Z = 0.080982, p-value = 0.9458
alternative hypothesis: true mu is not equal to 0
90 percent confidence interval:
 -0.5  0.5
sample estimates:
difference in location
                     0

The approximation is pretty good. However, I prefer exact tests if ever possible.

See also this lengthy thread.
BTW, it is a common mis­con­ception that the test compares medians. With both setups the shift is assessed based on the Hodges-Lehmann estimates (2.75 for T and R), not the medians (T 3, R 2.5) and dis­tri­butions of samples are assumed to be identical. I can live with it. In parametric methods for ABE (and EMA’s ABEL) we assume IIDs as well.

As you please:

x        <- 3
ci       <- as.numeric(confint(wt2)$conf.int)
pe       <- as.numeric(confint(wt2)$estimate)
bp       <- boxplot(tmax$t ~ tmax$treat, plot=FALSE)
desc     <- bp$stats
colnames(desc) <- bp$names
desc     <- rbind(
  c(min(tmax$t[tmax$treat == bp$names[1]]),
    min(tmax$t[tmax$treat == bp$names[2]])), desc)
desc    <- rbind(desc,
  c(max(tmax$t[tmax$treat == bp$names[1]]),
    max(tmax$t[tmax$treat == bp$names[2]])))
rownames(desc) <- c("minimum", "lower whisker", "1st quartile", "median",
  "3rd quartile","upper whisker", "maximum")
bxp(bp, xlab="treatment", ylab="tmax",
  xlim=c(0.5, 3.5), ylim=c(ci[1], max(tmax$t)),
  main="Boxplot of R and T; HL estimate and ~90% CI of T \u2013 R",
  cex.main=1, las=1, boxwex=0.5, boxfill="bisque",
  outline=TRUE, pars=list(outpch=8, outcol="red", outcex=0.8))
mtext("T \u2013 R", 1, line=1, at=3)
lines(1+c(-1, 1)*0.0625, rep(desc["minimum", "R"], 2), col="blue")
lines(2+c(-1, 1)*0.0625, rep(desc["minimum", "T"], 2), col="blue")
lines(1+c(-1, 1)*0.0625, rep(desc["maximum", "R"], 2), col="blue")
lines(2+c(-1, 1)*0.0625, rep(desc["maximum", "T"], 2), col="blue")
axis(4, las=1)
abline(h=0, lty=3)
abline(v=2.5)
points(1:2, c(HL(T), HL(R)), col="blue", pch=18, cex=1.5)
points(x, pe, col="red", pch=18, cex=1.5)
rect(x+0.15, pe-0.1, x+0.55, pe+0.1, col="white", lty=0)
text(x+0.125, pe, labels=sprintf("%+.3f", pe), cex=0.8, adj=c(-0.2, 0.3))
lines(rep(x, 2), ci, col="red", lty=2)
lines(x+c(-1, 1)*0.125, rep(ci[1], 2), col="red")
lines(x+c(-1, 1)*0.125, rep(ci[2], 2), col="red")
text(rep(x+0.125, 2), ci, labels=sprintf("%+.3f", ci), cex=0.8, adj=c(-0.2, 0.3))
leg.txt <- c("Hodges-Lehmann estimate", "Median", "Minimum, maximum")
leg.lty <- c(0, 1, 1)
leg.lwd <- c(0, 3, 1)
leg.col <- c("blue", "black", "blue")
leg.pch <- c(18, NA, NA)
leg.pt.cex <- c(1.5, NA, NA)
if (length(bp$out) >= 1) {
  leg.txt[4] <- "Outlier"
  leg.lty[4] <- 0
  leg.lwd[4] <- 0
  leg.col[4] <- "red"
  leg.pch[4] <- 8
  leg.pt.cex[4] <- 0.8
}
legend("bottomleft", legend=leg.txt, lty=leg.lty, lwd=leg.lwd,
  col=leg.col, pch=leg.pch, pt.cex=leg.pt.cex, bty="o", bg="white",
  box.lty=0, inset=c(0.02, 0))
legend("topright", legend=c("HL estimate", "~90% CI"),
  lty=c(0, 1), lwd=c(0, 1), col=c(rep("red", 2)),
  pch=c(18, NA), pt.cex=c(1.5, NA), bty="n")
print(desc)

                 R    T
minimum       1.50 1.50
lower whisker 1.50 1.50
1st quartile  2.50 2.25
median        2.50 3.00
3rd quartile  3.25 3.00
upper whisker 4.00 3.50
maximum       4.00 4.50

---

Don’t get confused by terminology.

• Independent samples
  • Wilcoxon Rank Sum Test ≡
    
  • Wilcoxon-Mann-Whitney Test ≡
    
  • Mann–Whitney–Wilcoxon Test ≡
    
  • Mann-Whitney-U-Test.
• Dependent samples
  • Wilcoxon Signed-Rank Test ≡
    
  • Wilcoxon T Test.