`design_mfa()` fits a trial-level multivariate design model and decomposes each estimable design effect into low-rank spatial components. Subjects are represented by fixed intercepts, with optional ridge-pooled term deviations, and trial weights are normalized within subject by default so unequal trial counts do not silently change population weighting.
Usage
design_mfa(Y, ...)
# Default S3 method
design_mfa(Y, ...)
# S3 method for class 'multidesign'
design_mfa(Y, design = NULL, formula, subject, ...)
# S3 method for class 'list'
design_mfa(
Y,
design,
formula,
subject = NULL,
nuisance = NULL,
contrasts = NULL,
terms = NULL,
reference_grid = NULL,
reference_weight = c("equal_cells", "observed"),
trial_weights = NULL,
block_weight = c("equal_subject", "precision"),
feature_names = NULL,
feature_map = NULL,
rank = NULL,
ss_type = c("partial", "sequential"),
term_scope = NULL,
alignment = c("shared", "soft"),
subject_model = c("fixed_intercept", "ridge_slopes"),
subject_penalty = 1,
random_terms = NULL,
spatial_penalty = 0,
spatial_laplacian = NULL,
alignment_penalty = 0,
store = c("sufficient", "fitted", "residuals"),
control = design_mfa_control(),
...
)Arguments
- Y
A list of subject matrices, or a `multidesign` object. Rows are trials. Columns are shared features under shared alignment and may be subject-specific under soft alignment.
- ...
Reserved for future methods.
- design
A list of data frames with one row per corresponding trial. Omit this when `Y` is a `multidesign` object.
- formula
A one-sided formula defining design effects.
- subject
Optional unique subject identifiers, one per list element.
- nuisance
Optional one-sided formula defining nuisance columns. Terms must not also appear in `formula`.
- contrasts
Optional named list of factor contrast specifications. The aliases `"poly"`, `"sum"`, `"treatment"`, and `"helmert"` are supported.
- terms
Optional character vector selecting formula terms to return.
- reference_grid
Optional data frame on which effects are evaluated. The default is the Cartesian product of declared factor levels and modest numeric grids.
- reference_weight
Either `"equal_cells"`, `"observed"`, or a nonnegative numeric vector with one value per reference-grid row.
- trial_weights
Optional list of nonnegative trial-precision vectors.
- block_weight
Population weighting: `"equal_subject"`, `"precision"`, or a numeric subject-weight vector.
- feature_names
Optional shared feature schema, required when `Y` has no column names. Under soft alignment these name the reference features.
- feature_map
For soft alignment, a list of sparse `p_s x p_0` maps taking reference fields into each subject's feature space.
- rank
Optional retained component rank, either a scalar or values named by term. By default all numerically nonzero effect components are retained.
- ss_type
`"partial"` for containing-versus-reduced effects or `"sequential"` for formula-order effects.
- term_scope
For partial effects, remove only assigned term columns or the term and all descendant interactions. The default is `"hierarchical"` for partial effects and `"assigned"` for sequential effects.
- alignment
`"shared"` for exact feature correspondence or `"soft"` for user-supplied sparse reference-to-subject maps.
- subject_model
`"fixed_intercept"` or centered `"ridge_slopes"`.
- subject_penalty
Positive ridge penalty for centered subject-specific deviations when `subject_model = "ridge_slopes"`.
- random_terms
Optional formula terms receiving subject deviations. The default is every term in `formula`.
- spatial_penalty
Nonnegative loading-smoothness strength.
- spatial_laplacian
Optional sparse feature Laplacian. Required when `spatial_penalty > 0`; named dimensions are aligned to `feature_names`.
- alignment_penalty
Optional reference-coefficient ridge stabilization for the soft map solver. It must be zero for shared alignment.
- store
`"sufficient"` stores compact component representations; `"fitted"` or `"residuals"` additionally stores reference effect matrices.
- control
A named list returned by [design_mfa_control()].
Value
A `design_mfa` object containing the weighted design fit, an estimability report, and one `design_mfa_effect` object per requested term. Training responses are retained in `fit_spec$refit$data` for subject-clustered refits; the stored model does not keep the dead `q x p` `XtY` product.
Details
For every requested term, the method compares a containing weighted model with a reduced model, evaluates their prediction difference on an explicit reference grid, and performs a weighted thin SVD of that effect matrix. Thus, the components maximize fitted variance associated with a named design term rather than total variance in the activation fields.
Partial sums of squares need not add across terms in an unbalanced design. When `ss_type = "sequential"`, results are explicitly order-dependent. Ordered factors should use planned contrasts, such as `"poly"`. Effects of observed response variables such as confidence are conditional associations, not causal experimental effects.
Shared alignment supports either fixed subject intercepts or centered ridge-pooled subject slopes. Soft alignment uses explicit sparse maps from a reference field into each subject's feature space; it currently supports the fixed-intercept subject model.
The exported function name is retained for API continuity, but the statistical method is an ASCA+/partial-RDA-style design decomposition. It does not perform the block normalization that defines multiple factor analysis (MFA).
Fitted objects retain the training responses in `fit_spec$refit$data` so subject-clustered refits do not need the caller to pass data again. `saveRDS(fit)` therefore writes that dataset. Peak working memory during the fit is approximately `O(p q)` plus the data already in memory.
Examples
set.seed(1)
subjects <- paste0("S", 1:6)
design <- lapply(subjects, function(id) {
data.frame(
saliency = ordered(rep(c(20, 40, 60, 80, 100), each = 6)),
probe = factor(rep(c("old", "new", "lure"), 10))
)
})
pattern <- seq(-1, 1, length.out = 8)
Y <- lapply(design, function(d) {
outer(as.numeric(d$saliency) - 3, pattern) +
matrix(rnorm(nrow(d) * 8, sd = 0.2), nrow(d), 8)
})
names(Y) <- names(design) <- subjects
feature_names <- paste0("voxel_", 1:8)
fit <- design_mfa(
Y,
design,
~ saliency + probe,
contrasts = list(saliency = "poly", probe = "sum"),
feature_names = feature_names
)
components(fit, "saliency")
#> component singular_value explained
#> 1 1 2.61248542 0.9992295230
#> 2 2 0.05361734 0.0004208894
#> 3 3 0.03571102 0.0001867076
#> 4 4 0.03335455 0.0001628800
continuous_design <- lapply(design, function(d) {
d$trial_index <- seq_len(nrow(d))
d
})
continuous_fit <- design_mfa(
Y, continuous_design, ~ trial_index,
feature_names = feature_names
)
head(effect_scores(continuous_fit, "trial_index"))
#> Component_1
#> [1,] -4.285504
#> [2,] -3.989952
#> [3,] -3.398848
#> [4,] -3.103296
#> [5,] -2.512192
#> [6,] -2.216640