Covariance and nonsphericity

Name the within-subject cells

Cell names combine within-factor levels with underscores. For within = c(time = 3, condition = 2), the expected names are time1_condition1, time1_condition2, time2_condition1, time2_condition2, time3_condition1, and time3_condition2.

Correlation names join a pair of cells with :.

covariance <- within_covariance(
  sd = 1,
  default_correlation = 0.5,
  correlations = c(
    "time1_condition1:time1_condition2" = 0.6,
    "time2_condition1:time2_condition2" = 0.7
  )
)

Every measurement uses the common sd; unlisted pairs use default_correlation. Pair-specific correlations can produce nonsphericity even though all diagonal variances are equal.

The package warns whenever these defaults are actually needed:

• Omitting covariance from a balanced simulation uses sd = 1 and a correlation of 0.5 for every within-subject pair.
• Calling within_covariance() without sd uses the common sd = 1.
• If correlations does not name every pair, the warning reports how many pairs are undefined. default_correlation fills only those undefined pairs; correlations explicitly listed by the user are never replaced.

Balanced power functions do not accept raw covariance matrices. Always use within_covariance() so pair names identify the intended cells explicitly; this prevents an unnamed matrix from being silently interpreted in the package’s generated within-cell order.

Use the covariance specification

Pass the object to any balanced simulation function:

power_n(
  within = c(time = 3, condition = 2),
  term = "time:condition",
  target_pes = 0.14,
  power = 0.90,
  n_sims = 5000,
  covariance = covariance,
  seed = 123
)

The custom covariance structure calibrates the target effect and generates every simulated dataset. The package projects the covariance into the tested term’s contrast space and calculates a population Greenhouse–Geisser epsilon. power_calc applies this epsilon to its degrees of freedom and noncentrality.

The exported manual helpers are currently different: design_term_means() and simulate_design_dataset() do not accept a within_covariance() object. They always use the compound-symmetric covariance defined by their sd and r arguments. Consequently, design_term_means() can return a different mean scale from a power-function call with custom covariance, because covariance changes the reference residual sum of squares used for calibration. The manual helpers therefore cannot be used to reproduce or inspect a custom-covariance power call exactly.

The sim_correction argument makes the simulated test explicit:

• "auto" is the default and preserves the package’s original behavior. It uses a Greenhouse–Geisser-corrected test when the population epsilon is below 1 - 1e-8 and ss_type is "II" or "III", and otherwise uses the uncorrected test.
• "GG" prespecifies a Greenhouse–Geisser-corrected test. Each simulated dataset uses its own sample-estimated epsilon from car::Anova().
• "none" prespecifies the uncorrected univariate test. Under nonsphericity, the package warns that the resulting excess rejection reflects inflated Type I error rather than real power.

Type I tests cannot provide GG-corrected p-values, so sim_correction = "GG" requires ss_type = "II" or "III". A between-only term or a one-degree-of-freedom within component needs no sphericity correction; in those cases a request for "GG" resolves silently to the uncorrected test.

sim_correction changes only power_sim. power_calc always applies the population epsilon and therefore always models the corrected test. If "GG" is forced under a truly spherical population, power_sim can be slightly below power_calc: correction by a sample-estimated epsilon is mildly conservative even when the population epsilon is one.

These functions estimate power for a prespecified corrected or uncorrected test. They do not simulate a conditional workflow that first applies Mauchly’s test and corrects only after rejecting sphericity. Huynh–Feldt corrections are not currently implemented.

Why the relative means pattern matters

For a within-subject term with more than one degree of freedom, the effect is a vector, theta, in a multidimensional contrast space. Its covariance in that space is Psi. Under sphericity, Psi is rotationally symmetric, so rotating a fixed-length theta does not change omnibus power. Under nonsphericity, the components have unequal weights: power depends on how theta aligns with the eigenspaces of Psi.

This answers an initially surprising question: If target_pes and the covariance are identical, how can power differ? The conventional noncentral-F calculation uses a scalar aggregate signal magnitude related to ||theta||^2 / tr(Psi). It therefore reports the same power_calc for every fixed-length direction. The actual univariate repeated-measures statistic has weighted noncentral components under nonsphericity, so power_sim is conditional on the relative mean direction. Greenhouse–Geisser correction changes the reference degrees of freedom and rejection threshold; it does not remove the underlying direction sensitivity. The uncorrected type-I path can also be direction-sensitive.

Use means_pattern() to state that relative shape. Sparse cells not listed have raw value zero. Factors outside the tested term are broadcast, and the raw values are projected exactly onto the requested term before one positive scale factor calibrates them to target_pes.

curve_pattern <- means_pattern(
  time = 1, value =  1,
  time = 2, value = -1,
  time = 3, value = -1,
  time = 4, value =  1
)

level_pattern <- means_pattern(time = "time4", value = 1)

The index time = 4 and generated level name time = "time4" are equivalent. Multiplying every raw value by the same positive constant gives the same calibrated means. Reversing all signs gives the same omnibus power, and adding an intercept or lower-order effect does not change the projected target-term component. Projection may make final target-term means nonzero in cells omitted from the sparse input.

The following validation audit used the same four-level AR(1) covariance, target_pes = 0.20, population epsilon 0.7244823, n = 25, and alpha = 0.05 for both patterns:

ar_covariance <- within_covariance(
  sd = 1,
  correlations = c(
    "time1:time2" = 0.800,
    "time1:time3" = 0.640,
    "time1:time4" = 0.512,
    "time2:time3" = 0.800,
    "time2:time4" = 0.640,
    "time3:time4" = 0.800
  )
)

power_achieved(
  within = c(time = 4), term = "time", target_pes = 0.20, n = 25,
  covariance = ar_covariance, means_pattern = curve_pattern,
  n_sims = 10000, seed = 55
)

power_achieved(
  within = c(time = 4), term = "time", target_pes = 0.20, n = 25,
  covariance = ar_covariance, means_pattern = level_pattern,
  n_sims = 10000, seed = 55
)

The observed corrected simulation powers were power_sim = 0.945 and power_sim = 0.896, respectively, while power_calc was identical. The gap is the motivation for means_pattern; it is materially larger than Monte Carlo instability at this simulation count.

The default linear/Kronecker direction

When means_pattern = NULL, the package uses a deterministic convention. For a factor with L generated levels in design order, it centers the scores 1, ..., L and normalizes them to unit Euclidean length. Before normalization:

Interaction directions are Kronecker products of the component vectors. In the package’s documented cell order, the unnormalized 2-by-2 and 2-by-4 patterns are

 1 -1       3  1 -1 -3
-1  1      -3 -1  1  3

The complete pattern is normalized once after broadcasting across factors outside the term. This convention is readable for ordered factors such as time or dose, but it is not scientifically neutral, particularly for nominal multi-level factors. Supply an explicit pattern whenever the expected shape is known. A balanced simulation warns when the implicit default is consequential: a multi-df within component and population epsilon below one. Calculation-only functions do not simulate a direction and do not issue this warning.

The two sparse-cell interfaces

means_pattern() ends each sparse definition in value, while cell_design() ends one in both n and m. Balanced designs store factor counts and generate names such as time1, so means_pattern() accepts either one-based indices or those exact generated names. In contrast, cell_design() defines an unbalanced design directly and can therefore retain arbitrary user labels. Literal cell_design() means already specify an effect direction, so power_unbalanced() needs no means_pattern argument or direction warning.

Calculation-only and unbalanced designs

Calculation-only functions accept a planned epsilon rather than a within_covariance() object; see Calculated power.

For power_unbalanced(), the common SD and correlations belong in unbalanced_covariance(); individual cells specify only sample sizes and means. See Power for unbalanced designs.