I would like to produce some power curves with effect size on the x axis, and power on the y axis. I am only slightly modifying the code from this source: https://moderndata.plot.ly/power-curves-r-plotly-ggplot2/ to create a power curve. The author is using a two sided t-test, while I am using a two proportions test.
I ran the author's code for testing purposes, and it works fine. When I run my own code with the slight modifications (I only modified the part in the loop to do a two proportions test instead of a t test), nothing prints in the plot. I cannot figure out what I am doing wrong. Please note that I kept the author's column names for now, as I just want to test my code.
library(pwr) # for power calcs
library(dplyr) # for data manipulation
library(tidyr) # for data manipulation
library(ggplot2) # for plotting power curves
library(plotly) # for interactive power curves
# Generate power calculations
ptab <- cbind(NULL, NULL)
for (i in seq(0,1, length.out = 200)){
pwrt1 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.05),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt2 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.10),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt3 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.15),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt4 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.20),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt5 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.25),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt6 <- pwr.2p.test(h = ES.h(p1 = 0.3, p2 = 0.30),n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
ptab <- rbind(ptab, cbind(pwrt1$h, pwrt1$power,
pwrt2$h, pwrt2$power,
pwrt3$h, pwrt3$power,
pwrt4$h, pwrt4$power,
pwrt5$h, pwrt5$power,
pwrt6$h, pwrt6$power))
}
ptab <- cbind(seq_len(nrow(ptab)), ptab)
colnames(ptab) <- c("id","n1=28, n2=1406.effect size","n1=28,
n2=1406.power",
"n1=144, n2=1290.effect size","n1=144, n2=1290.power",
"n1=287, n2=1147.effect size","n1=287, n2=1147.power",
"n1=430, n2=1004.effect size","n1=430, n2=1004.power",
"n1=574, n2=860.effect size","n1=574, n2=860.power",
"n1=717, n2=717.effect size","n1=717, n2=717.power")
# get data into right format for ggplot2
temp <- ptab %>%
as.data.frame() %>%
gather(key = name, value = val, 2:13) %>%
separate(col = name, into = c("group", "var"), sep = "\\.") %>%
spread(key = var, value = val)
# factor group
temp$group <- factor(temp$group,
levels = c("n1=28, n2=1406", "n1=144, n2=1290",
"n1=287, n2=1147", "n1=430, n2=1004",
"n1=574, n2=860", "n1=717, n2=717"))
# plot
p <- ggplot(temp, aes(x = size, y = power, color = group))+
geom_line(size=2) +
theme_bw() +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=14),
legend.text=element_text(size=14)) +
geom_vline(xintercept = .54, linetype = 2) +
geom_hline(yintercept = 0.80, linetype = 2)
p
Edit: I figured out what I was doing wrong, I was specifying the effect size option in when I didn't have to since the loop was going to go through all the effect sizes between 0 and 1. The following modification to the loop fixed my issue:
for (i in seq(0,1, length.out = 200)){
pwrt1 <- pwr.2p.test(h = i,n=16,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt2 <- pwr.2p.test(h = i,n=32,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt3 <- pwr.2p.test(h = i,n=48,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt4 <- pwr.2p.test(h = i,n=64,power = NULL,
sig.level = 0.01,alternative=c("greater"))
pwrt5 <- pwr.2p.test(h = i,n=80,power = NULL,
sig.level = 0.01,alternative=c("greater"))
summary(ptab). You'll see that for each column, the values are just the same number repeated 200 times. If you make a scatterplot oftemp, you just have 5 points that are overplotted, plus a bunch ofNAvalues