loosing specificity due to low sensitivity [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2017-05-09 02:55 (2622 d 15:34 ago) – Posting: # 17329
Views: 30,215

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,099 posts in 4,857 threads, 1,646 registered users;
93 visitors (0 registered, 93 guests [including 8 identified bots]).
Forum time: 18:30 CEST (Europe/Vienna)

Imagine if every Thursday your shoes exploded
if you tied them the usual way.
This happens to us all the time with computers,
and nobody thinks of complaining.    Jef Raskin

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