Skip to content

Develop - #2

Merged
RWKrause merged 12 commits into
mainfrom
develop
Sep 8, 2026
Merged

RWKrause merged 12 commits into
mainfrom
develop

Conversation

@RWKrause

@RWKrause RWKrause commented Sep 7, 2026

Copy link
Copy Markdown
Owner

No description provided.

RWKrause and others added 12 commits August 26, 2026 16:18
.Rproj.user held 12 tracked RStudio state files. test.RData is a 41 MB
local artefact that .Rbuildignore kept out of the tarball but nothing kept
out of a commit. Also de-duplicates the doubled ^doc$/^Meta$ entries.
…thods

QAPglm() and QAPcss() shared 19 of 21 arguments and ~90% duplicated
permutation code. They are replaced by a single QAP() that detects CSS
data from the dimensions of the dependent variable.

- multi-mode data: each dimension its own node set, permuted independently
  (RMPerm, make_qap_data, make_css_data, residuals_to_array)
- direct OLS path bypassing the formula interface for gaussian fits,
  ~9x faster, verified bit-identical against lm() by the golden master
- batched GPU permutations via torch, replacing the per-permutation loop
- standard extractors: coef, summary, confint, vcov, nobs, fitted,
  residuals, logLik, as.data.frame, all on a shared "QAP" parent class

Fixes that change results: mode = "undirected" was accepted and never
read, so every dyad entered twice; robust SEs fed deviance residuals into
an OLS HC3 formula for every non-gaussian family. array_to_vector() looped
over dimension 1 while slicing dimension 3, truncating non-cubic arrays.

BREAKING CHANGE: QAPglm() and QAPcss() are removed; nullhyp defaults to
"qapspp" everywhere.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
- DESCRIPTION said "MIT + file LICENSE" while LICENSE held the full GPL-3
  text. Set License: GPL (>= 3) and drop the LICENSE file, which CRAN
  ships centrally and must not accompany a GPL declaration.
- inst/CITATION.txt -> inst/CITATION so citation("MrQAP") finds it
- add method references with DOIs, URL and BugReports, shorten Title
- drop LazyData (there is no data/ directory)
- NEWS.md and an expanded README covering the removed functions and the
  two fixes that change results
- vignettes/mrqap.Rmd: one narrative across matrix, CSS and 2-mode data

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Multinomial:

QAP(family = "multinom") errored outright on its default nullhyp.
qap_init_pmats() sized the result matrices with length(coefficients) --
(ncat-1)*k for a matrix -- instead of ncol(), and read names() of a matrix,
which is NULL, so the per-predictor assignment out$lower[, xi] failed.

Where it did run (nullhyp = "qapy"), the p-values were computed and then
lost: qap_coef_table() returned NA unless the matrix had exactly two rows,
so summary() and as.data.frame() reported nothing, and print.QAPGLM()
hardcoded row 2 -- the second CATEGORY's coefficient comparisons -- and
labelled it Pr(<=t), recycling a length-k vector into a (ncat-1)*k table.

The p-value matrices stack a block of b-comparisons above a block of
t-comparisons. Consumers now take the half, not the row, which is correct
for the ordinary two-row case as well. Rows are labelled b:<cat>/t:<cat>
and both nullhyp paths produce identically shaped matrices. confint() now
says why permutation intervals do not exist for multinomial fits instead
of blaming less_mem.

GMM:

gmm() was started from rnorm(), redrawn on every permutation. On a
well-scaled design this was harmless, but on a wider one the moment
conditions have several local optima: 25 random starts on the same data
gave x1 estimates from -80.6 to +184.8 (sd 24.6) where the truth was 0.6.
That noise went straight into the null distribution. gmm_start() supplies
the corresponding GLM estimate, with a neutral nuisance start for negbin
and zip; the same design now gives 0.528.

Adds test-multinom.R (48 assertions) and test-gmm-models.R (107), the
latter covering recovery, seed reproducibility, effect-vs-null separation,
qapspp, MDS errors and the extractors for all four GMM families.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
combine_qap_estimates():

- checked only class, nullhyp, comparison names and dimnames, so two fits
  to DIFFERENT networks pooled without a word. It now also requires the
  family, estimator, mode, diagonal, data-shape flags and reference to
  match, and compares the coefficients: the baseline fit does not depend
  on the permutations, so coefficients that differ prove the runs are not
  poolable whatever their settings claim.
- never pooled $null_dist, so confint() on a combined object used run 1's
  draws while reporting the summed reps. The draws are now stacked.
- weighted by reps, but the proportions have denominator n_valid, which is
  smaller whenever permutations failed. Fits now record n_valid -- a
  scalar under qapy, one entry per predictor under qapspp -- and pooling
  weights column by column.
- a single fit passed alone is itself a list, so it was read as a list of
  models and failed incomprehensibly well past the length check.

validate_qap_input() compared dim() but never dimnames(), so matrices
built from separate sources in different node orders produced a silently
meaningless model. Where both sides label a dimension the labels must now
agree; unlabelled data is unaffected.

Also folds in the single-pass aggregation from the performance phase:
aggregate_perm_results() built a 4 * n_valid element list and compared its
names three times to organise three running sums.

The combine fixtures in test-extended.R and test-helpers.R had drifted
from the real object -- no "QAP" parent class, one class name misspelt
"QAPregression", none of the fields the checks read. Updated, with the
existing drift guard extended to catch it next time.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…user

Under qapspp only one column of the design changes per permutation, and
only that column's statistics are ever read -- compare_perm_to_baseline()
discards the rest. The fast path nonetheless ran a fresh qr() over the
whole design every replication, computing p-1 coefficients to drop them.

qap_fwl_template() decomposes the untested columns once; qap_fwl_solve()
then gets the tested coefficient from two inner products, and the full
model's residual sum of squares as rss_y - b^2 * sxx. For robust errors it
uses the leverage decomposition h = h_Z + x_t^2/sxx, so HC3 specialises to
sum(x_t^2 * omega)/sxx^2. All exact algebra, held to lm(), to the full
solver and to HC3() at 1e-10 in test-fwl.R.

Measured, n = 60, 3 predictors, 400 reps, best of 3:

  qapspp plain    1.882 -> 0.480 ms/rep   (3.9x)
  qapspp robust   2.928 -> 0.524 ms/rep   (5.6x)

which puts qapspp at parity with qapy per replication. The golden master
is unchanged, as it must be: the RNG stream is untouched.

Progress reporting: the package called progressr::with_progress() itself,
making it compulsory and nesting badly when a user wrapped QAP() in their
own. It now only emits a progressor, as progressr intends; users switch
reporting on with progressr::with_progress(QAP(...)). Signals are batched
at ~1% of the run rather than one condition throw per permutation.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
predict() reshapes back to the network. type = "matrix" returns the
predictions in the shape of the input with NA in every structurally
excluded cell, which is what you want to plot a predicted network. The
valid-cell rules are factored out of qap_ols_template() into
qap_valid_mask() and shared, rather than restated.

plot() draws the permutation distribution per coefficient with the
observed statistic marked -- the picture the p-value summarises. Free, now
that null_dist is retained.

tidy()/glance() register with broom at load time if broom is installed, so
it stays out of Imports. std.error is deliberately absent: a QAP fit has
no single standard error that the permutation p-value corresponds to.

nobs() no longer returns NA under less_mem: the count is recorded when the
data are vectorised. Its dead null_dist line is gone.

weights and offset are new arguments, supported for every family and for
BOTH estimators. They travel as reserved columns of the vectorised frame,
so they are masked, stacked and NA-dropped exactly like a predictor but
never enter the model matrix, get residualised, or get permuted -- a
weight and an offset are properties of the cell, not of the actors, and
that is documented in ?QAP rather than left implicit.

For GMM they enter the moment conditions themselves: gmm_w()/gmm_off()
default to 1 and 0 so an unweighted fit is bit-identical, and weights
scale each observation's contribution while the offset joins the linear
predictor. Verified against the equivalent GLM -- for poisson and logit
the moment conditions are the score equations and the system is
just-identified, so GMM lands on the weighted, offset MLE to 0.00e+00.

Weighted residualisation uses the same weights for the semi-partialling;
the offset does not, having no place in a predictor-on-predictor
regression. Weighted fits take the general path: qap_ols_eligible()
returns FALSE, since the FWL decomposition and the closed-form HC3 are
both unweighted. The degeneracy guard in QAPPermEst() now skips the
reserved columns, or a constant weight vector would reject every
permutation.

Adds test-weights-offset.R (43 assertions) and test-fwl.R.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Input adapters (qap_coerce_data), run before anything else reads `data`:

- a dyadic or triadic data frame, handed to df_to_mat() with sender=,
  receiver=, perceiver= and split_by= -- the function that already does
  exactly this, rather than a second implementation
- igraph and network objects anywhere in the list, converted to adjacency
  matrices, reading edge_attr= when given

Both are checked against the equivalent matrix input for an identical fit,
because an adapter that quietly changes the model is worse than none. The
class check has to precede the list check: igraph and network objects are
themselves lists, so testing for a list of networks first iterates their
internal structure.

GPU: torch and CUDA turn out to be available on this machine, so the torch
algebra is now exercised rather than assumed. test-gpu-torch.R runs
gpu_solve_fixed_x() and gpu_solve_varying_x() against qap_ols_solve() on
both the CPU and CUDA backends, checks CUDA and CPU torch against each
other, covers partial final batches, multi-network, grouped and CSS data
through QAP(use_gpu = TRUE), and confirms the ineligible cases warn and
fall back. It skips cleanly where torch is absent.

The GPU draws permutations in batches and the CPU one at a time, so the
two see different permutations; the end-to-end test holds the baseline fit
to 1e-10 and the p-values to Monte Carlo error, which is the honest
comparison rather than a tolerance chosen to pass.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
less_mem promised not to store the fitted model but skipped only the
top-level copy: the same object survived in fit$base$base_model. At n = 120
the "smaller" result was 4.68 MB of 4.69 MB, while fitted(), vcov(),
logLik() and residuals() all refused to run and advised refitting without
it. Dropping it for real takes the fit from 9.37 MB to 0.01 MB. coef(),
nobs() and summary() still work, and that is now what ?QAP says.

This came out of measuring the model frame, which the plan gated on
evidence rather than assumption: it is 0.96 MB of a 12.64 MB lm at n = 200,
so model = FALSE was not worth its risk to predict() and update(). The
real waste was elsewhere.

Infrastructure:
- R-CMD-check on Windows, macOS and Linux across release, devel and
  oldrel-1; test-coverage via covr
- a pkgdown site with a grouped reference index
- NEWS.md and the vignette cover predict(), plot(), weights and offsets and
  their permutation semantics, the new input shapes, opt-in progress,
  pooling, FWL, and where the GPU is and is not worth using

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
predict.lm() has no "link" type, so asking a gaussian fit for one errored
rather than returning the identity-link values it obviously means. Mapped
to the response scale, which for a linear model is the same thing.

Adds test-predict-plot.R: the shapes predict() restores for directed,
undirected, CSS and multi-network fits, link vs response for a GLM, plot()
including its refusals, and the broom methods.

README and ?QAP's use_gpu entry now describe what the GPU path actually
does and where it is worth using, verified on CUDA rather than assumed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The structural mask says which cells the model COULD use -- not which it
did. Any NA in the outcome, a predictor, the weights or the offset drops
its cell too, so the fitted vector is shorter than the mask has slots.
Filling the first n slots put every value after the first gap in the wrong
cell: on an 8-node network with two missing entries, the two cells that
should have been NA held numbers and two others were blanked instead.

Fits now record which structurally-valid cells survived (fit$kept, absent
when none were lost), and the reshaper places by that. newdata is
unaffected -- it carries its own missingness, not the fit's.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
torch installs in two halves -- the R package, and the LibTorch binaries
it downloads on first use -- and every guard here checked only the first.
On a machine with the package but no back end, requireNamespace("torch")
returns TRUE, so gpu_batch_ols() and the whole torch test suite walked
straight into "Lantern is not loaded".

That is the ordinary state of a CI runner: `needs: check` installs
Suggests, torch among them, and nothing fetches LibTorch. Both workflows
carried a comment claiming torch was not installed at all, which is what
hid this.

torch_ready() checks both halves. gpu_available() reports FALSE there
instead of erroring, use_gpu = TRUE says which half is missing, and the
torch tests skip rather than fail. Verified against a stub torch package
standing in for a runner: 15 failures before, 0 failures and 20 skips
after, with the full suite still green where torch really works.

Nothing touches the permutation draws, so the golden-master fixtures are
unaffected.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@RWKrause
RWKrause merged commit 7e46837 into main Sep 8, 2026
12 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant