Logo

Exercises – Monte Carlo Power Simulation for Mixed-Effects Logistic Regression, and Some Parallelization

Basics of R for Data Science

(Invented) Scenario: Does Smartphone Distraction Reduce Trial-Level Accuracy?

Imagine you are studying sustained attention in PhD students during a computerized Go/No-Go task.
Each participant completes trials in two conditions:

  • Control: no notifications (baseline).
  • Distraction: intermittent smartphone pings while performing the task.

The response is binary (1 = correct, 0 = incorrect). Each participant goes through many trials and provides several responses. Therefore, a mixed-effects logistic regression is appropriate:

  • Random intercepts account for individual differences in general correctness (they can be set in any case in which there is more than one single response per participant);
  • Random slopes account for individual differences in the effect/treatment (some participants are more susceptible to the treatment, other are less; random slopes can be set in cases where each participant is observed in more than one condition and provide more than one response in each condition… which is indeed the present case!)

We will:

  • Part 1: Simulate one single dataset for this scenario (trial-level, repeated within participants) to understand how it works;
  • Part 2: Estimate power by simulation for a single design (a given N, a given effect size) using the traditional significance rule (p < .05);
  • Part 3: Repeat the power estimate using parallelized iterations with the future.apply package.

Part 1: Simulate one single dataset

Assumptions for the data-generating process:

  • Participants: n_participants = 60
  • Trials per condition per participant: n_trials = 40 (so 80 trials per person)
  • Baseline accuracy (Control) ≈ 0.73 → beta0 ≈ logit(0.73) ≈ 1.0
  • Distraction reduces accuracy: beta1 = -0.30 (moderate impairment)
  • Between-participant variability in baseline accuracy: sd_participant = 0.8 (random intercept SD)

First of all, install (if needed) and load the lme4 package for fitting mixed-effects models:

# install.packages("lme4") # if needed
library(lme4)

Second, let’s define all design parameters:

n_participants  = 60
n_trials        = 40
beta0           = 1.0    # baseline log-odds in Control
beta1           = -0.30   # effect of Distraction (negative -> lower accuracy)
sd_participant  = 0.8    # random intercept SD
alpha           = 0.05   # significance threshold

Third, let’s build one dataset:

# Initialize empty dataset, one row per participant per condition
df = expand.grid(
  participant = factor(1:n_participants),
  condition   = c("Control", "Distraction")
)
# Repeat each row times trials
df = df[rep(1:nrow(df),each=n_trials) , ]
df$trial = rep(1:n_trials, times=n_participants*2)

# Generate random intercept per participant
rand_int = rnorm(n_participants, 0, sd_participant)
rand_slo = rnorm(n_participants, 0, sd_participant)
df$randInt_i    = rand_int[df$participant]
df$randSlo_i    = rand_slo[df$participant]

# Dummy coding for simulation
df$cond_dummy = ifelse(df$condition == "Distraction", 1, 0)

# Log-odds and probabilities
df$logit     = beta0 + beta1*df$cond_dummy + 
                df$randInt_i + df$randSlo_i*df$cond_dummy
df$prob       = plogis(df$logit)

# Simulate binary responses
df$response = rbinom(nrow(df), size = 1, prob = df$prob)

Inspect the first few rows of the dataset:

head(df)
  participant condition trial randInt_i randSlo_i cond_dummy logit  prob
1           1   Control     1      1.01     0.287          0  2.01 0.882
2           1   Control     2      1.01     0.287          0  2.01 0.882
3           1   Control     3      1.01     0.287          0  2.01 0.882
4           1   Control     4      1.01     0.287          0  2.01 0.882
5           1   Control     5      1.01     0.287          0  2.01 0.882
6           1   Control     6      1.01     0.287          0  2.01 0.882
  response
1        1
2        1
3        1
4        0
5        1
6        1

Fit the planned analysis model:

fit = glmer(response ~ condition + (condition | participant),
                  data = df, family = binomial(link="logit"),
                  control = glmerControl(optimizer = "bobyqa",
                                         optCtrl = list(maxfun = 1e5)))

See model summary:

summary(fit)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: response ~ condition + (condition | participant)
   Data: df
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  5574.9   5607.2  -2782.4   5564.9     4795 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5079 -0.9407  0.4845  0.6603  1.8236 

Random effects:
 Groups      Name                 Variance Std.Dev. Corr
 participant (Intercept)          0.4710   0.6863       
             conditionDistraction 0.2374   0.4872   0.50
Number of obs: 4800, groups:  participant, 60

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.96864    0.10122   9.570  < 2e-16 ***
conditionDistraction -0.24201    0.09359  -2.586  0.00971 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
cndtnDstrct 0.041

Interpretation: The coefficient for condition estimates how much "Distraction" affects the response accuracy relative to "Control" on the logit scale. A negative, significant estimate supports the “distraction impairs accuracy” hypothesis.

Part 2: Monte Carlo Power Simulation (with replicate)

We keep the same design and true parameters as above. We repeatedly:

  1. Simulate a new dataset;
  2. Fit the same glmer(...) model;
  3. Record the p-value for the (fixed) effect of interest.

Below we will run only \(10\) iterations to avoid too long computations… But in real life you will certainly need many more iterations for precision.

First of all, let’s wrap the whole simulation of one dataset inside a custom function, which returns the p-value associated with the effect of interest, which is condition:

simulate_once = function(n_participants = 60,
                         n_trials = 40,
                         beta0 = 1,
                         beta1 = -0.30,
                         sd_participants = 0.8) {
  
  # Random intercept per participant (new draw each experiment)
  rand_int = rnorm(n_participants, 0, sd_participant)
  rand_slo = rnorm(n_participants, 0, sd_participant)

  # Build the design
  df = expand.grid(
    participant = factor(1:n_participants),
    condition   = c("Control", "Distraction"),
    trial       = 1:n_trials
  )
  df$cond_dummy = ifelse(df$condition == "Distraction", 1, 0)
  df$randInt_i    = rand_int[df$participant]
  df$randSlo_i    = rand_slo[df$participant]
  
  df$logit   = beta0 + beta1*df$cond_dummy + 
                df$randInt_i + df$randSlo_i*df$cond_dummy
  
  df$prob        = plogis(df$logit)
  df$response = rbinom(nrow(df), 1, df$prob)
  
  # Fit model; use try() to handle possible non-convergences
  fit = try(glmer(response ~ condition + (condition | participant),
                   data = df, family = binomial(link="logit"),
                   control = glmerControl(optimizer = "bobyqa",
                                          optCtrl = list(maxfun = 1e5))),
             silent = TRUE)
  if (inherits(fit, "try-error")) return(NA)
  
  coefs = summary(fit)$coefficients
  pval = coefs["conditionDistraction", "Pr(>|z|)"]

  return(pval)
}

Before we run the Monte Carlo simulation, some best practices:

  • set a random seed for reproducibility;
  • define and store the number of iterations;
  • initialize vectors of results.
set.seed(0)
nsim = 10
pvals = rep(NA, nsim)

RUN!!! let’s launch replicate() times nsim and get all the p-values we want!

pvals = replicate(nsim,simulate_once())

Finally, some post-processing: remove NAs if any (non-converged runs), and see results.

We define power as the proportion of times where \(p < 0.05\):

pvals_clean = pvals[is.finite(pvals)]
prop_converged = length(pvals_clean) / nsim
power_pvalue   = mean(pvals_clean < 0.05)

list(
  simulations_requested = nsim,
  simulations_converged = length(pvals_clean),
  convergence_rate      = round(prop_converged, 3),
  power_p_less_0_05     = round(power_pvalue, 3)
)
$simulations_requested
[1] 10

$simulations_converged
[1] 10

$convergence_rate
[1] 1

$power_p_less_0_05
[1] 0.7

Note. Mixed-effects models can occasionally fail to converge under certain parameterizations, that’s why the use the try() function, and also track and report the proportion of successful fits. Not using try() is dangerous: Think what happens if after \(1,000\) iterations and a full day of computation the whole script breaks because of one single error due to non-convergence!

Part 3: Speeding up with parallel computing

We keep the same simulate_once() and parameters as in Part 2.

To parallelize in a way that should work on Windows, macOS, and Linux, we use the future.apply package with a multisession plan (each worker is its own R session)

⚠️ The below code was actually run on Windows… it needs to be actually checked on other operating systems.

First, some setup for parallel power simulations using future.apply:

# preliminarily install packages if needed
library(future)
library(future.apply)
library(parallel)   # for detectCores()

# Choose workers, leave one core free for safety 🙂
workers = max(1, detectCores(logical = TRUE) - 1)

# Cross-platform backend: multisession should work on Windows+macOS+Linux
plan(multisession, workers = workers)

# Define a number of iterations
nsim = 10

RUN!!! in parallel 🙂

pvals_parallel = future_sapply(
  1:nsim,
  function(i) simulate_once(),   # uses objects defined above (e.g., beta0, beta1, etc.)
  future.seed = 0        # reproducible, independent RNG streams per worker
)

Finally, some post-processing identical and see results, just like in Part 2:

pvals_clean     = pvals_parallel[is.finite(pvals_parallel)]
prop_converged  = length(pvals_clean) / nsim
power_pvalue    = mean(pvals_clean < 0.05)

list(
  backend               = sprintf("future.apply::multisession (%d workers)", workers),
  simulations_requested = nsim,
  simulations_converged = length(pvals_clean),
  convergence_rate      = round(prop_converged, 3),
  power_p_less_0_05     = round(power_pvalue, 3)
)
$backend
[1] "future.apply::multisession (11 workers)"

$simulations_requested
[1] 10

$simulations_converged
[1] 10

$convergence_rate
[1] 1

$power_p_less_0_05
[1] 0.7