loosing specificity due to low sensitivity [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2017-05-09 02:55 (2395 d 04:53 ago) – Posting: # 17329
Views: 29,391

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
22,811 posts in 4,783 threads, 1,641 registered users;
18 visitors (0 registered, 18 guests [including 8 identified bots]).
Forum time: 06:49 CET (Europe/Vienna)

Every man gets a narrower and narrower field of knowledge
in which he must be an expert in order to compete with other people.
The specialist knows more and more about less and less
and finally knows everything about nothing.    Konrad Lorenz

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