Tmax evaluation for parallel groups in R [Software]

posted by Helmut Homepage – Vienna, Austria, 2015-10-03 17:26 (3424 d 02:54 ago) – Posting: # 15519
Views: 5,802

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 [image].

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)


[image]

                 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.

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

Complete thread:

UA Flag
Activity
 Admin contact
23,381 posts in 4,914 threads, 1,663 registered users;
51 visitors (0 registered, 51 guests [including 9 identified bots]).
Forum time: 19:20 CET (Europe/Vienna)

Science is built up with facts, as a house is with stones.
But a collection of facts is no more a science
than a heap of stones is a house.    Henri Poincaré

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