Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement