2 conditions pupil size example (bad defaults)

demo.Rmd

---
title: 'Example'
output: html_document
---



```{r, message=FALSE}
# load required packages
library("lme4")        # model specification / estimation
library("afex")        # anova and deriving p-values from lmer
library("broom.mixed") # extracting data from model fits 
library("faux")        # data simulation
# NOTE: to install the 'faux' package, use:
# devtools::install_github("debruine/faux")
library("tidyverse")   # data wrangling and visualisation

# ensure this script returns the same results on each run
set.seed(8675309)
```


### Data simulation function


```{r}
# set up the custom data simulation function
my_sim_data <- function(nsubj  = 100, # number of subjects
                        ntrials = 130, # number of trials per condition
                        b0     = 800, # grand mean
                        b1     =  50, # effect of condition
                        S0s_sd = 100, # by-subject random intercept sd
                        S1s_sd =  40, # by-subject random slope sd
                        scor   = 0.2, # correlation between intercept and slope
                        err_sd = 200  # residual (standard deviation)
                        ) {
  # simulate subjects
  subjects <- faux::sim_design(
    within = list(effect = c(S0s = "By-subject random intercepts", 
                             S1s = "By-subject random slopes")), 
    n = nsubj,
    sd = c(S0s_sd, S1s_sd), 
    r = scor,
    id = "subj_id",
    plot = FALSE
  )
  
  # simulate trials
  dat_sim <- crossing(subj_id = subjects$subj_id,
                      trial_id = 1:ntrials,
                      condition = c("A", "B")) %>%
    inner_join(subjects, "subj_id") %>%
    mutate(err = rnorm(nrow(.), mean = 0, sd = err_sd),
           cond = recode(condition, "A" = -0.5, "B" = 0.5)) %>%
    mutate(pupil_size = b0 + S0s + (b1 + S1s) * cond + err) %>%
    select(subj_id, condition, cond, pupil_size)
  
  dat_sim
}
```

## Analyzing the simulated data

```{r}
# fit a linear mixed-effects model to data
mod_sim <- lmer(pupil_size ~ 1 + cond + (1 + cond | subj_id),
                data = my_sim_data(), REML = TRUE)

summary(mod_sim, corr = FALSE)
```



## Calculate Power

```{r}
# set up the power function
my_lmer_power <- function(...) {
  # ... is a shortcut that forwards any arguments to my_sim_data()
  dat_sim <- my_sim_data(...)
  mod_sim <- lmer(pupil_size ~ 1 + cond + (1 + cond | subj_id),
                data = my_sim_data(), REML = TRUE)
  
  broom.mixed::tidy(mod_sim)
}
```


```{r}
# run one model with default parameters
my_lmer_power()
```


```{r}
# run one model with new parameters
my_lmer_power(ntrials = 50, b1 = 20)
```


```{r, eval = FALSE}
# run simulations and save to a file
reps <- 100
sims <- purrr::map_df(1:reps, ~my_lmer_power())
write_csv(sims, "sims/sims.csv")
```

```{r, message=FALSE}
# read saved simulation data
sims <- read_csv("sims/sims.csv")
```

```{r}
# calculate mean estimates and power for specified alpha
alpha <- 0.05

sims %>% 
  filter(effect == "fixed") %>%
  group_by(term) %>%
  summarise(
    mean_estimate = mean(estimate),
    mean_se = mean(std.error),
    power = mean(p.value < alpha)
  )
```