loosing specificity due to low sensitivity [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2017-05-09 02:55 (2911 d 18:57 ago) – Posting: # 17329
Views: 31,578

Hi mittyri,

❝ you've made a great work!


THX!

❝ Won't it be published?


I hope so. It is on my to-do-list since last summer…

❝ In your examples (simulations/practice) you showed the TxG test is not a good idea.

❝ I was impressed by this: […]

❝

❝ I see that the sensitivity is really low, but I think it is not a good idea to compensate it with low specificity (high false positive).


Right. I have no idea where this idea come from. :confused:

❝ Wouldn't you mind to publish the code of data building for simulations?


I used code developed by Martin and me many years ago to simulate 2×2 designs followed by a simple rbind(part1, part2). I improved it to speed things up (which worked). Unfortunately I screwed up a couple of days ago (no version control, saving over). Shit.
In the meantime you can (mis)use the code Detlew distributed last year to simulate replicate designs. I duplicated his functions. Start with the RTRT|TRTR and set CVWT = CVwR = CVbT = CVbR. The code below shows the relevant changes after his example and how to sim. For the plot you have to attach package lattice.
set.seed(123456)
G1 <- 0.95 # T/R in group 1
G2 <- 1/G1 # T/R in group 2
mvc1 <- mean_vcov(c("TRTR", "RTRT"), muR=log(100), ldiff=log(G1),
                 sWT=CV2se(0.3), sWR=CV2se(0.3),
                 sBT=CV2se(0.3), sBR=CV2se(0.3), rho=1)
mvc2 <- mean_vcov(c("TRTR", "RTRT"), muR=log(100), ldiff=log(G2),
                 sWT=CV2se(0.3), sWR=CV2se(0.3),
                 sBT=CV2se(0.3), sBR=CV2se(0.3), rho=1)
# get the data
ow    <- options()
nsims <- 1e4
sig   <- pass.2 <- pass.3 <- pass.3.1 <- pass.3.2 <- pass.3.a <- 0
PE2   <- PE3.1 <- PE3.2 <- PE3 <- p.GxT <- numeric(0)
alpha <- 0.05
p.level <- 0.1
L <- 80
U <- 125
sub.seq <- 8 # subjects / sequence (within each group)
for (j in 1:nsims) {
  part1 <- prep_data(seqs=c("TRTR", "RTRT"), rep(sub.seq, 2),
                     metric="PK", dec=3, mvc_list=mvc1)
  part1 <- part1[, !(names(part1) %in% c("seqno", "logval"))]
  part1 <- part1[!part1$period %in% c(3, 4), ]
  part1$sequence[part1$sequence == "TRTR"] <- "TR"
  part1$sequence[part1$sequence == "RTRT"] <- "RT"
  part1$group <- 1
  part2 <- prep_data(seqs=c("TRTR", "RTRT"), rep(sub.seq, 2),
                     metric="PK", dec=3, mvc_list=mvc2)
  part2 <- part2[, !(names(part2) %in% c("seqno", "logval"))]
  part2 <- part2[!part2$period %in% c(3, 4), ]
  part2$sequence[part2$sequence == "TRTR"] <- "TR"
  part2$sequence[part2$sequence == "RTRT"] <- "RT"
  part2$subject <- part2$subject+sub.seq*2
  part2$group <- 2
  study <- rbind(part1, part2)
  study$subject   <- factor(study$subject)
  study$period    <- factor(study$period)
  study$sequence  <- factor(study$sequence)
  study$treatment <- factor(study$treatment)
  study$group     <- factor(study$group)
  options(contrasts=c("contr.treatment", "contr.poly"), digits=12)
  # model 1 of pooled data
  model1 <- lm(log(PK)~group+sequence+treatment+group*sequence+
                       group*period+group*treatment+subject%in%group*sequence,
                       data=study)
  p.GxT[j] <- anova(model1)[["group:treatment", "Pr(>F)"]]

  if (p.GxT[j] >= p.level) { # if no sign. interaction: model 2 of pooled data
    model2 <- lm(log(PK)~group+sequence+treatment+group*sequence+
                         group*period+subject%in%group*sequence,
                         data=study)
    CI2 <- round(100*exp(confint(model2, "treatmentT", level=1-2*alpha)), 2)
    if (CI2[1] >= L & CI2[2] <= U) {
      pass.2 <- pass.2 + 1 # count passing studies
      PE2[pass.2] <- as.numeric(exp(coef(model2)["treatmentT"]))
    }
  } else { # sign. interaction (otherwise): model 3 of both groups
    sig <- sig + 1 # count studies with significant interaction
    # first group (use part1 data)
    model3.1 <- lm(log(PK)~sequence+treatment+period+subject%in%sequence,
                           data=part1)
    CI3.1 <- round(100*exp(confint(model3.1, "treatmentT", level=1-2*alpha)), 2)
    if (CI3.1[1] >= L & CI3.1[2] <= U) {
      pass.3.1 <- pass.3.1 + 1 # count passing studies
      PE3.1[pass.3.1] <- exp(coef(model3.1)["treatmentT"])
    }
    # second group (use part2 data)
    model3.2 <- lm(log(PK)~sequence+treatment+period+subject%in%sequence,
                           data=part2)
    CI3.2 <- round(100*exp(confint(model3.2, "treatmentT", level=1-2*alpha)), 2)
    if (CI3.2[1] >= L & CI3.2[2] <= U) {
      pass.3.2 <- pass.3.2 + 1 # count passing studies
      PE3.2[pass.3.2] <- as.numeric(exp(coef(model3.2)["treatmentT"]))
    }
    # check whether /both/ groups pass (haha)
    if ((CI3.1[1] >= L & CI3.1[2] <= U) &
        (CI3.2[1] >= L & CI3.2[2] <= U)) pass.3.a <- pass.3.a + 1
  }
  # model 3 of pooled data (simple 2x2x2 crossover)
  model3 <- lm(log(PK)~sequence+treatment+period+subject%in%sequence,
                       data=study)
  CI3 <- round(100*exp(confint(model3, "treatmentT", level=1-2*alpha)), 2)
  if (CI3[1] >= L & CI3[2] <= U) {
    pass.3 <- pass.3 + 1 # count passing studies
    PE3[pass.3] <- as.numeric(exp(coef(model3)["treatmentT"]))
  }
} # end of sim loop
PE2est <- prod(PE2)^(1/length(PE2)) # geom. mean of PEs (passing with model 2)
if (length(PE3.1) > 0) { # geom. mean of PEs (passing with model 3; group 1)
  PE3.1est <- prod(PE3.1)^(1/length(PE3.1))
} else {
  PE3.1est <- NA
}
if (length(PE3.2) > 0) { # geom. mean of PEs (passing with model 3; group 2)
  PE3.2est <- prod(PE3.2)^(1/length(PE3.1))
} else {
  PE3.2est <- NA
}
PE3est <- prod(PE3)^(1/length(PE3)) # geom. mean of PEs (passing with model 3)
options(ow) # restore options
x <- c(0.25, 0.75)
y <- as.numeric(quantile(p.GxT, probs=x))
b <- diff(y)/diff(x)
a <- y[2]-b*x[2]
numsig <- length(which(p.GxT < p.level))
MajorInterval <- 5 # interval for major ticks
MinorInterval <- 4 # interval within major
Major <- seq(0, 1, 1/MajorInterval)
Minor <- seq(0, 1, 1/(MajorInterval*MinorInterval))
labl  <- sprintf("%.1f", Major)
ks    <- ks.test(x=p.GxT, y="punif", 0, 1)
if (G1 != G2) {
  main <- list(label=paste0("\nSimulation of \'true\' interaction\nT/R (G1) ",
                            G1, ", T/R (G2) ", round(G2, 4)), cex=0.9)
} else {
  main <- list(label=paste0("\nSimulation of no interaction\nT/R (G1, G2) ",
                            G1), cex=0.9)
}
if (ks$p.value == 0) {
  sub <- list(label=sprintf("Kolmogorov-Smirnov test: p <%1.5g",
              .Machine$double.eps), cex=0.8)
} else {
  sub <- list(label=sprintf("Kolmogorov-Smirnov test: p %1.5g",
              ks$p.value), cex=0.8)
}
trellis.par.set(layout.widths=list(right.padding=5))
qqmath(p.GxT, distribution=qunif,
  prepanel=NULL,
  panel=function(x) {
    panel.grid(h=-1, v=-1, lty=3)
    panel.abline(h=p.level, lty=2)
    panel.abline(c(0, 1), col="lightgray")
    panel.abline(a=a, b=b)
    panel.qqmath(x, distribution=qunif, col="blue", pch=46) },
    scales=list(x=list(at=Major), y=list(at=Major),
                tck=c(1, 0), labels=labl, cex=0.9),
    xlab="uniform [0, 1] quantiles",
    ylab="p (Group-by-Treatment Interaction)",
    main=main, sub=sub, min=0, max=1)
trellis.focus("panel", 1, 1, clip.off=TRUE)
panel.axis("bottom", check.overlap=TRUE, outside=TRUE, labels=FALSE,
           tck=0.5, at=Minor)
panel.axis("left", check.overlap=TRUE, outside=TRUE, labels=FALSE,
           tck=0.5, at=Minor)
panel.polygon(c(0, 0, numsig/nsims, numsig/nsims, 0),
              c(0, rep(p.level, 2), 0, 0), lwd=1, border="red")
trellis.unfocus()
cat("Model 1: p(G\u00D7T) <0.1 in", sprintf("%.2f%%", 100*sig/nsims),
"of studies.",
"\nEvaluation of studies with p(G\u00D7T) <0.1 (Groups):",
"\n  passed model 3 (1)      :",
  sprintf("%5.2f%%", 100*pass.3.1/sig), "(of tested); PE",
  sprintf("%6.2f%%", 100*PE3.1est),
"\n                                  range of PEs:",
  sprintf("%5.2f%% to %6.2f%%", 100*range(PE3.1)[1], 100*range(PE3.1)[2]),
"\n  passed model 3 (2)      :",
  sprintf("%5.2f%%", 100*pass.3.2/sig), "(of tested); PE",
  sprintf("%6.2f%%", 100*PE3.2est),
"\n                                  range of PEs:",
  sprintf("%5.2f%% to %6.2f%%", 100*range(PE3.2)[1], 100*range(PE3.2)[2]),
"\n  passed model 3 (1 and 2):",
  sprintf("%5.2f%%", 100*pass.3.a/sig), "(of tested)",
"\nEvaluation of studies with p(G\u00D7T) \u22650.1 (pooled):",
"\n  passed model 2          :",
  sprintf("%5.2f%%", 100*pass.2/nsims), "(overall)",
"\n                           ",
  sprintf("%5.2f%%", 100*pass.2/(nsims-sig)), "(of tested); PE",
  sprintf("%6.2f%%", 100*PE2est),
"\n                                  range of PEs:",
  sprintf("%5.2f%% to %6.2f%%", 100*range(PE2)[1], 100*range(PE2)[2]),
"\nStudies passing any of model 2 or 3:",
  sprintf("%5.2f%%", 100*(pass.3.1/nsims+pass.3.2/nsims+pass.2/nsims)),
"\nCriteria for simple model fulfilled:",
"\n  passed model 3          :",
  sprintf("%5.2f%%;             PE", 100*pass.3/nsims),
  sprintf("%6.2f%%", 100*PE3est),
"\n                                  range of PEs:",
  sprintf("%5.2f%% to %6.2f%%", 100*range(PE3)[1], 100*range(PE3)[2]), "\n")
round(summary(p.GxT), 7)
—
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,424 posts in 4,927 threads, 1,673 registered users;
102 visitors (0 registered, 102 guests [including 6 identified bots]).
Forum time: 21:52 CEST (Europe/Vienna)

There are two possible outcomes: if the result confirms the
hypothesis, then you’ve made a measurement. If the result is
contrary to the hypothesis, then you’ve made a discovery.    Enrico Fermi

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