# How to do power simulations for structural equation models in R
# R code from https://philippmasur.de/blog/2020-03-11-power-simulations-sem/

library(tidyverse)
library(paramtest)

# The simulations take a while. Their results are saved next to this post
# (folder sims/) and re-used; delete the files to run them again.
cached <- function(file, expr) {
  if (file.exists(file)) return(readRDS(file))
  dir.create(dirname(file), showWarnings = FALSE)
  res <- expr
  saveRDS(res, file)
  res
}

# Power per sample size, with a 95% (Wilson) interval for the Monte Carlo error
summarise_power <- function(data, significant) {
  data |>
    group_by(n = n.test) |>
    summarise(k = sum({{ significant }}, na.rm = TRUE),
              iter = sum(!is.na({{ significant }})),
              .groups = "drop") |>
    mutate(power = k / iter,
           z = qnorm(.975),
           centre = (power + z^2 / (2 * iter)) / (1 + z^2 / iter),
           halfw = z * sqrt(power * (1 - power) / iter + z^2 / (4 * iter^2)) / (1 + z^2 / iter),
           lower = centre - halfw, upper = centre + halfw) |>
    select(n, k, iter, power, lower, upper)
}

# A smooth power curve: power is approximately Phi(a + b * sqrt(n))
power_fit <- function(pw) glm(cbind(k, iter - k) ~ sqrt(n), family = binomial("probit"), data = pw)
n_for_power <- function(fit, target = .80) unname(((qnorm(target) - coef(fit)[1]) / coef(fit)[2])^2)
power_curve <- function(fit, n) tibble(n = n, power = predict(fit, tibble(n = n), type = "response"))
fmt_n <- function(n) format(round(n, -1), big.mark = ",")

# Plot style: transparent background, muted text (readable in light and dark mode)
pal <- c(blue = "#1497d6", orange = "#d9772b", green = "#1baf7a", grey = "#8a8f98")
theme_blog <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = "transparent", colour = NA),
          panel.background = element_rect(fill = "transparent", colour = NA),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#8a8f9840", linewidth = .4),
          axis.text  = element_text(colour = "#7a7f88"),
          axis.title = element_text(colour = "#7a7f88"),
          legend.position = "none")
}

# Define the function
model_bivar <- function(simNum, n, b1) {
  
  # compute factor matrix 
  factor <- matrix(c(1, 0,    # x (no measurement error)
                     0, 1),   # y (no measurement error)
                   nrow = 2, ncol = 2, 
                   byrow = TRUE, 
                   dimnames = list(c('x', 'y'),  # rows
                                   c('x', 'y'))) # cols 
  # compute effect matrix
  y_res <- sqrt(1 - b1^2) 
  effect <- matrix(c(1, b1,    # b1 = effect size of the relationship
                     0, y_res), 
                    nrow = 2, ncol = 2, 
                    byrow = TRUE, 
                    dimnames = list(c('x', 'y'),    # rows
                                    c('x', 'y')))   # columns
  
  # generate data based on both matrices
  data <- paramtest::gen_data(factor, effect, n_cases = n)
  
  # fit lavaan model
  model <- '
    y ~ b1*x
  '
  fit <- lavaan::sem(model, data = data)
  
  # extract relevant parameters
  ests <- subset(lavaan::parameterEstimates(fit), label == "b1") 
  return(ests)
}

# Run simulations (L'Ecuyer-CMRG makes the parallel random numbers reproducible)
RNGkind("L'Ecuyer-CMRG")
set.seed(2020)
power_bivar <- cached("sims/power_bivar.rds", grid_search(
  func = model_bivar,
  params = list(n = seq(100, 1500, by = 50)),
  b1 = .1,
  n.iter = 1000,
  output = 'data.frame',
  parallel = 'snow',
  ncpus = parallel::detectCores() - 1
  ))

# result table
results(power_bivar) %>% as_tibble

# compute power
(result_bivar <- results(power_bivar) %>%
  group_by(n.test) %>%
  mutate(sig = (pvalue < .05)) %>%
  summarise(power = mean(sig, na.rm=TRUE)))

# Extract sample size for 80% power
result_bivar %>%
  filter(power >= 0.80)

# Power with Monte Carlo intervals, fitted curve and the n for 80% power
pw_bivar  <- summarise_power(results(power_bivar), pvalue < .05)
fit_bivar <- power_fit(pw_bivar)
n80_bivar <- n_for_power(fit_bivar)

alt_bivar <- paste0("Simulated power curve for a manifest correlation of .10, with 95% intervals around each estimate and a fitted curve. ",
                    "Power reaches 80% at about n = ", fmt_n(n80_bivar), ".")

ggplot(pw_bivar, aes(n, power)) +
  geom_hline(yintercept = .80, linetype = "dashed", colour = pal["grey"]) +
  annotate("segment", x = n80_bivar, xend = n80_bivar, y = 0, yend = .80,
           linetype = "dotted", colour = pal["grey"]) +
  geom_line(data = power_curve(fit_bivar, 100:1500), colour = pal["blue"], linewidth = 1) +
  geom_pointrange(aes(ymin = lower, ymax = upper), colour = pal["blue"], size = .2) +
  annotate("text", x = n80_bivar + 30, y = .74, hjust = 0, colour = "#7a7f88",
           label = paste0("80% power at n ≈ ", fmt_n(n80_bivar))) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(x = "Sample size", y = "Power") +
  theme_blog()

library(pwr)
pwr.r.test(r = c(.1),
           sig.level = .05,
           power = .80,
           alternative = "two.sided")

# define function
model_sem <- function(simNum,
                      n, 
                      b1, 
                      x_ind1, x_ind2, x_ind3, 
                      y_ind1, y_ind2, y_ind3) {
  
  # compute factor matrix (in this case with measurement error)
  factor <- matrix(c(x_ind1, 0,   # y1
                     x_ind2, 0,   # y2
                     x_ind3, 0,   # y3
                     0, y_ind1,   # y1
                     0, y_ind2,   # y2
                     0, y_ind3),  # y3
                   nrow=6, ncol=2, 
                   byrow=TRUE, 
                   dimnames=list(c('x1', 'x2', 'x3', 'y1', 'y2', 'y3'),  # rows (observed)
                                 c('x', 'y')))                           # columns (latent)
  
  # matrix of the effects structure (variance-covariance)
  y_res <- sqrt(1 - b1^2)
  effect <- matrix(c(1, b1,     # x
                     0, y_res), # y
                    nrow=2, ncol=2, 
                    byrow=TRUE, 
                    dimnames=list(c('x', 'y'),    # rows
                                  c('x', 'y')))   # columns
  
  # generate data based on both matrices
  data <- paramtest::gen_data(factor, effect, n_cases = n)

  # fit lavaan model
  model <- '
   # measurement model
   x =~ x1 + x2 + x3
   y =~ y1 + y2 + y3

   # structural regression model
   y ~ b1*x
  '
  fit <- lavaan::sem(model, data = data)
  ests <- lavaan::parameterEstimates(fit)
  
  # extract relevant parameters
  est <- ests[ests$label == 'b1', 'est']
  p <- ests[ests$label == 'b1', 'pvalue']
  return(c(est = est, 
           p = p, 
           sig = (p < .05)))
}

# Measurement error specified by factor loading = .7
RNGkind("L'Ecuyer-CMRG")
set.seed(2021)
power_sem1 <- cached("sims/power_sem1.rds", grid_search(
  func = model_sem, 
  params = list(n = seq(500, 4000, by = 100)),
  b1 = .1,
  n.iter = 1000, 
  output = 'data.frame',
  x_ind1 = .7,    
  x_ind2 = .7, 
  x_ind3 = .7,      
  y_ind1 = .7, 
  y_ind2 = .7, 
  y_ind3 = .7, 
  parallel = 'snow',
  ncpus = parallel::detectCores() - 1
  ))

results(power_sem1) %>% as_tibble

# Let's also test a more realistic example (factor loadings = .5)
set.seed(2022)
power_sem2 <- cached("sims/power_sem2.rds", grid_search(
  func = model_sem, 
  params = list(n = seq(500, 4000, by = 100)),
  b1 = .1,
  n.iter = 1000, 
  output = 'data.frame',
  x_ind1 = .5,    # changed to .5 
  x_ind2 = .5, 
  x_ind3 = .5,      
  y_ind1 = .5, 
  y_ind2 = .5, 
  y_ind3 = .5, 
  parallel = 'snow',
  ncpus = parallel::detectCores() - 1
  ))

# Extract result table
(result_sem1 <- results(power_sem1) %>%
    mutate(rel = .7) %>%
    bind_rows(results(power_sem2) %>%
                mutate(rel = .5)) %>%
  group_by(rel, n.test) %>%
  summarise(power=mean(sig, na.rm=TRUE)))

pw_sem <- bind_rows(
  summarise_power(results(power_sem1), sig == 1) |> mutate(loadings = ".70"),
  summarise_power(results(power_sem2), sig == 1) |> mutate(loadings = ".50"))
fits_sem <- lapply(split(pw_sem, pw_sem$loadings), power_fit)
n80_sem  <- sapply(fits_sem, n_for_power)
curves_sem <- bind_rows(lapply(names(fits_sem), \(l) power_curve(fits_sem[[l]], 500:4000) |> mutate(loadings = l)))

alt_sem <- paste0("Simulated power curves for a correlation of .10 between two latent variables, for item factor loadings of .70 and .50, with 95% intervals and fitted curves. ",
                  "80% power is reached at about n = ", fmt_n(n80_sem[[".70"]]), " with loadings of .70 and only at about n = ", fmt_n(n80_sem[[".50"]]), " with loadings of .50.")

ggplot(pw_sem, aes(n, power, colour = loadings)) +
  geom_hline(yintercept = .80, linetype = "dashed", colour = pal["grey"]) +
  geom_line(data = curves_sem, linewidth = 1) +
  geom_pointrange(aes(ymin = lower, ymax = upper), size = .2) +
  annotate("text", x = 4000, y = c(.20, .12), hjust = 1, colour = pal[c("blue", "orange")],
           label = c("loadings .70 (ω ≈ .74)", "loadings .50 (ω = .50)")) +
  scale_colour_manual(values = c(".70" = pal[["blue"]], ".50" = pal[["orange"]])) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(x = "Sample size", y = "Power") +
  theme_blog()

model_med <- function(simNum,
                      n, 
                      b1, 
                      b2,
                      b3,
                      x_ind1, x_ind2, x_ind3,
                      m_ind1, m_ind2, 
                      y_ind1, y_ind2, y_ind3) {
  
  # compute factor matrix (in this case with measurement error)
  factor <- matrix(c(x_ind1, 0, 0,   # x1
                     x_ind2, 0, 0,   # x2
                     x_ind3, 0, 0,   # x3
                     0, m_ind1, 0,   # m1
                     0, m_ind2, 0,   # m2
                     0, 0, y_ind1,   # y1
                     0, 0, y_ind2,   # y2
                     0, 0, y_ind3),  # y3
                   nrow=8, ncol=3, 
                   byrow=TRUE, 
                   dimnames=list(c('x1', 'x2', 'x3', 'm1', 'm2', 'y1', 'y2', 'y3'),  # rows (observed)
                                 c('x', 'm', 'y')))                           # columns (latent)
  
  # matrix of the effects structure: gen_data() computes latent = Z %*% effect,
  # with Z independent standard normal. To get m = b1*x + e_m and
  # y = b2*x + b3*m + e_y (all latent variances = 1), the weights must be:
  m_res <- sqrt(1 - b1^2)          # m's own (residual) part
  y_x   <- b2 + b1 * b3            # weight of x in y = direct + indirect effect
  y_m   <- b3 * m_res              # weight of m's residual part in y
  y_res <- sqrt(1 - y_x^2 - y_m^2) # y's own (residual) part
  effect <- matrix(c(1, b1,    y_x,   # x
                     0, m_res, y_m,   # m
                     0, 0,     y_res),# y
                    nrow=3, ncol=3, 
                    byrow=TRUE, 
                    dimnames=list(c('x', 'm', 'y'),    # rows
                                  c('x', 'm', 'y')))   # columns
  
  # generate data based on both matrices
  data <- paramtest::gen_data(factor, effect, n_cases = n)

  # fit lavaan model
  model <- '
   # measurement model
   x =~ x1 + x2 + x3
   y =~ y1 + y2 + y3
   m =~ m1 + m2

   # structural regression model
   m ~ b1*x
   y ~ b2*x + b3*m

   ind := b1*b3
   total := b2 + ind
  '
  fit <- lavaan::sem(model, data = data)
  ests <- lavaan::parameterEstimates(fit)
  
  # extract relevant parameters
  b1_est     <- ests[ests$label == 'b1', 'est']
  b1_p       <- ests[ests$label == 'b1', 'pvalue']
  b2_est     <- ests[ests$label == 'b2', 'est']
  b2_p       <- ests[ests$label == 'b2', 'pvalue']
  b3_est     <- ests[ests$label == 'b3', 'est']
  b3_p       <- ests[ests$label == 'b3', 'pvalue']
  ind_est    <- ests[ests$label == 'ind', 'est']
  ind_p      <- ests[ests$label == 'ind', 'pvalue']
  total_est  <- ests[ests$label == 'total', 'est']
  total_p    <- ests[ests$label == 'total', 'pvalue']

  return(c(b1_est   = b1_est,    b1_p    = b1_p,
           b2_est   = b2_est,    b2_p    = b2_p,
           b3_est   = b3_est,    b3_p    = b3_p,
           ind_est  = ind_est,   ind_p   = ind_p,
           total_est= total_est, total_p = total_p))
}

RNGkind("L'Ecuyer-CMRG")
set.seed(2023)
power_med <- cached("sims/power_med.rds", grid_search(
  func = model_med, 
  params = list(n = seq(700, 4500, by = 100)),
  b1 = .1,  # assuming a small effect for all three paths       
  b2 = .1,
  b3 = .1,
  n.iter = 1000, 
  output = 'data.frame',
  x_ind1 = .7,    
  x_ind2 = .7, 
  x_ind3 = .7, 
  m_ind1 = .7,    
  m_ind2 = .7,
  y_ind1 = .7, 
  y_ind2 = .7, 
  y_ind3 = .7, 
  parallel = 'snow',
  ncpus = parallel::detectCores() - 1
  ))

# Extract result table
(result_med <- results(power_med) %>%
  group_by(n.test) %>%
    mutate(b1_sig = (b1_p < .05),
           b2_sig = (b2_p < .05),
           b3_sig = (b3_p < .05),
           ind_sig = (ind_p < .05),
           total_sig = (total_p < .05)) %>%
  summarise(b1_power = mean(b1_sig, na.rm=TRUE),
            b2_power = mean(b2_sig, na.rm=TRUE),
            b3_power = mean(b3_sig, na.rm=TRUE),
            ind_power = mean(ind_sig, na.rm=TRUE),
            total_power = mean(total_sig, na.rm=TRUE)))

pw_med <- results(power_med) |>
  pivot_longer(ends_with("_p"), names_to = "effect", values_to = "p") |>
  mutate(effect = sub("_p$", "", effect)) |>
  group_by(effect) |>
  group_modify(~ summarise_power(.x, p < .05)) |>
  ungroup()
fits_med <- lapply(split(pw_med, pw_med$effect), power_fit)
n80_med  <- sapply(fits_med, n_for_power)
curves_med <- bind_rows(lapply(names(fits_med), \(e) power_curve(fits_med[[e]], 700:4500) |> mutate(effect = e)))

alt_med <- paste0("Simulated power curves for the latent mediation model with 95% intervals and fitted curves. ",
                  "The three direct paths (all .10) reach 80% power at about n = ", fmt_n(min(n80_med[c("b1", "b2", "b3")])), " to ", fmt_n(max(n80_med[c("b1", "b2", "b3")])),
                  ", the total effect earlier, and the indirect effect only at about n = ", fmt_n(n80_med[["ind"]]), ".")

# Three direct paths in one colour, the indirect and total effect highlighted
type <- c(b1 = "direct", b2 = "direct", b3 = "direct", ind = "indirect", total = "total")
cols <- c(direct = pal[["blue"]], indirect = pal[["orange"]], total = pal[["green"]])
labs_med <- tibble(type = c("direct", "total", "indirect"),
                   label = c("direct paths a, b, c′", "total effect c′ + a × b", "indirect effect a × b"),
                   y = c(.30, .22, .14))

ggplot(mutate(pw_med, type = type[effect]), aes(n, power, colour = type, group = effect)) +
  geom_hline(yintercept = .80, linetype = "dashed", colour = pal["grey"]) +
  geom_line(data = mutate(curves_med, type = type[effect]), linewidth = .9) +
  geom_pointrange(aes(ymin = lower, ymax = upper), size = .15, alpha = .6) +
  geom_text(data = labs_med, aes(x = 4500, y = y, label = label, colour = type),
            inherit.aes = FALSE, hjust = 1) +
  scale_colour_manual(values = cols) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(x = "Sample size", y = "Power") +
  theme_blog()

sessionInfo()