Power Simulation for Longitudinal Mixed Model

This script simulates longitudinal data for a mixed effects model with repeated ratings over time and delay, fits the target model with lme4, and reports the effect size for the time effect in several metrics.

The target model is:

rating ~ time * delay +
  (1 + time | subj) +        # random intercept and slope for time by subject
  (1 | subj:mem_id)          # random intercept for memories nested within subjects

The script provides:

  1. A data simulation function: simulate_longitudinal()
  2. A helper function to extract effect sizes from the fitted model: report_time_effect()

Dependencies

The script uses the following R packages:

library(lme4)
library(MASS)
library(effectsize)

Install them if needed:

install.packages(c("lme4", "MASS", "effectsize"))

Functions

1. simulate_longitudinal()

This function simulates data for a longitudinal design with:

  • Subjects (subj)
  • Three delay conditions (short, mid, long)
  • Multiple memories per delay within each subject
  • Two time points: initial and followup

Signature:

simulate_longitudinal(
  n_subj = 30,
  n_memories_per_delay = 14,
  beta_0 = 0,
  beta_time = 0.5,
  beta_delay_mid = 0.0,
  beta_delay_long = 0.0,
  beta_time_delay_mid = 0.0,
  beta_time_delay_long = 0.0,
  sd_subj_intercept = 0.8,
  sd_subj_time = 0.4,
  cor_subj = 0.2,
  sd_mem = 0.5,
  sd_resid = 1.0
)

Arguments (key ones):

  • n_subj: Number of participants.
  • n_memories_per_delay: Number of memories per delay per participant.
  • beta_0: Fixed intercept.
  • beta_time: Main effect of time (followup vs initial at short delay).
  • beta_delay_mid, beta_delay_long: Main effects of delay at initial time.
  • beta_time_delay_mid, beta_time_delay_long: Time by delay interaction terms.
  • sd_subj_intercept: SD of subject random intercepts.
  • sd_subj_time: SD of subject random slopes for time.
  • cor_subj: Correlation between random intercept and random slope.
  • sd_mem: SD of random intercepts for memories nested within subjects.
  • sd_resid: Residual SD.

What it does:

  • Prints the expected within subject effect size (Cohen’s dz) for the time effect, given beta_time, sd_subj_time, and sd_resid:

    Expected within-subject time effect size (Cohen's dz) ~ ...

  • Returns a data frame dat containing:

    • subj, delay, mem_id, time
    • Dummy variables (time_num, delay_mid, delay_long)
    • Random effects (b0, b_time, b_mem)
    • Linear predictor mu
    • Simulated outcome rating

2. report_time_effect()

This function takes a fitted model and prints a short summary of the time effect in three metrics:

  • Raw fixed effect (unstandardized coefficient)
  • Standardized coefficient
  • Partial eta squared

Signature:

report_time_effect <- function(mod) {
  ...
}

Input:

  • mod: An lmerMod object fitted to the simulated data with the same naming structure for time:

    • The fixed effect for followup is named timefollowup.
    • The standardized parameter table has a row time [followup].
    • The partial eta squared table has a row time.

What it prints:

Example output:

Your effect size is: raw = 0.556 , standardized = 0.39 , partial eta² = 0.63

This is:

  • raw: the unstandardized estimate for timefollowup
  • standardized: the standardized coefficient for the time contrast (followup vs initial)
  • partial eta²: partial eta squared for the time effect from eta_squared(mod, partial = TRUE)

Example Workflow

# Load packages
library(lme4)
library(MASS)
library(effectsize)

# Define helper function
report_time_effect <- function(mod) {
  # raw fixed effect for time (followup vs initial)
  raw <- fixef(mod)["timefollowup"]

  # standardized coefficient for time [followup]
  std_tab <- standardize_parameters(mod)
  std_col <- grep("Std", names(std_tab), value = TRUE)  # e.g., "Std. Coef."
  std <- as.numeric(std_tab[std_tab$Parameter == "time [followup]", std_col])

  # partial eta² for time
  eta_tab <- eta_squared(mod, partial = TRUE)
  eta_col <- grep("Eta2", names(eta_tab), value = TRUE) # e.g., "Eta2 (partial)"
  eta <- as.numeric(eta_tab[eta_tab$Parameter == "time", eta_col])

  cat(
    "Your effect size is:",
    "raw =", round(raw, 3), ",",
    "standardized =", round(std, 3), ",",
    "partial eta² =", round(eta, 3), "
"
  )
}

# Define simulation function (as above)
# simulate_longitudinal <- function(...) { ... }

# 1. Simulate data for 30 participants
set.seed(123)
dat <- simulate_longitudinal(n_subj = 30)

# 2. Fit the mixed model
m_long <- lmer(
  rating ~ time * delay +
    (1 + time | subj) +
    (1 | subj:mem_id),
  data = dat
)

# 3. Inspect model and effect sizes
summary(m_long)
eta_squared(m_long, partial = TRUE)
standardize_parameters(m_long)

# 4. Report time effect in one line
report_time_effect(m_long)

Modifying Sample Size and Effect Size

To explore different scenarios:

  • Change the number of participants:

    dat <- simulate_longitudinal(n_subj = 50)

  • Change the strength of the time effect:

    dat <- simulate_longitudinal(n_subj = 30, beta_time = 0.3)

Then refit the model and call:

m_long <- lmer(
  rating ~ time * delay +
    (1 + time | subj) +
    (1 | subj:mem_id),
  data = dat
)

report_time_effect(m_long)

The function will always report the effect sizes for the time contrast in a consistent format, so you can quickly compare scenarios across different sample sizes and population parameters.