An a priori power analysis of a mixed model design from the terHorst lab. Just for fun.

casey_mixed_mod_power.R

#Some libraries that will help
library(nlme)
library(dplyr)
library(tidyr)
library(mvtnorm)

######This controls it all!
n_sims <- 5

#A function to make a sigma matrix
#but, assumes a standard normal - sd of 1, 
#and thus correlation = covriance. You can monkey here
make_sigma <- function(cor_effect, len){
  ret <- matrix(rep(cor_effect, len^2), ncol=len)
  diag(ret) <- 1
  
  return(ret)
}

#OK, start by making a data frame with effect sizes and correlation in tanks
#then use that to add data
test_df <- expand.grid(effect_size = rep(1:5, n_sims), 
                              cor_tank = seq(0,0.3, length.out=4)) %>%
  mutate(sim = 1:n()) %>%
   #add tanks 
   crossing(tank = 1:4) %>%
   #add subsample
   crossing(subsample = 1:10) %>%
  #check that we're ok
 # group_by(sim, tank, effect_size, cor_tank) %>%
  # summarize(len = n())
  
  #add data
  group_by(sim, tank) %>%
  mutate(sim_obs = rmvnorm(1, mean = effect_size*tank, 
                           sigma=make_sigma(cor_tank[1], 10))) %>%
  
  ungroup()

#Let's just run one analysis
one_test_df <- test_df %>% filter(sim==8) #schoose a number with a nonzero cor_tank


one_lme <- lme(sim_obs ~ tank, random = ~1|factor(tank), data=one_test_df)
summary(one_lme)

broom::tidy(one_lme)
broom::tidy(one_lme, "fixed")


#let's run a lot of models
fit_df <- test_df %>%
  group_by(sim, effect_size, cor_tank) %>%
  nest() %>%
  mutate(fit_mod = purrr::map(data, ~lme(sim_obs ~ tank, random = ~1|factor(tank), data=.))) %>%
  mutate(coef = purrr:::map(fit_mod, ~broom::tidy(., "fixed"))) %>%
  unnest(coef) %>%
  filter(term != "(Intercept)")

#let's see what this looks like

library(ggplot2)
ggplot(data=fit_df, mapping = aes(x=cor_tank, y=p.value, color=factor(effect_size))) +
  geom_point(position = position_dodge(width=0.1), size=2, alpha=0.5)


ggplot(data=fit_df, mapping = aes(x=cor_tank, y=p.value, color=factor(effect_size))) +
  geom_jitter( size=2, alpha=0.5) +
  facet_wrap(~effect_size) +
  scale_y_log10() +
  geom_hline(yintercept=0.05, lty=2) #show an alpha


#### OK - power time for different alphas!
alphas = seq(0.01, 0.1, length.out=10)

pow_df <- fit_df %>% 
  crossing(alpha = alphas) %>%
  group_by(alpha, effect_size, cor_tank) %>%
  summarize(type_2_error_rate = sum(p.value>alpha)/n()) %>%
  ungroup() %>%
  mutate(power = 1-type_2_error_rate)


#visualize power
ggplot(data = pow_df, mapping=aes(x=cor_tank, y=power, color=factor(alpha))) +
  geom_line() + geom_point() +
  facet_wrap(~effect_size)