Skip to content

Instantly share code, notes, and snippets.

@eduardszoecs
Created September 9, 2014 15:08
Show Gist options
  • Save eduardszoecs/5bffbd9face790ac4d55 to your computer and use it in GitHub Desktop.
Save eduardszoecs/5bffbd9face790ac4d55 to your computer and use it in GitHub Desktop.
power
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment