CODEWITHSUNDEEP
rggplot23d
user3022875
user3022875

Reputation: 9018

3d surface plot with function call

I'd like to do a surface plot where z is a call to the pwr.t.test function evaluated at x an y

d = rep( seq(.05,.1,.01), length(seq(.2,.26,.01)) )
s = rep(seq(.2,.26,.01), length(seq(.05,.1,.01)))
M =mesh(d,s)
surf3D(x = M$x,
       y = M$y,
       z = pwr.t.test(d=M$x/M$y,power=.8,sig.level=.1,type="two.sample",alternative="two.sided")$n        ,
       colkey=FALSE,
       bty="b2",
       main="test")

how can that be done?

Upvotes: 1

Views: 383

Answers (1)

Mark Peterson
Mark Peterson

Reputation: 9570

If you are just trying to plot required sample sizes, a 3D plot may be a bit of overkill. So, here first is an approach to just making a 2D heatmap.

First, set your X and Y values. Note that these do not have to be entered multiple times (e.g. with rep)

x <- seq(.05,.1,.01)
y <- seq(.2,.26,.01)

Then, we can create a custom function to calculate the needed sample size for each combination (using outer). Of note: the effect size you had calculated before is not correct. Effect size is {mean(y) - mean(x) } / {sd(x)} (the sd should/could include y; see here for more on effect sizes). So, here, I am assuming that your standard deviation is 1, though I noted where you can put that. The pwr functions also don't appear to like playing with vectors, so I am looping over each value individually.

getN <- function(x,y){
  sapply(1:length(x), function(idx){
    pwr::pwr.t.test(d= (y[idx] - x[idx] ) # / SD???
                    , power=.8, sig.level=.1
                    , type="two.sample", alternative="two.sided"
    )$n
  })
}

Then, use outer to calculate the sample sizes for each combination. I also like to set the row/column names to help keep track of things

z <- outer(x, y, getN)
row.names(z) <- x
colnames(z) <- y

gives

           0.2      0.21     0.22     0.23     0.24     0.25     0.26
0.05  550.2098  483.6650 428.5144 382.2977 343.1846 309.7905 281.0524
0.06  631.5180  550.2098 483.6650 428.5144 382.2977 343.1846 309.7905
0.07  732.3030  631.5180 550.2098 483.6650 428.5144 382.2977 343.1846
0.08  859.3212  732.3030 631.5180 550.2098 483.6650 428.5144 382.2977
0.09 1022.5344  859.3212 732.3030 631.5180 550.2098 483.6650 428.5144
0.1  1237.1243 1022.5344 859.3212 732.3030 631.5180 550.2098 483.6650

Which is the matrix that I believe you want to plot.

Here, I am using dplyr to turn it into long format and pass it in to ggplot (and using viridis for a nice color palette, though others work fine as well)

data.frame(z, check.names = FALSE) %>%
  mutate(xVal = row.names(.)) %>%
  gather(yVal, `Sample Size`, -xVal) %>%
  ggplot(
    aes(x = xVal
        , y = yVal
        , fill = `Sample Size`)) +
  geom_tile() +
  viridis::scale_fill_viridis()

gives

enter image description here

Now, if you really want a 3D version, you can use plotly to accomplish that:

plot_ly(
  x = x
  , y = y
  , z = z
  , type = 'surface')

gives

enter image description here

And note that the produced version is interactive.

Or with plot3D:

M <- mesh(x,y)

surf3D(
  x = M$x
  , y = M$y
  , z = z
)

gives

enter image description here

And note that M is now much smaller because it is not unneccesarily duplicated.

Upvotes: 2

Related Questions