rm(list = ls())
library(dplyr)
library(reshape2)

#--store results
mat_arry <- array(dim = c(5000, 2, 8))

#-litter icc
lit_icc <- seq(0, 0.7, .1)

for(k in 1:length(lit_icc)) {
  for(i in 1:5000){
    icc_loop <- lit_icc[k]
      v_overall <- 10
      n_litters <- 8
      pups_litter <- 4
      v_litter <- icc_loop * v_overall
      v_error <- v_overall - v_litter
      litter <- rep(1:n_litters, each = pups_litter)
      # two treatments
      treat <- rep(0:1, each = pups_litter * n_litters / 2)
      treat <- factor(treat, labels = c('C', 'T'))
      # litter effect
      litter_eff <- rnorm(n_litters, 0, sqrt(v_litter))
      # residual
      residual <- rnorm(n_litters * pups_litter, 0, sqrt(v_error))
      # the outcome measure
      y <- 5 + 0 * (treat == 'T') + litter_eff[litter] + residual
      litter <- factor(paste0('l', litter))
      my_data <- data.frame(litter, treat, y)
      mat_arry[i, 1, k] <- lmerTest::rand(lmerTest::lmer(y ~ treat + (1|litter),
                                         data = my_data))$rand.table$p.value
      m_lm <- stats::lm(y ~ treat, data = my_data)
      mat_arry[i, 2, k] <- anova(m_lm)[1, 5]
      }
}

mt_df <- melt(mat_arry)

mt_1 <- mt_df %>% filter(Var2 == 1)
mt_2 <- mt_df %>% filter(Var2 == 2)

temp <- data.frame(icc = mt_1$Var3, u_0 = mt_1$value,
                   p_lm = mt_2$value)

df_non <- temp %>% filter(u_0 < 0.05)
df_sig <- temp %>% filter(u_0 > 0.05)


t1_non <- df_non %>% group_by(icc) %>%  summarise(mean(p_lm < 0.05))
t1_sig <- df_sig %>% group_by(icc) %>%  summarise(mean(p_lm < 0.05))

results <- data.frame(cond = rep(c("sig", "non_sig"), each = 8),
                      icc = rep(seq(0, .7, .1), 2),
                      t1 = c(t1_sig$`mean(p_lm < 0.05)`,
                             t1_non$`mean(p_lm < 0.05)`))

write.csv(results, "conditional_t1_results.csv")