Power computations (help)

[privefl]

I'm trying to do power calculations. I tried comparing

gap::pbsize2(
  N     = 200,
  fc    = 0.5,
  alpha = 0.05,
  gamma = 1.5,
  p     = 0.1,
  kp    = 0.1,
  model = "additive"
)

with the online tool available on http://csg.sph.umich.edu//abecasis/cats/gas_power_calculator/index.html (100 cases, 100 controls, etc.).

I do not get the same results. Maybe I'm using one of the tools wrong. Have you any input on this?

[jinghuazhao]

I need to look at it again as it was implemented quite a while ago. However, from R you can see the code as this, pbsize2 KCC The closest comparison really is through ?tscc which should be comparable.

[privefl]

OK, comparing 2 functions + online tool:

tmp <- structure(list(
  V1 = c(1, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 1.5, 1.6, 1.7, 
         1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.5, 4), 
  V2 = c(0.05, 0.062, 0.076, 0.095, 0.119, 0.147, 0.179, 0.214, 0.252, 
         0.292, 0.377, 0.463, 0.547, 0.625, 0.695, 0.756, 0.807, 0.85, 
         0.884, 0.912, 0.933, 0.95, 0.963, 0.972, 0.979, 0.996, 0.999)), 
  class = "data.frame", row.names = c(NA, -34L))

rbind(
  sapply(setNames(nm = tmp$V1), function(rr) {
    gap::tscc(
      model = "additive",
      GRR = rr,
      p1 = 0.1,
      n1 = 100,
      n2 = 100,
      M = 1,
      alpha.genome = 0.05,
      pi.samples = 1,
      pi.markers = 1,
      K = 0.1
    )$power
  })[-3, ],

  sapply(setNames(nm = tmp$V1), function(rr) {
    gap::pbsize2(
      N = 200,
      fc = 0.5,
      alpha = 0.05,
      gamma = rr,
      p = 0.1,
      kp = 0.1,
      model = "additive"
    )
  }),

  tmp$V2
)

gap::tscc gives the same result as the online tool. What am I doing wrong with gap::pbsize2 then?

[jinghuazhao]

I had a look at the pbsize2/KCC -- it is simple normal z-test of two proportions based on the same model specification as in tscc which is in line with CaTS. I gather pbsize2 would be more appropriate with common variants (AF nearer to 0.5) whereas CaTS does have genomeiwide association element -- does this agree with your findings?

[privefl]

I don't know about this. You say that both are good? But in different settings?

[jinghuazhao]

I had a bit of time from yesterday looking into issues listed in the package (the other one was done now). I was able to visit the GAS website but felt slightly puzzled if lines 400-401 of the .js (https://github.com/jenlij/GAS-power-calculator/blob/master/gas_power_calculator.js) there should have the 0.5, see for instance the paper by Evans and Purcell (2012), doi:10.1101/pdb.top069559. I had a bit of comparison with R/genetics and Bioconductor (GeneticsDesign) but was forced to drop the reference since neither is recommended by the authors nor available.