Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Type I Error
### 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.