0

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"))
3
  • 2
    Call summary(ptab). You'll see that for each column, the values are just the same number repeated 200 times. If you make a scatterplot of temp, you just have 5 points that are overplotted, plus a bunch of NA values Commented Jul 9, 2018 at 20:58
  • @camille thanks for this tip, I just figured out what I was doing wrong, and everything is printing now. Commented Jul 10, 2018 at 17:19
  • 2
    Cool, feel free to post an answer to your own question if you figured out how to solve it Commented Jul 10, 2018 at 17:23

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.