loosing specificity due to low sensitivity [Regulatives / Guidelines]
Hi mittyri,
THX!
I hope so. It is on my to-do-list since last summer…
Right. I have no idea where this idea come from.
I used code developed by Martin and me many years ago to simulate 2×2 designs followed by a simple
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
❝ 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 🖖🏼 Довге життя Україна!![[image]](https://static.bebac.at/pics/Blue_and_yellow_ribbon_UA.png)
Helmut Schütz
![[image]](https://static.bebac.at/img/CC by.png)
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
![[image]](https://static.bebac.at/pics/Blue_and_yellow_ribbon_UA.png)
Helmut Schütz
![[image]](https://static.bebac.at/img/CC by.png)
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Complete thread:
- Russian «Экспертами» and their hobby Helmut 2017-04-29 00:46 [Regulatives / Guidelines]
- Low power of Group-by-Treatment interaction mittyri 2017-04-29 22:57
- Let’s forget the Group-by-Treatment interaction, please! Helmut 2017-04-30 13:54
- Let’s forget the Group-by-Treatment interaction, please! ElMaestro 2017-05-01 16:19
- Some answers Helmut 2017-05-02 01:10
- Some answers ElMaestro 2017-05-02 09:04
- Example Helmut 2017-05-02 12:35
- Sensitivity of term? mittyri 2017-05-02 18:29
- Simulations Helmut 2017-05-05 14:38
- loosing specificity due to low sensitivity mittyri 2017-05-08 23:28
- loosing specificity due to low sensitivityHelmut 2017-05-09 00:55
- loosing specificity due to low sensitivity mittyri 2017-05-08 23:28
- Loss in power Helmut 2017-05-06 17:31
- Interval between groups Helmut 2017-05-08 19:02
- IMP handling mittyri 2017-05-08 23:40
- IMP handling Helmut 2017-05-09 01:08
- IMP handling mittyri 2017-05-08 23:40
- Loss in power Helmut 2017-05-14 17:22
- Simulations Helmut 2017-05-05 14:38
- Some answers ElMaestro 2017-05-02 09:04
- No convergence in JMP and Phoenix WinNonlin Helmut 2017-05-25 15:26
- Ouch?!??? ElMaestro 2017-05-25 16:24
- Some answers Helmut 2017-05-02 01:10
- Let’s forget the Group-by-Treatment interaction, please! ElMaestro 2017-05-01 16:19
- Let’s forget the Group-by-Treatment interaction, please! Helmut 2017-04-30 13:54
- Russian «Экспертами» and their hobby Artem Gusev 2017-05-02 16:13
- be careful with mixed models mittyri 2017-05-02 17:53
- be careful with mixed models Artem Gusev 2017-05-03 11:02
- p-value(s) in model 2 Helmut 2017-05-05 14:48
- be careful with mixed models mittyri 2017-05-02 17:53
- Russian «Экспертами» following the EEU GLs Helmut 2017-05-24 20:17
- Russian «Экспертами» following the EEU GLs Beholder 2017-05-24 22:37
- Penalty for carelessness mittyri 2017-05-25 08:52
- Russian «Экспертами» following the EEU GLs Beholder 2017-05-25 10:43
- Russian «Экспертами» following the EEU GLs Mikalai 2018-01-04 10:43
- Belarus = member of the EEU Helmut 2018-01-04 13:08
- Belarus = member of the EEU Mikalai 2018-01-04 19:49
- Trying your model for EEU mittyri 2018-01-04 22:04
- Trying your model for EEU Helmut 2018-01-05 00:06
- help us to stop it, please... Astea 2018-01-10 12:09
- help us to stop it, please... Beholder 2018-01-10 12:49
- regulators convinced by science? d_labes 2018-01-10 15:15
- regulators convinced by science? Beholder 2018-01-10 17:14
- Чёрт побери! d_labes 2018-01-10 18:53
- regulators convinced by science? Astea 2018-01-10 19:10
- regulators convinced by science? Beholder 2018-01-10 17:14
- help us to stop it, please... Astea 2018-01-10 12:09
- Trying your model for EEU Helmut 2018-01-05 00:06
- Belarus = member of the EEU Helmut 2018-01-04 13:08
- Russian «Экспертами» following the EEU GLs Beholder 2017-05-24 22:37
- Low power of Group-by-Treatment interaction mittyri 2017-04-29 22:57