The Population

The residual variances below are chosen so every variable has unit variance, so the paths read as standardized effects: $a = 0.4$, $b = 0.35$, $c' = 0.15$, and therefore $ab = 0.14$.

pop_med <- "
  x ~~ 1*x
  m ~ 0.4*x
  m ~~ 0.84*m
  y ~ 0.35*m + 0.15*x
  y ~~ 0.813*y
"
cov_sem(pop_med)$sigma_theta
#>      m    y    x
#> m 1.00 0.41 0.40
#> y 0.41 1.00 0.29
#> x 0.40 0.29 1.00

The Analysis Model and the Parameters of Interest

The analysis model is free; the labels name the parameters of interest, and ab := a*b defines the indirect effect (its standard error comes from the delta method, as in lavaan itself).

med_model <- "
  m ~ a*x
  y ~ b*m + cp*x
  ab := a*b
"

Composite Power at a Candidate Sample Size

Suppose $N = 100$ is under consideration. parameters selects the labeled set; G is the number of Monte Carlo replications, kept at 25 here so the document knits quickly, and the seed makes the result reproducible. As the note at the top says, a plan you intend to defend is worth G = 1000 or more.

med_at_100 <- ss_power_composite_sem(
  model = med_model, pop_model = pop_med,
  parameters = c("a", "b", "ab"),
  N = 100, G = 25, seed = 113)
med_at_100

compare(med_at_100, ref$med_at_100,
        c("composite_power", "power_a", "power_b", "power_ab"))
#>                 G = 25 G = 10000
#> composite_power   0.84    0.8692
#> power_a           0.96    0.9883
#> power_b           0.96    0.9367
#> power_ab          0.84    0.8692

The power_a, power_b, and power_ab rows are the marginal powers, each the proportion of the G replications in which that parameter was significant. The composite_power row is the proportion in which all three were significant in the same replication: here 0.84, at most the smallest marginal power, and estimated with a simulation standard error of about 0.073. Because the three tests share one fitted model they are dependent, so the composite need not equal the product of the marginals; the simulation gets the joint probability right without any independence assumption.

The Necessary Sample Size

Planning replaces N with desired_power:

med_plan <- ss_power_composite_sem(
  model = med_model, pop_model = pop_med,
  parameters = c("a", "b", "ab"),
  desired_power = 0.80, G = 25, seed = 113)
med_plan

compare(med_plan, ref$med_plan,
        c("necessary_N", "composite_power", "power_a", "power_b", "power_ab"))
#>                 G = 25 G = 10000
#> necessary_N      76.00   88.0000
#> composite_power   0.80    0.8052
#> power_a           1.00    0.9762
#> power_b           0.96    0.9067
#> power_ab          0.80    0.8052

A sample of $N = 88$ is the smallest at which the estimated composite power reaches 0.80, reading the reference column. The result carries the same broom summary as the rest of the ss_power_* family:

generics::tidy(med_plan)
#>          term estimate power
#> 1 sample_size       76   0.8

The Weakest Parameter Governs the Design

The direct path $c' = 0.15$ was deliberately left out of the set above. Adding it shows why the choice of the set is a substantive decision, not a formality:

med_plan_cp <- ss_power_composite_sem(
  model = med_model, pop_model = pop_med,
  parameters = c("a", "b", "cp", "ab"),
  desired_power = 0.80, G = 25, seed = 113)
med_plan_cp

compare(med_plan_cp, ref$med_plan_cp,
        c("necessary_N", "composite_power", "power_cp"))
#>                 G = 25 G = 10000
#> necessary_N     353.00  343.0000
#> composite_power   0.84    0.8143
#> power_cp          0.84    0.8143

Requiring the small direct path to be significant as well moves the necessary sample size from $N = 88$ to $N = 343$. The composite is bounded by its weakest member, so the set should contain exactly the parameters the paper’s conclusion requires, and each addition is a design commitment with a visible price.

Accuracy for the Set: AIPE

When the research questions concern magnitudes, the goal is a sufficiently narrow confidence interval for every parameter of interest. desired_width states the full width per parameter, as a named vector so a width can never silently attach to the wrong parameter; the indirect effect, on its smaller scale, is held to a narrower interval here.

med_aipe <- ss_aipe_composite_sem(
  model = med_model, pop_model = pop_med,
  parameters = c("a", "b", "ab"),
  desired_width = c(a = 0.25, b = 0.25, ab = 0.15),
  G = 25, seed = 113)
med_aipe

Confidence level: 95%

compare(med_aipe, ref$med_aipe,
        c("necessary_N", "composite_assurance", "mean_width_a",
          "mean_width_b", "mean_width_ab"))
#>                       G = 25 G = 10000
#> necessary_N         245.0000  238.0000
#> composite_assurance   0.5200    0.3816
#> mean_width_a          0.2271    0.2330
#> mean_width_b          0.2495    0.2499
#> mean_width_ab         0.1254    0.1301

With no assurance, the criterion is the expected width: the returned $N = 238$ is the smallest at which the mean simulated width of every interval is within its target. Widths vary from sample to sample, so a study of that size obtains all three sufficiently narrow intervals in only about 38 percent of its realizations (the composite_assurance row). Supplying an assurance plans against that joint event directly:

med_aipe_80 <- ss_aipe_composite_sem(
  model = med_model, pop_model = pop_med,
  parameters = c("a", "b", "ab"),
  desired_width = c(a = 0.25, b = 0.25, ab = 0.15),
  assurance = 0.80, G = 25, seed = 113)
med_aipe_80

Confidence level: 95%

compare(med_aipe_80, ref$med_aipe_80,
        c("necessary_N", "composite_assurance"))
#>                     G = 25 G = 10000
#> necessary_N         262.00  268.0000
#> composite_assurance   0.92    0.8082

A planning summary an author could report: with $N = 268$, all three intervals are simultaneously no wider than their targets (0.25, 0.25, and 0.15) in an estimated 92 percent of studies, under the stated population model and a 95% confidence level.

The Population Growth Process

The posited population: average status 5 at wave 1, average growth 0.3 per year, intercept variance 1, slope variance 0.2, intercept-slope covariance $-0.15$ (a correlation of about $-0.34$: those who start higher grow less), and residual variance 0.5 at every wave. Every parameter, including every intercept and latent mean, is fixed.

pop_lgm <- "
  i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  i ~~ 1*i
  s ~~ 0.2*s
  i ~~ -0.15*s
  t1 ~~ 0.5*t1; t2 ~~ 0.5*t2; t3 ~~ 0.5*t3; t4 ~~ 0.5*t4
  t1 ~ 0*1; t2 ~ 0*1; t3 ~ 0*1; t4 ~ 0*1
  i ~ 5*1
  s ~ 0.3*1
"
lgm_pop <- cov_sem(pop_lgm)
lgm_pop$mu_theta
#>  t1  t2  t3  t4 
#> 5.0 5.3 5.6 5.9
lgm_pop$sigma_theta
#>      t1   t2   t3   t4
#> t1 1.50 0.85 0.70 0.55
#> t2 0.85 1.40 0.95 1.00
#> t3 0.70 0.95 1.70 1.45
#> t4 0.55 1.00 1.45 2.40

cov_sem() returns both moments the population implies: the wave means rise by 0.3 per year, and the wave variances and covariances follow from the growth factor variances, their covariance, and the residuals.

The population is worth seeing. Each light line below is one person’s model implied trajectory, drawn from the population distribution of $(i, s)$; the bold line is the average trajectory.

set.seed(113)
growth_factors <- MASS::mvrnorm(
  n = 40, mu = c(5, 0.3),
  Sigma = matrix(c(1, -0.15, -0.15, 0.2), 2, 2))
traj <- data.frame(
  person = rep(seq_len(40), each = 4),
  wave = rep(1:4, times = 40),
  wellbeing = growth_factors[rep(seq_len(40), each = 4), 1] +
    growth_factors[rep(seq_len(40), each = 4), 2] * rep(0:3, times = 40))

ggplot(traj, aes(wave, wellbeing, group = person)) +
  geom_line(color = "grey70", linewidth = 0.4) +
  geom_line(data = data.frame(person = 0, wave = 1:4,
                              wellbeing = 5 + 0.3 * 0:3),
            color = unname(grDevices::palette.colors(2)[2]),
            linewidth = 1.4) +
  labs(title = "The Population the Design Is Planned On",
       subtitle = "Individual latent trajectories; the bold line is the average trajectory",
       x = "Wave", y = "Well-being")

plot of chunk lgm-trajectories

The negative intercept-slope covariance is visible: trajectories that start high tend to tilt down relative to the average, so the fan narrows over time.

The Analysis Model

The analysis model is the same growth structure with its parameters free. Written for lavaan::sem(), the growth parameterization fixes the observed intercepts to zero and frees the latent means (this is exactly what lavaan::growth() does internally). The two parameters of interest carry labels.

lgm_model <- "
  i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  i ~~ cov_is*s
  t1 ~ 0*1; t2 ~ 0*1; t3 ~ 0*1; t4 ~ 0*1
  i ~ 1
  s ~ mu_s*1
"

Composite Power for the Growth Questions

First, what does a candidate $N = 150$ deliver?

lgm_at_150 <- ss_power_composite_sem(
  model = lgm_model, pop_model = pop_lgm,
  parameters = c("mu_s", "cov_is"),
  N = 150, G = 25, seed = 113)
lgm_at_150

compare(lgm_at_150, ref$lgm_at_150,
        c("composite_power", "power_mu_s", "power_cov_is"))
#>                 G = 25 G = 10000
#> composite_power   0.64      0.67
#> power_mu_s        1.00      1.00
#> power_cov_is      0.64      0.67

The average growth of 0.3 per year is easy to detect (power_mu_s is 1), and the composite is governed almost entirely by the covariance question: power_cov_is is 0.64, and the composite power is 0.64. Planning for the pair:

lgm_plan <- ss_power_composite_sem(
  model = lgm_model, pop_model = pop_lgm,
  parameters = c("mu_s", "cov_is"),
  desired_power = 0.80, G = 25, seed = 113)
lgm_plan

compare(lgm_plan, ref$lgm_plan,
        c("necessary_N", "composite_power", "power_cov_is"))
#>                 G = 25 G = 10000
#> necessary_N     216.00  201.0000
#> composite_power   0.88    0.8008
#> power_cov_is      0.88    0.8008

About $N = 201$ participants are needed for both growth questions to be answered affirmatively in the same study with probability 0.80, under the stated population. A researcher who planned only for the slope mean, the headline effect, would have chosen a far smaller study and then usually failed to detect the covariance that the theory also requires.

Accuracy for the Growth Parameters

The magnitude of growth and of its dependence on initial status are usually the quantities of substantive interest, so the AIPE plan holds the slope mean to an interval no wider than 0.15 (that is, $\pm 0.075$ around an effect of 0.3 per year) and the covariance to 0.25, jointly, in 80 percent of studies:

lgm_aipe <- ss_aipe_composite_sem(
  model = lgm_model, pop_model = pop_lgm,
  parameters = c("mu_s", "cov_is"),
  desired_width = c(mu_s = 0.15, cov_is = 0.25),
  assurance = 0.80, G = 25, seed = 113)
lgm_aipe

Confidence level: 95%

compare(lgm_aipe, ref$lgm_aipe,
        c("necessary_N", "composite_assurance",
          "width_within_desired_mu_s", "width_within_desired_cov_is"))
#>                             G = 25 G = 10000
#> necessary_N                 220.00  220.0000
#> composite_assurance           0.92    0.8197
#> width_within_desired_mu_s     0.92    0.8197
#> width_within_desired_cov_is   1.00    0.9956

The accuracy goal needs $N = 220$, and the width_within_desired_* rows show which target binds: the slope mean’s interval is the harder one to keep narrow at this $N$. When both existence and magnitude matter, the defensible design uses the larger of the power and AIPE sample sizes, here $N = 220$.