power

gistfile1.txt

require(plyr)
require(reshape2)
require(MASS)
require(ggplot2)

### Power
### Simulate data
# no of simulated datasets
nsims <- 100
# sample sizes
N <- c(3,6,9,12)
ctrl <- 2^(1:10)
todo <- expand.grid(N = N, ctrl = ctrl)
theta  <- rep(3.91, 6)  

# function to create simulated data
dosim <- function(N, mu, theta, nsims = 100){
  Nj     <- rep(N, time = 6)                # number of groups
  mus    <- rep(mu, times = Nj)             # vector of mus
  thetas <- rep(theta, times=Nj)            # vector of thetas
  x      <- factor(rep(1:6, times=Nj))      # factor
  y      <- replicate(nsims, rnegbin(sum(Nj), mus, thetas))
  return(list(x = x, y = y))
}

# simulate data
sims <- NULL
set.seed(1234)
for(i in seq_len(nrow(todo))){
  N <- todo[i, 'N']
  takectrl <- todo[i, 'ctrl']
  # reduction of 50 %
  taketrt <- takectrl * 0.5
  mu <- c(rep(takectrl, each = 2), rep(taketrt, each = 4))
  sims[[i]] <- dosim(N = N, mu = mu, nsims = nsims, theta = theta)
}

# compare methods
res <- llply(sims, function(z){
    ana <- function(y, x){
      A <- 1/min(y[y!=0])         # ln(Ax + 1) transformation
      yt <- log(A*y + 1)
      
      df <- data.frame(x, y, yt)
      
      # models
      modlm <- glm(yt ~ x, data = df)
      modglm <- glm.nb(y ~ x, data = df)
      
      # global test
      plm <- drop1(modlm, test = 'Chisq')["x", "Pr(>Chi)"]
      pglm <- drop1(modglm, test = 'Chisq')["x", "Pr(>Chi)"]
      
      return(list(A = A, 
                  plm=plm, pglm=pglm
      ))
    }
    # run models on simulated data
    res <- apply(z$y, 2, ana, x = z$x)
    res
  }, .progress = 'text'
)



##### compare methods
# global power
pow <- function(z){
  ps <- ldply(z, function(w) c(lm = w$plm, glm = w$pglm))
  apply(ps, 2, function(z) sum(z < 0.05)) / length(z)
}
pows <- ldply(res, pow)
pows$muc <- todo$ctrl
pows$N <- todo$N
powsm <- melt(pows, id.vars = c('muc', 'N'))

p2 <- ggplot(powsm) +
  geom_line(aes(y = value, x = muc, group = variable)) +
  geom_point(aes(y = value, x = muc, fill = variable), size = 4, pch = 21, color = 'black') +
  # maby try log2 transformation?
  coord_trans(xtrans = 'log2') +
  # use here
  scale_x_continuous(breaks = round(ctrl, 0)) +
  facet_wrap(~N) + 
  # axes
  labs(x = expression(mu[C]), 
       y = expression(paste('Power (global test , ', alpha, ' = 0.05)'))) +
  # appearance
  theme_bw(base_size = 12, 
           base_family = "Helvetica") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        text = element_text(size=14),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  # legend
  scale_fill_grey(name = '', 
                  breaks = c('lm', 'glm', 'pk'), 
                  labels = c('LM + log(Ay+1)', 'GLM (neg. bin.)', 'Kruskal'),
                  start = 0, end = 1) +
  theme(legend.position="bottom", legend.key = element_blank())
p2