## loosing specificity due to low sensitivity [Regulatives / Guidelines]

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.

❝ 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 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes