normal response

niter <- 2000
n <- 9
treatment_levels <- c("cn", "high", "high_bdnf")
insulin <- data.table(treatment = rep(treatment_levels, each=n))
X <- model.matrix(~ treatment, data=insulin)
beta <- c(0,0,0) # no effects

# the three responses are taken from the same cluster of cells and so have expected
# correlation rho. This means we need a correlated error
rho <- 0.8
Rho <- matrix(c(c(1, rho, rho), c(rho, 1, rho), c(rho, rho, 1)), nrow=3)
mu <- c(beta[1], beta[1]+beta[2], beta[1]+beta[3])

prob <- matrix(-9999, nrow=niter, ncol=3)
colnames(prob) <- c("raw", "norm", "oneway")

for(iter in 1:niter){
  fd_wide <- rmvnorm(n, mu, sigma = Rho) %>%
    data.table()
  setnames(fd_wide, old=colnames(fd_wide), new=treatment_levels)
  fd_wide[, cn.norm := cn/cn]
  fd_wide[, high.norm := high/cn]
  fd_wide[, high_bdnf.norm := high_bdnf/cn]
  fd_wide[, id := factor(paste0("id_",1:.N))]
  
  fd <- melt(fd_wide, id.vars = "id",
             measure.vars = list(treatment_levels,
                                          paste0(treatment_levels, ".norm")),
             value.name = c("insulin", "insulin.norm"))
  fd[, treatment := treatment_levels[variable]]
  fit <- lmer(insulin ~ treatment + (1|id), data = fd)
  prob[iter, "raw"] <- coef(summary(fit))["treatmenthigh", "Pr(>|t|)"]

  fit <- lm(insulin.norm ~ treatment, data = fd)
  prob[iter, "norm"] <- coef(summary(fit))["treatmenthigh", "Pr(>|t|)"]
  
  fit <- t.test(fd[treatment == "high",insulin.norm], mu=1)
  prob[iter, "oneway"] <- fit$p.value
}

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00233038
## (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00289207
## (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0100862
## (tol = 0.002, component 1)

apply(prob, 2, function(x) sum(x < 0.05)/niter)

##    raw   norm oneway 
## 0.0500 0.0165 0.0345

neg binomial response

niter <- 2000
n <- 9
treatment_levels <- c("cn", "high", "high_bdnf")
insulin <- data.table(treatment = rep(treatment_levels, each=n))
X <- model.matrix(~ treatment, data=insulin)
beta <- c(0,0,0) # no effects

# the three responses are taken from the same cluster of cells and so have expected
# correlation rho. This means we need a correlated error
rho <- 0.8
Rho <- matrix(c(c(1, rho, rho), c(rho, 1, rho), c(rho, rho, 1)), nrow=3)
mu <- c(beta[1], beta[1]+beta[2], beta[1]+beta[3])

prob <- matrix(-9999, nrow=niter, ncol=3)
colnames(prob) <- c("raw", "norm", "oneway")

for(iter in 1:niter){
  fd_wide <- rmvnorm(n, mu, sigma = Rho) %>%
    data.table()
  setnames(fd_wide, old=colnames(fd_wide), new=treatment_levels)
  fd_wide[, cn.norm := cn/cn]
  fd_wide[, high.norm := high/cn]
  fd_wide[, high_bdnf.norm := high_bdnf/cn]
  fd_wide[, id := factor(paste0("id_",1:.N))]
  
  fd <- melt(fd_wide, id.vars = "id",
             measure.vars = list(treatment_levels,
                                          paste0(treatment_levels, ".norm")),
             value.name = c("insulin", "insulin.norm"))
  fd[, treatment := treatment_levels[variable]]
  fit <- lmer(insulin ~ treatment + (1|id), data = fd)
  prob[iter, "raw"] <- coef(summary(fit))["treatmenthigh", "Pr(>|t|)"]

  fit <- lm(insulin.norm ~ treatment, data = fd)
  prob[iter, "norm"] <- coef(summary(fit))["treatmenthigh", "Pr(>|t|)"]
  
  fit <- t.test(fd[treatment == "high",insulin.norm], mu=1)
  prob[iter, "oneway"] <- fit$p.value
}

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00281995
## (tol = 0.002, component 1)

apply(prob, 2, function(x) sum(x < 0.05)/niter)

##    raw   norm oneway 
## 0.0450 0.0225 0.0360