|
| 1 | +function C = extractCategoryFactors(C, opts) |
| 2 | +%EXTRACTCATEGORYFACTORS Reduce each model/category to a few factors and store them. |
| 3 | +% C = EXTRACTCATEGORYFACTORS(C) takes the FULL expanded corpus from |
| 4 | +% READANNOTATIONCORPUSFULL and runs factor analysis SEPARATELY within each category: |
| 5 | +% each of the 13 multivariate model outputs (SigLIP embedding, DINOv2 embedding, |
| 6 | +% CLAP embedding, AudioSet tags, action posteriors, EmoNet, GoEmotions text-emotion, |
| 7 | +% SigLIP/CLAP probes, MFCC, chroma, facial affect, text sentiment) AND the block of |
| 8 | +% interpretable scalar features. For every category it extracts factor-score time |
| 9 | +% series (one per factor, over all timepoints) and attaches them to C as |
| 10 | +% C.extracted_factors, so the ~2.7k raw variables collapse to a compact, labeled set |
| 11 | +% of per-model factors usable as fMRI regressors. |
| 12 | +% |
| 13 | +% Method: features are z-scored; the factor model is fit from the pairwise |
| 14 | +% correlation matrix (nearest-positive-definite repaired), so audio-only clips still |
| 15 | +% contribute their audio/language factors. FACTORAN (max-likelihood, varimax) is used |
| 16 | +% where the block's degrees of freedom allow; otherwise it falls back to PCA (typical |
| 17 | +% for the high-dimensional opaque embeddings). Per-timepoint scores use the regression |
| 18 | +% method with mean-imputation, so every timepoint gets a score. |
| 19 | +% |
| 20 | +% Name-value options: |
| 21 | +% "NumFactors" (10) target factors per category (capped by block dof/rank) |
| 22 | +% "Rotate" ("varimax") rotation for the FACTORAN path |
| 23 | +% "Method" ("auto") "auto" | "factoran" | "pca" |
| 24 | +% "Save" ("") if set, save the extracted_factors struct to this .mat |
| 25 | +% "Verbose" (true) |
| 26 | +% |
| 27 | +% Adds C.extracted_factors with fields: |
| 28 | +% .scores [T x F] all factor scores concatenated across categories |
| 29 | +% .labels [F x k] table: FactorName, Category, Model, Class, Method, VarExpl |
| 30 | +% .byCategory struct array, one per category, with: |
| 31 | +% name, model, class, kind, nVars, nFactors, method, |
| 32 | +% loadings [nVars x m], varExplained [1 x m], scores [T x m], |
| 33 | +% factorNames [1 x m], vars (FEATUREINFO subset) |
| 34 | +% .stim, .time_sec, .ids row keys (same rows as C.X) |
| 35 | +% |
| 36 | +% Example: |
| 37 | +% C = readAnnotationCorpusFull("annotations/corpus"); |
| 38 | +% C = extractCategoryFactors(C, "NumFactors", 10, "Save", "analysis/extracted_factors.mat"); |
| 39 | +% % SigLIP visual factors as a time series for one clip: |
| 40 | +% r = C.stim=="BigBuckBunny"; s = C.extracted_factors.byCategory; |
| 41 | +% si = s(strcmp({s.name},"siglip_embedding")); |
| 42 | +% plot(C.time_sec(r), si.scores(r,1:3)); |
| 43 | +% |
| 44 | +% See also READANNOTATIONCORPUSFULL, FACTORANALYSISCORPUS, FEATUREINFO. |
| 45 | + |
| 46 | +arguments |
| 47 | + C (1,1) struct |
| 48 | + opts.NumFactors (1,1) double {mustBePositive, mustBeInteger} = 10 |
| 49 | + opts.Rotate (1,1) string = "varimax" |
| 50 | + opts.Method (1,1) string {mustBeMember(opts.Method, ["auto","factoran","pca"])} = "auto" |
| 51 | + opts.Save (1,1) string = "" |
| 52 | + opts.Verbose (1,1) logical = true |
| 53 | +end |
| 54 | + |
| 55 | +info = C.info; X = C.X; T = size(X,1); |
| 56 | + |
| 57 | +% ---- define categories: the 13 multivariate channels + the interpretable-scalar block |
| 58 | +vecLeaves = unique(info.Leaf(info.Dtype == "vector"), "stable"); |
| 59 | +cats = struct("name",{}, "model",{}, "class",{}, "kind",{}, "cols",{}); |
| 60 | +for i = 1:numel(vecLeaves) |
| 61 | + idx = find(info.Leaf == vecLeaves(i)); |
| 62 | + kind = "interpretable-vector"; |
| 63 | + if info.IsEmbedding(idx(1)), kind = "embedding"; |
| 64 | + elseif ismember(vecLeaves(i), ["audioset_tags","action_posteriors"]), kind = "taxonomy"; end |
| 65 | + cats(end+1) = struct("name", vecLeaves(i), "model", info.Model(idx(1)), ... |
| 66 | + "class", info.Class(idx(1)), "kind", kind, "cols", idx(:)'); %#ok<AGROW> |
| 67 | +end |
| 68 | +sIdx = find(info.Dtype ~= "vector" & info.Numeric); |
| 69 | +cats(end+1) = struct("name","interpretable_scalars", "model","(mixed)", ... |
| 70 | + "class","(mixed)", "kind","interpretable-scalar", "cols", sIdx(:)'); |
| 71 | + |
| 72 | +% ---- run FA/PCA per category |
| 73 | +byCat = struct([]); allScores = []; FN=strings(0,1); Cat=strings(0,1); |
| 74 | +Mdl=strings(0,1); Cls=strings(0,1); Meth=strings(0,1); VE=zeros(0,1); |
| 75 | +for i = 1:numel(cats) |
| 76 | + ca = cats(i); |
| 77 | + cols = ca.cols; |
| 78 | + sub = X(:, cols); |
| 79 | + vinfo = info(cols, :); |
| 80 | + |
| 81 | + % z-score; drop zero-variance columns |
| 82 | + mu = mean(sub,1,"omitnan"); sd = std(sub,0,1,"omitnan"); |
| 83 | + keep = sd > 0 & isfinite(sd); |
| 84 | + sub = sub(:, keep); vinfo = vinfo(keep,:); mu = mu(keep); sd = sd(keep); |
| 85 | + Z = (sub - mu) ./ sd; |
| 86 | + d = size(Z,2); |
| 87 | + if d == 0, continue; end |
| 88 | + |
| 89 | + % pairwise correlation -> nearest PD; drop fully-collinear columns |
| 90 | + R = corr(Z, "rows","pairwise"); R(~isfinite(R)) = 0; R(1:d+1:end) = 1; |
| 91 | + [R, vinfo, Z] = localDropCollinear(R, vinfo, Z); |
| 92 | + R = localNearestPD(R); |
| 93 | + d = size(R,1); |
| 94 | + nEff = size(Z,1); |
| 95 | + |
| 96 | + m = min(opts.NumFactors, max(1, d-1)); |
| 97 | + method = opts.Method; |
| 98 | + if method == "auto" |
| 99 | + method = "factoran"; |
| 100 | + if d < 3 || m > localMaxFactoranM(d), method = "pca"; end |
| 101 | + end |
| 102 | + |
| 103 | + Zi = Z; Zi(isnan(Zi)) = 0; % standardized mean-imputation |
| 104 | + if method == "factoran" && m > localMaxFactoranM(d) |
| 105 | + m = max(1, localMaxFactoranM(d)); |
| 106 | + end |
| 107 | + try |
| 108 | + if method == "factoran" |
| 109 | + [L, ~, ~, ~] = factoran(R, m, "Xtype","covariance", "Nobs",nEff, "Rotate",char(opts.Rotate)); |
| 110 | + W = R \ L; |
| 111 | + scores = Zi * W; |
| 112 | + ve = 100 * sum(L.^2,1) / d; |
| 113 | + else |
| 114 | + error("use pca"); % jump to catch/pca path |
| 115 | + end |
| 116 | + catch |
| 117 | + method = "pca"; |
| 118 | + [V, ev] = localTopEig(R, m); |
| 119 | + L = V .* sqrt(max(ev,0))'; % component loadings |
| 120 | + scores = Zi * V; % PC scores on standardized data |
| 121 | + ve = 100 * ev' / d; |
| 122 | + end |
| 123 | + |
| 124 | + % order factors by variance explained (desc) |
| 125 | + [ve, si] = sort(ve, "descend"); L = L(:,si); scores = scores(:,si); |
| 126 | + fnames = ca.name + "_F" + string(1:m); |
| 127 | + |
| 128 | + e = numel(byCat) + 1; |
| 129 | + byCat(e).name = ca.name; byCat(e).model = ca.model; byCat(e).class = ca.class; |
| 130 | + byCat(e).kind = ca.kind; byCat(e).nVars = d; byCat(e).nFactors = m; |
| 131 | + byCat(e).method = string(method); byCat(e).loadings = L; byCat(e).varExplained = ve; |
| 132 | + byCat(e).scores = scores; byCat(e).factorNames = fnames; byCat(e).vars = vinfo; |
| 133 | + |
| 134 | + allScores = [allScores, scores]; %#ok<AGROW> |
| 135 | + FN=[FN; fnames(:)]; Cat=[Cat; repmat(ca.name,m,1)]; Mdl=[Mdl; repmat(string(ca.model),m,1)]; |
| 136 | + Cls=[Cls; repmat(string(ca.class),m,1)]; Meth=[Meth; repmat(string(method),m,1)]; VE=[VE; ve(:)]; |
| 137 | + |
| 138 | + if opts.Verbose |
| 139 | + fprintf(" %-22s %-9s d=%4d -> %2d factors (%s) var%%=%s\n", ... |
| 140 | + ca.name, "["+ca.kind+"]", d, m, method, join(compose("%.0f",ve),"/")); |
| 141 | + end |
| 142 | +end |
| 143 | + |
| 144 | +EF.scores = allScores; |
| 145 | +EF.labels = table(FN, Cat, Mdl, Cls, Meth, VE, 'VariableNames', ... |
| 146 | + {'FactorName','Category','Model','Class','Method','VarExplained'}); |
| 147 | +EF.byCategory = byCat; |
| 148 | +EF.stim = C.stim; EF.time_sec = C.time_sec; EF.ids = C.ids; |
| 149 | +C.extracted_factors = EF; |
| 150 | + |
| 151 | +if opts.Verbose |
| 152 | + fprintf("Extracted %d factors across %d categories (%d timepoints).\n", ... |
| 153 | + size(allScores,2), numel(byCat), T); |
| 154 | +end |
| 155 | +if opts.Save ~= "" |
| 156 | + extracted_factors = EF; |
| 157 | + save(opts.Save, "extracted_factors", "-v7.3"); |
| 158 | + fprintf("wrote %s\n", opts.Save); |
| 159 | +end |
| 160 | +end |
| 161 | + |
| 162 | +% ------------------------------------------------------------------------- |
| 163 | +function [R, info, Z] = localDropCollinear(R, info, Z) |
| 164 | +p = size(R,1); keep = true(p,1); |
| 165 | +for i = 2:p |
| 166 | + if any(abs(R(i,1:i-1)) > 0.999 & keep(1:i-1)'), keep(i) = false; end |
| 167 | +end |
| 168 | +R = R(keep,keep); info = info(keep,:); Z = Z(:,keep); |
| 169 | +end |
| 170 | + |
| 171 | +function A = localNearestPD(A) |
| 172 | +A = (A + A')/2; [V,D] = eig(A); d = diag(D); d(d < 1e-6) = 1e-6; |
| 173 | +A = V*diag(d)*V'; A = (A + A')/2; s = sqrt(diag(A)); A = A ./ (s*s'); |
| 174 | +A = (A + A')/2; A(1:size(A,1)+1:end) = 1; |
| 175 | +end |
| 176 | + |
| 177 | +function [V, ev] = localTopEig(R, m) |
| 178 | +[V,D] = eig((R+R')/2); ev = diag(D); |
| 179 | +[ev, ix] = sort(ev, "descend"); V = V(:, ix); |
| 180 | +m = min(m, size(V,2)); V = V(:,1:m); ev = ev(1:m); |
| 181 | +end |
| 182 | + |
| 183 | +function mmax = localMaxFactoranM(d) |
| 184 | +% largest m>=1 satisfying FACTORAN's dof rule (d-m)^2 >= d+m |
| 185 | +mmax = 0; |
| 186 | +for m = 1:d-1 |
| 187 | + if (d-m)^2 >= d+m, mmax = m; end |
| 188 | +end |
| 189 | +end |
0 commit comments