This is the R Script referenced in the blog post LINK
# ratio of true to not-true relationships
tnott <- 0.01
# sample sizes for which PPV and power are calculated
n <- c(6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44)
alpha <- 0.05
power <- vector()
waste <- vector()
wasted <- vector()
fdr <- vector()
ppv <- vector()
wastedf <- vector()
tg <- vector()
for (i in 1:20) {
power[i] <- power.t.test(n[i], delta = 1, sd = 1, sig.level = alpha, power = NULL,
type = "two.sample", alternative = "two.sided")$power
waste[i] <- alpha*2*n[i]*(1-tnott) + (1-power[i])*2*n[i]*tnott
wasted[i] <- waste[i]*100/(2*n[i])
fdr[i] <- alpha*(1-tnott)/(power[i]*tnott + alpha*(1-tnott))
ppv[i] <- 100*(1-fdr[i])
wastedf[i] <- fdr[i]*2*n[i]
}
for (i in 2:20) {
tg[i] <- (ppv[i]-ppv[i-1])/(n[i]-n[i-1])
}
grafik <- data.frame(n,power,fdr,wastedf)
# Value used to transform the data
coeff <- 0.01
# A few constants
ppvColor <- "#69b3a2"
powerColor <- rgb(1, 0, 0, 1)
ggplot(head(grafik, 80), aes(x=n)) +
geom_bar( aes(y=ppv), stat="identity", size=.1, fill=ppvColor, color="black", alpha=.4) +
geom_line( aes(y=power / coeff), size=2, color=powerColor) +
scale_y_continuous(
# Features of the first axis
name = "Positive Predictive Value (%)",
# Add a second axis and specify its features
sec.axis = sec_axis(~.*coeff, name="Power (1-ß)")
) +
xlab("Sample Size") +
theme(
axis.title.y = element_text(color = ppvColor, size=24, vjust = 2.5),
axis.title.y.right = element_text(color = powerColor, size=24, vjust = 3.5)
) +
ggtitle("Low odds (1:100)") +
# theme_minimal() +
theme(
plot.title=element_text(face='bold', colour='black', size=26)
) +
theme(axis.title.x = element_text(face="bold", size = 20, vjust = 0.1)) +
theme(axis.text.x= element_text(face="bold", size=16)) +
theme(axis.text.y= element_text(face="bold", size=16)) +
coord_fixed(ratio = 0.5)
0 Comments
Leave A Comment