Type I Error

gistfile1.txt

### Type1 Error
require(plyr)
require(reshape2)
require(MASS)
require(ggplot2)
nsims <- 100
# sample sizes
N <- c(3,6,9,12)
ctrl <- 2^(1:10)
# both as grid
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(12345)
for(i in seq_len(nrow(todo))){
  N <- todo[i, 'N']
  takectrl <- todo[i, 'ctrl']
  # control = treatment
  taketrt <- takectrl * 1
  mu <- c(rep(takectrl, each = 2), rep(taketrt, each = 4))
  sims[[i]] <- dosim(N = N, mu = mu, nsims = nsims, theta = theta)
}

resfoo <- 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
}
res <- llply(sims, resfoo, .progress = 'text')



t1 <- 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)
}
t1s <- ldply(res, t1)
t1s$muc <- todo$ctrl
t1s$N <- todo$N
t1sm <- melt(t1s, id.vars = c('muc', 'N'))

pt1 <- ggplot(t1sm) +
  geom_line(aes(y = value, x = muc, group = variable)) +
  geom_point(aes(y = value, x = muc, fill = variable), size = 4, pch = 21, color = 'black') +
  coord_trans(xtrans = 'log2') +
  scale_x_continuous(breaks = round(ctrl, 0)) +
  facet_wrap(~N) + 
  # axes
  labs(x = expression(mu[C]), 
       y = expression(paste('Type 1 error (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())
pt1