From bc96783a3ed6c427aa216c86be622da8a81bb257 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:19:22 +0200 Subject: [PATCH 01/12] Renamed the engine to the front end's vocabulary, and fixed #4 Fixed `net_regression()` failing on a two-mode network with more columns than rows (closing #4). The validity mask was built as rows-by-rows, so a wider predictor extended it with NA and the dyad count came back as NA. Renamed `reps` to `times`, `nullhyp`/`method` to `permute` with the values `"predictor"` and `"outcome"`, and `mode` to a logical `directed`. Retired `data` as an identifier: `matlist` is the list of matrices the engine fits and `net` is one coerced network in the formula front end. Added reporting of every default the model resolves for itself, and test-qap_reporting.R, which runs with snet_verbosity = "verbose" so that a message cli cannot parse cannot stay invisible. Co-Authored-By: Claude Opus 5 --- .github/CONTRIBUTING.md | 53 ++++++- DESCRIPTION | 2 +- NEWS.md | 45 ++++++ R/model_regression.R | 197 +++++++++++++++------------ R/model_tests.R | 12 +- R/qap_css.R | 169 +++++++++++------------ R/qap_engine.R | 109 ++++++++------- R/qap_gpu.R | 60 ++++---- R/qap_misc.R | 61 +++++---- R/qap_utils.R | 26 ++-- man/regression.Rd | 16 ++- tests/testthat/helper-infernet.R | 22 ++- tests/testthat/test-model_tests.R | 6 +- tests/testthat/test-net_regression.R | 6 +- tests/testthat/test-qap_control.R | 33 +++-- tests/testthat/test-qap_reporting.R | 74 ++++++++++ tests/testthat/test-qap_shapes.R | 26 +++- 17 files changed, 581 insertions(+), 336 deletions(-) create mode 100644 tests/testthat/test-qap_reporting.R diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index 859486f..23320c1 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -104,6 +104,28 @@ as well as the burden on users to understand all of the options. Use sensible defaults instead. Function and argument names should also follow the house rules (see below). +One word means one thing, on both sides of the seam between the formula front +end and the engine. The engine was ported from `MrQAP` and used its own +vocabulary; the front end's words won, since those are the ones users read: + +| Word | Means | Not | +|---|---|---| +| `times` | how many permutations | `reps` | +| `directed` | logical, whether i→j differs from j→i | `mode`, `"digraph"`/`"graph"` | +| `permute` | what the null distribution permutes: `"predictor"` or `"outcome"` | `nullhyp`, `method`, `"qapspp"`/`"qapy"` | +| `.data` | the network the user passes in | — | +| `matlist` | the named list of matrices the engine fits | `data` | +| `net` | one coerced network, inside the formula front end | `data` | + +`mode` is reserved for a nodeset, as in one-mode and two-mode, which is what it +means everywhere else in the ecosystem. Do not use it for directedness. +`permute` replaced `method` because "method" says nothing about what differs; +`"predictor"` and `"outcome"` name the thing that is actually shuffled. +`data` is retired as an identifier: it named the network in one half of +[R/model_regression.R](../R/model_regression.R) and the matrix list in the +other, one letter away from `.data`. Reserve `data =` for the argument a model +fitter takes. + When writing documentation or NEWS items, prefer breaking lines at punctuation. Make it clear when you are referring to functions by adding backticks and parentheses, @@ -212,17 +234,21 @@ regression entry point. Its control flow is: network, drops the ones that are missing a predictor with a warning, and pools the rest. 3. Resolve `family = "auto"` against the dependent variable - (binomial for a 0/1 outcome, gaussian otherwise), and resolve `mode` and + (binomial for a 0/1 outcome, gaussian otherwise), and resolve `directed` and `diag` from the network with `manynet::is_directed()` and `manynet::is_complex()`. + Report each resolution with `snet_info()`: a model the user did not state is + one they cannot describe in a paper. 4. Call `QAPglm()`, which parses the formula, fits the baseline model once - via `fit_qap_model()`, then runs `reps` permutations and aggregates them. + via `fit_qap_model()`, then runs `times` permutations and aggregates them. 5. Attach a probabilistic confusion matrix where the outcome is binary, and class the result `net_regression`. -Inside `QAPglm()` the null hypothesis decides the permutation scheme: -`"qapy"` permutes the dependent matrix only, while `"qapspp"` implements Dekker -et al.'s double semi-partialling, running one permutation set per main predictor -after residualising it against the others. +Inside `QAPglm()` the `permute` control names what the null distribution +permutes: `"outcome"` permutes the dependent matrix only, while `"predictor"` +implements Dekker et al.'s double semi-partialling, running one permutation set +per main predictor after residualising it against the others. +With one predictor there is nothing to residualise against, so `"predictor"` +falls back to `"outcome"` and says so. Permuted coefficients and test statistics are then compared against the baseline by `compare_perm_to_baseline()` and reduced to `lower`/`larger`/`abs` p-value matrices by `aggregate_perm_results()`. @@ -328,6 +354,12 @@ Users opt in with e.g. `options(snet_verbosity = "verbose")`. These wrappers pass their input to `{cli}`, so: - Braces interpolate, replacing `paste()`: `snet_abort("{.val {dep}} is not in the data.")`. +- A brace expression beginning with a dot is read as a *style*, not as code, so + `{.val {.directed_label(x)}}` aborts with "Invalid cli literal". Resolve a + call to a dot-prefixed function into a local variable first. +- `snet_info()` pastes its arguments, so pass separate strings for a longer + message rather than a named `c()` vector: the names are dropped and the + strings run together without a space. - Use `{cli}` inline classes to mark up what you refer to — `{.fn}` for functions, `{.arg}`/`{.var}` for arguments and variables, `{.val}` for values, `{.pkg}` for packages, `{.url}` for links. @@ -335,6 +367,15 @@ These wrappers pass their input to `{cli}`, so: `snet_warn("Dropped {length(dropped)} network{?s}.")`. Prefer "`{.arg times}` must be a positive whole number" over "invalid input". + +Report every default the model resolves for itself, with `snet_info()`: the +family read from the outcome's values, the directedness read from the network, +and any fallback such as `permute = "predictor"` reducing to `"outcome"`. +A model the user did not state is one they cannot describe in a paper. +Because this output is silent by default, a broken message is invisible in +every other test, so cover it in +[tests/testthat/test-qap_reporting.R](../tests/testthat/test-qap_reporting.R), +which runs with `snet_verbosity = "verbose"`. Where a function needs a package from `Suggests`, name it and say how to get it: `snet_abort(c("The {.pkg lme4} package is required for random effects.", i = "Install it with {.run install.packages(\"lme4\")}."))`. diff --git a/DESCRIPTION b/DESCRIPTION index d9a01f2..8740e49 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: infernet Title: Inferential Models for Many Different Types of Networks -Version: 0.1.1 +Version: 0.2.0 Description: A set of tools for testing networks. It includes functions for univariate and multivariate conditional uniform graph and quadratic assignment procedure testing, diff --git a/NEWS.md b/NEWS.md index 61cc9c9..b914f08 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,48 @@ +# infernet 0.2.0 + +## Package + +- Renamed the engine's vocabulary to the front end's, so one word means one + thing on both sides of the seam + - `reps` is now `times`, everywhere including on the returned fit + - `mode` is now `directed`, a logical, and `"digraph"`/`"graph"` are gone; + `mode` is reserved for a nodeset, as in one-mode and two-mode + - `nullhyp` is now `permute`, and its values name what is shuffled: + `"predictor"` for Dekker's double semi-partialling, `"outcome"` for + permuting the dependent matrix alone + - `data` is retired as an identifier: it named the network in one half of + `R/model_regression.R` and the matrix list in the other, one letter away + from `.data` + - `matlist` is the named list of matrices the engine fits + - `net` is one coerced network, inside the formula front end + - `.data` remains the network the user passes in +- Updated CONTRIBUTING with the vocabulary table and the reporting rule + +## Regression + +- Fixed `net_regression()` failing on a two-mode network with more columns than + rows (closing #4) + - The validity mask was built as rows-by-rows, so a wider predictor extended + it with `NA` and the dyad count came back as `NA` + - The reported 448 by 12489 network now fits, on all 5,595,072 dyads +- Renamed the `method` control to `permute` + - `method = "qap"` is now `permute = "predictor"`, and `method = "qapy"` is + now `permute = "outcome"` +- Renamed the `mode` control to `directed` + - `mode = "undirected"` is now `directed = FALSE` +- Added reporting of every default the model resolves for itself + - The family chosen from the outcome's values + - The directedness read from the network + - `permute = "predictor"` falling back to `"outcome"` with one predictor + - These use `snet_info()`, so `options(snet_verbosity = "verbose")` shows them + +## Tests + +- Added a wide two-mode fixture and two regression tests for #4 +- Added `test-qap_reporting.R`, which runs with `snet_verbosity = "verbose"` + - Informational output is silent in every other test, so a message that + `{cli}` cannot parse was invisible until it aborted; two shipped that way + # infernet 0.1.1 ## Package diff --git a/R/model_regression.R b/R/model_regression.R index 5b639d9..023457f 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -12,8 +12,8 @@ #' #' - gaussian, binomial, poisson, negbin, zero-inflated Poisson, and #' multinomial families; -#' - `"qap"` (Dekker's double semi-partialling plus) and `"qapy"` -#' (permute-y-only) null hypotheses; +#' - two permutation schemes: `"predictor"` (Dekker's double semi-partialling) +#' and `"outcome"`; #' - random intercepts (lme4 / glmmTMB) and fixed effects (fixest); #' - robust (HC3) standard errors; #' - optional torch-based batch OLS on the GPU; @@ -40,16 +40,20 @@ #' 1000 is the default; publication-ready work usually needs 1000-10000. #' @param control Named list of additional controls; unspecified entries fall #' back to the defaults below. -#' - `method`: `"qap"` (double semi-partialling plus, default) or `"qapy"` -#' (permute y only). +#' - `permute`: what the null distribution permutes. `"predictor"` (the +#' default) residualises each main predictor against the others and +#' permutes that residual, following Dekker et al. (2007). `"outcome"` +#' permutes the dependent matrix and leaves the predictors alone. With one +#' predictor there is nothing to residualise against, so `"predictor"` +#' reduces to `"outcome"` and says so. #' - `strategy`: future plan, e.g. `"sequential"` (default), `"multisession"`. #' - `family`: `"auto"` (default; gaussian for weighted networks, binomial #' for binary), `"gaussian"`, `"binomial"`, `"poisson"`, `"negbin"`, #' `"zip"`, or `"multinom"`. #' - `estimator`: `"standard"` (default) or `"gmm"` (binomial/poisson/ #' negbin/zip). -#' - `mode`: `"directed"` / `"undirected"` (default auto-detected from -#' `.data`). +#' - `directed`: logical, whether a tie from i to j differs from one from j +#' to i. Read from `.data` unless given, and reported when read. #' - `diag`: logical, include loops (default auto-detected). #' - `seed`, `groups`, `ncores`: passed through to the engine. #' - `use_robust_errors`: HC3 standard errors. @@ -97,42 +101,49 @@ net_regression <- function(formula, if (.is_list_of_graphs(.data)) { prepared <- .prepare_list_of_graphs(formula, .data) - data <- prepared$data + matlist <- prepared$matlist formula <- prepared$formula first_graph <- prepared$first_graph } else { ml <- convertToMatrixList(formula, .data) - data <- ml$mydata + matlist <- ml$mydata formula <- ml$formula first_graph <- manynet::as_tidygraph(.data) } dep <- .dep_name(formula) + # Each of these is resolved from the data rather than stated by the user, so + # each is reported. A defaulted family or a fallback that nobody sees is a + # model the user cannot describe in a paper. if (identical(ctrl$family, "auto")) { - ctrl$family <- if (.is_binary_outcome(data[[dep]])) "binomial" else "gaussian" + ctrl$family <- if (.is_binary_outcome(matlist[[dep]])) "binomial" else "gaussian" + manynet::snet_info( + "Treating the outcome as {.val {ctrl$family}}, from its values.") } user_requested_gaussian_binary <- - identical(ctrl$family, "gaussian") && .is_binary_outcome(data[[dep]]) + identical(ctrl$family, "gaussian") && .is_binary_outcome(matlist[[dep]]) - if (is.null(ctrl$mode)) { - ctrl$mode <- if (manynet::is_directed(first_graph)) "directed" else "undirected" + if (is.null(ctrl$directed)) { + ctrl$directed <- manynet::is_directed(first_graph) + # `{cli}` reads a brace expression beginning with a dot as a style, so a + # call to a dot-prefixed function has to be resolved before interpolation. + direction <- .directed_label(ctrl$directed) + manynet::snet_info("Reading the network as {.val {direction}}.") } if (is.null(ctrl$diag)) { ctrl$diag <- isTRUE(manynet::is_complex(first_graph)) } - nullhyp <- if (ctrl$method == "qap") "qapspp" else "qapy" - fit <- QAPglm( formula = formula, - data = data, + matlist = matlist, family = ctrl$family, - mode = ctrl$mode, + directed = ctrl$directed, diag = ctrl$diag, - nullhyp = nullhyp, + permute = ctrl$permute, estimator = ctrl$estimator, - reps = times, + times = times, seed = ctrl$seed, groups = ctrl$groups, strategy = ctrl$strategy, @@ -167,7 +178,7 @@ net_regression <- function(formula, .resolve_control <- function(control = list()) { ctrl <- .default_control() if (length(control) == 0) { - ctrl$method <- match.arg(ctrl$method, choices = c("qap", "qapy")) + ctrl$permute <- match.arg(ctrl$permute, choices = .permute_schemes()) return(ctrl) } if (is.null(names(control)) || any(!nzchar(names(control)))) { @@ -187,17 +198,33 @@ net_regression <- function(formula, manynet::snet_abort(msg) } ctrl[names(control)] <- control - ctrl$method <- match.arg(ctrl$method, choices = c("qap", "qapy")) + ctrl$permute <- match.arg(ctrl$permute, choices = .permute_schemes()) ctrl } +# What the permutation scheme permutes, which is the only thing that separates +# the two. "predictor" is Dekker et al's double semi-partialling: each main +# predictor is residualised against the others and that residual is permuted. +# "outcome" permutes the dependent matrix and leaves the predictors alone. +#' @keywords internal +#' @noRd +.permute_schemes <- function() c("predictor", "outcome") + +#' @keywords internal +#' @noRd +.permute_label <- function(permute) { + switch(permute, + predictor = "each predictor's residuals (Dekker's double semi-partialling)", + outcome = "the outcome only") +} + .default_control <- function() { list( - method = c("qap", "qapy"), + permute = .permute_schemes(), strategy = "sequential", family = "auto", estimator = "standard", - mode = NULL, + directed = NULL, diag = NULL, seed = NULL, groups = NULL, @@ -280,17 +307,17 @@ net_regression <- function(formula, kept_mls <- ml_list[keep] ref_names <- names(kept_mls[[1]]$mydata) - data <- vector("list", length(ref_names)) - names(data) <- ref_names + matlist <- vector("list", length(ref_names)) + names(matlist) <- ref_names for (nm in ref_names) { - data[[nm]] <- lapply(kept_mls, function(ml) ml$mydata[[nm]]) + matlist[[nm]] <- lapply(kept_mls, function(ml) ml$mydata[[nm]]) } first_graph <- manynet::as_tidygraph(glist[[keep[1]]]) specificationAdvice(getRHSNames(formula)$IVnames, first_graph) list( - data = data, + matlist = matlist, formula = kept_mls[[1]]$formula, first_graph = first_graph ) @@ -363,12 +390,12 @@ print.net_regression <- function(x, ..., if (!is.null(x$groups)) cat("Permutations were performed within groups only.\n") - if (x$nullhyp == "qapy") - cat("The outcome matrix Y was permuted", format(x$reps), "times.\n") - if (x$nullhyp == "qapspp") { + if (x$permute == "outcome") + cat("The outcome matrix Y was permuted", format(x$times), "times.\n") + if (x$permute == "predictor") { cat("Significance was estimated using Dekker's\n") cat(" 'semi-partialling plus' procedure with", - format(x$reps), "permutations.\n") + format(x$times), "permutations.\n") } if (x$diag) { @@ -376,7 +403,8 @@ print.net_regression <- function(x, ..., } else { cat("Diagonal values (loops) were ignored.\n") } - cat("The outcome was treated as", format(paste0(x$mode, ".")), "\n") + cat("The outcome was treated as", + format(paste0(.directed_label(x$directed), ".")), "\n") if (!is.null(x$r.squared)) { cat("\nR-squared: ", format(round(x$r.squared, 4))) @@ -391,7 +419,7 @@ print.net_regression <- function(x, ..., format(x$abs[1, ])) colnames(cmat) <- c("Estimate", "Pr(<=b)", "Pr(>=b)", "Pr(>=|b|)") rownames(cmat) <- names(x$coefficients) - if (x$nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (x$permute == "predictor") cmat[1, 2:4] <- "*" print.table(cmat) cat("\n--------------\n") } @@ -402,11 +430,11 @@ print.net_regression <- function(x, ..., format(x$abs[2, ])) colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- names(x$coefficients) - if (x$nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (x$permute == "predictor") cmat[1, 2:4] <- "*" print.table(cmat) - if (x$nullhyp == "qapspp") - cat("\n* Significance test for the intercept is not defined with qapspp.\n") + if (x$permute == "predictor") + cat("\n* The intercept has no significance test when predictors are permuted.\n") cat("\n--------------\n") @@ -445,12 +473,12 @@ print.net_regression <- function(x, ..., if (!is.null(x$groups)) cat("\nPermutations were performed within groups only.") - if (x$nullhyp == "qapy") - cat("\nThe outcome matrix Y was permuted", format(x$reps), "times.") - if (x$nullhyp == "qapspp") { + if (x$permute == "outcome") + cat("\nThe outcome matrix Y was permuted", format(x$times), "times.") + if (x$permute == "predictor") { cat("\nSignificance was estimated using Dekker's") cat("\n 'semi-partialling plus' procedure with", - format(x$reps), "permutations.") + format(x$times), "permutations.") } if (x$diag) { @@ -458,7 +486,8 @@ print.net_regression <- function(x, ..., } else { cat("\nDiagonal values (loops) were ignored.") } - cat("\nThe outcome was treated as", format(paste0(x$mode, "."))) + cat("\nThe outcome was treated as", + format(paste0(.directed_label(x$directed), "."))) cat("\nModel family:", format(x$family)) if (!is.null(x$comp)) { @@ -466,7 +495,7 @@ print.net_regression <- function(x, ..., cat("\n\n--- Comparison:", names(x$comp)[k], "---") cat("\n ", x$comp[[k]][1], "vs", x$comp[[k]][2]) .print_glm_table(x$base[[k]], x$lower[[k]], x$larger[[k]], x$abs[[k]], - x$nullhyp, print_b) + x$permute, print_b) } } else { cat("\n\nCoefficients:\n") @@ -477,7 +506,7 @@ print.net_regression <- function(x, ..., cmat[, 3] <- format(x$lower[1, ]) cmat[, 4] <- format(x$larger[1, ]) cmat[, 5] <- format(x$abs[1, ]) - if (x$nullhyp == "qapspp") cmat[1, 3:5] <- "*" + if (x$permute == "predictor") cmat[1, 3:5] <- "*" colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=b)", "Pr(>=b)", "Pr(>=|b|)") rownames(cmat) <- names(x$coefficients) print.table(cmat) @@ -490,14 +519,14 @@ print.net_regression <- function(x, ..., cmat[, 3] <- format(x$lower[2, ]) cmat[, 4] <- format(x$larger[2, ]) cmat[, 5] <- format(x$abs[2, ]) - if (x$nullhyp == "qapspp") cmat[1, 3:5] <- "*" + if (x$permute == "predictor") cmat[1, 3:5] <- "*" colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- names(x$coefficients) print.table(cmat) } - if (x$nullhyp == "qapspp") - cat("\n* Significance test for the intercept is not defined with qapspp.\n") + if (x$permute == "predictor") + cat("\n* The intercept has no significance test when predictors are permuted.\n") cat("--------------\n") @@ -523,7 +552,7 @@ print.net_regression <- function(x, ..., } -.print_glm_table <- function(base, lower, larger, abs_mat, nullhyp, print_b) { +.print_glm_table <- function(base, lower, larger, abs_mat, permute, print_b) { cat("\n\nCoefficients:\n") nc <- length(base$coefficients) cmat <- matrix(NA, nrow = nc, ncol = 4) @@ -531,7 +560,7 @@ print.net_regression <- function(x, ..., cmat[, 2] <- format(lower[2, ]) cmat[, 3] <- format(larger[2, ]) cmat[, 4] <- format(abs_mat[2, ]) - if (nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (permute == "predictor") cmat[1, 2:4] <- "*" colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- names(base$coefficients) print.table(cmat) @@ -545,8 +574,8 @@ print.net_regression <- function(x, ..., #' @keywords internal #' @noRd convertToMatrixList <- function(formula, .data, advise = TRUE) { - data <- manynet::as_tidygraph(.data) - DV <- manynet::as_matrix(data) + net <- manynet::as_tidygraph(.data) + DV <- manynet::as_matrix(net) # The sender and the receiver of a tie come from one nodeset in a one-mode # network and from two in a two-mode one, so a dyadic term must read the # attribute once per mode rather than once per network. @@ -560,52 +589,52 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { } list(rows = rows, cols = cols) } - twomode <- manynet::is_twomode(data) - node_type <- if (twomode) manynet::node_attribute(data, "type") else NULL + twomode <- manynet::is_twomode(net) + node_type <- if (twomode) manynet::node_attribute(net, "type") else NULL names_form <- getRHSNames(formula) - .check_formula_vars(names_form$IVnames, data) - if (advise) specificationAdvice(names_form$IVnames, data) + .check_formula_vars(names_form$IVnames, net) + if (advise) specificationAdvice(names_form$IVnames, net) IVs <- lapply(names_form$IVnames, function(IV) { out <- lapply(seq_along(IV), function(elem) { if (IV[[elem]][1] == "ego") { - vct <- manynet::node_attribute(data, IV[[elem]][2]) - if (manynet::is_twomode(data)) { - vct <- vct[!manynet::node_attribute(data, "type")] + vct <- manynet::node_attribute(net, IV[[elem]][2]) + if (manynet::is_twomode(net)) { + vct <- vct[!manynet::node_attribute(net, "type")] } out <- matrix(vct, nrow(DV), ncol(DV)) out <- list(out) names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "alter") { - vct <- manynet::node_attribute(data, IV[[elem]][2]) - if (manynet::is_twomode(data)) { - vct <- vct[manynet::node_attribute(data, "type")] + vct <- manynet::node_attribute(net, IV[[elem]][2]) + if (manynet::is_twomode(net)) { + vct <- vct[manynet::node_attribute(net, "type")] } out <- matrix(vct, nrow(DV), ncol(DV), byrow = TRUE) out <- list(out) names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "same") { - attrib <- manynet::node_attribute(data, IV[[elem]][2]) + attrib <- manynet::node_attribute(net, IV[[elem]][2]) if (manynet::is_twomode(.data)) { if (all(is.na(attrib[!manynet::node_is_mode(.data)]))) { attrib <- attrib[manynet::node_is_mode(.data)] out <- vapply(1:length(attrib), function(x) { - net <- manynet::as_matrix( + held_out <- manynet::as_matrix( manynet::delete_nodes(.data, manynet::net_dims(.data)[1] + x)) - rowSums(net * matrix((attrib[-x] == attrib[x]) * 1, - nrow(DV), ncol(DV) - 1, byrow = TRUE)) / - rowSums(net) + rowSums(held_out * matrix((attrib[-x] == attrib[x]) * 1, + nrow(DV), ncol(DV) - 1, byrow = TRUE)) / + rowSums(held_out) }, FUN.VALUE = numeric(nrow(DV))) out[is.nan(out)] <- 0 } else { attrib <- attrib[!manynet::node_is_mode(.data)] out <- t(vapply(1:length(attrib), function(x) { - net <- manynet::as_matrix(manynet::delete_nodes(.data, x)) - colSums(net * matrix((attrib[-x] == attrib[x]) * 1, - nrow(DV) - 1, ncol(DV))) / - colSums(net) + held_out <- manynet::as_matrix(manynet::delete_nodes(.data, x)) + colSums(held_out * matrix((attrib[-x] == attrib[x]) * 1, + nrow(DV) - 1, ncol(DV))) / + colSums(held_out) }, FUN.VALUE = numeric(ncol(DV)))) out[is.nan(out)] <- 0 } @@ -618,12 +647,12 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "dist") { - if (is.character(manynet::node_attribute(data, IV[[elem]][2]))) { + if (is.character(manynet::node_attribute(net, IV[[elem]][2]))) { manynet::snet_abort( c("{.fn dist} is undefined for a categorical attribute.", i = "Try {.fn same} instead.")) } - sides <- side_matrices(manynet::node_attribute(data, IV[[elem]][2]), + sides <- side_matrices(manynet::node_attribute(net, IV[[elem]][2]), DV, twomode, node_type) rows <- sides$rows cols <- sides$cols @@ -632,12 +661,12 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "sim") { - if (is.character(manynet::node_attribute(data, IV[[elem]][2]))) { + if (is.character(manynet::node_attribute(net, IV[[elem]][2]))) { manynet::snet_abort( c("{.fn sim} is undefined for a categorical attribute.", i = "Try {.fn same} instead.")) } - sides <- side_matrices(manynet::node_attribute(data, IV[[elem]][2]), + sides <- side_matrices(manynet::node_attribute(net, IV[[elem]][2]), DV, twomode, node_type) rows <- sides$rows cols <- sides$cols @@ -648,9 +677,9 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "tertius") { - vct <- manynet::node_attribute(data, IV[[elem]][2]) - if (manynet::is_twomode(data)) { - vct <- vct[!manynet::node_attribute(data, "type")] + vct <- manynet::node_attribute(net, IV[[elem]][2]) + if (manynet::is_twomode(net)) { + vct <- vct[!manynet::node_attribute(net, "type")] } val <- matrix(vct, nrow(DV), ncol(DV)) * DV if (is.na(IV[[elem]][3])) { @@ -675,8 +704,8 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { names(out) <- paste(IV[[elem]][1:2], collapse = " ") out } else { - if (IV[[elem]][1] %in% manynet::net_tie_attributes(data)) { - out <- manynet::as_matrix(manynet::to_uniplex(data, + if (IV[[elem]][1] %in% manynet::net_tie_attributes(net)) { + out <- manynet::as_matrix(manynet::to_uniplex(net, tie = IV[[elem]][1])) out <- list(out) names(out) <- IV[[elem]][1] @@ -791,14 +820,14 @@ getRHSNames <- function(formula) { #' @keywords internal #' @noRd -.check_formula_vars <- function(IVnames, data) { +.check_formula_vars <- function(IVnames, net) { node_fns <- c("ego", "alter", "same", "dist", "sim", "tertius") - node_attrs <- manynet::net_node_attributes(data) + node_attrs <- manynet::net_node_attributes(net) # The engine builds these columns itself in `make_qap_data()`: the sender, the # receiver, the network, and the perceiver index. They are what a user names # after a `|` to absorb sender or receiver fixed effects, so they are not # attributes of the network and must not be looked for among them. - tie_attrs <- c(manynet::net_tie_attributes(data), .structural_vars()) + tie_attrs <- c(manynet::net_tie_attributes(net), .structural_vars()) missing_node <- character(0) missing_tie <- character(0) @@ -823,7 +852,7 @@ getRHSNames <- function(formula) { i = "Available node attributes: {.val {node_attrs}}.")) } if (length(missing_tie) > 0) { - available <- manynet::net_tie_attributes(data) + available <- manynet::net_tie_attributes(net) structurals <- .structural_vars() manynet::snet_abort( c("Tie attribute or predictor{?s} {.val {unique(missing_tie)}} not found.", @@ -849,13 +878,13 @@ getDependentName <- function(formula) { #' @keywords internal #' @noRd -specificationAdvice <- function(formula, data) { +specificationAdvice <- function(formula, net) { formdf <- t(data.frame(formula)) if (any(formdf[, 1] %in% c("sim", "same"))) { vars <- formdf[formdf[, 1] %in% c("sim", "same"), 2] suggests <- vapply(vars, function(x) { incl <- unname(formdf[formdf[, 2] == x, 1]) - if (manynet::is_twomode(data)) { + if (manynet::is_twomode(net)) { excl <- setdiff(c("ego", "tertius"), incl) } else excl <- setdiff(c("ego", "alter"), incl) if (length(excl) > 0) { @@ -865,7 +894,7 @@ specificationAdvice <- function(formula, data) { } }, FUN.VALUE = character(1)) suggests <- suggests[!is.na(suggests)] - if (!manynet::is_directed(data)) { + if (!manynet::is_directed(net)) { suggests <- suggests[!grepl("ego\\(", suggests)] } if (length(suggests) > 0) { diff --git a/R/model_tests.R b/R/model_tests.R index 908b110..2738c1a 100644 --- a/R/model_tests.R +++ b/R/model_tests.R @@ -70,12 +70,12 @@ test_random <- function(.data, FUN, ..., out <- list(test = "CUG", testval = obsd, testdist = simd, - mode = manynet::is_directed(.data), + directed = manynet::is_directed(.data), diag = manynet::is_complex(.data), cmode = "edges", plteobs = mean(simd <= obsd), pgteobs = mean(simd >= obsd), - reps = times) + times = times) class(out) <- "network_test" out } @@ -115,12 +115,12 @@ test_configuration <- function(.data, FUN, ..., out <- list(test = "configuration", testval = obsd, testdist = simd, - mode = manynet::is_directed(.data), + directed = manynet::is_directed(.data), diag = manynet::is_complex(.data), cmode = "edges", plteobs = mean(simd <= obsd), pgteobs = mean(simd >= obsd), - reps = times) + times = times) class(out) <- "network_test" out } @@ -162,11 +162,11 @@ test_permutation <- function(.data, FUN, ..., out <- list(test = "QAP", testval = obsd, testdist = simd, - mode = manynet::is_directed(.data), + directed = manynet::is_directed(.data), diag = manynet::is_complex(.data), plteobs = mean(simd <= obsd), pgteobs = mean(simd >= obsd), - reps = times) + times = times) class(out) <- "network_test" out } diff --git a/R/qap_css.R b/R/qap_css.R index 7b14fe4..c003bee 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -6,10 +6,10 @@ #' @keywords internal #' @noRd -array_to_vector <- function(ar, mode., diag.) { +array_to_vector <- function(ar, directed., diag.) { v <- c() for (i in 1:nrow(ar)) { - if (mode. == 'undirected') { + if (!directed.) { v <- c(v, as.vector(ar[, , i][upper.tri(ar[, , i], diag = diag.)])) } else { v <- c(v, as.vector(ar[, , i])) @@ -21,7 +21,7 @@ array_to_vector <- function(ar, mode., diag.) { #' @keywords internal #' @noRd -make_css_data <- function(y, x, nets, diag, mode) { +make_css_data <- function(y, x, nets, diag, directed) { n <- dim(y)[1] nx <- length(x) valid <- array(TRUE, dim = c(n, n, n)) @@ -41,7 +41,7 @@ make_css_data <- function(y, x, nets, diag, mode) { valid[is.na(x[[var]])] <- FALSE } - if (mode == 'undirected') { + if (!directed) { for (i in 1:n) { y[, , i][lower.tri(y[, , i])] <- NA valid[, , i][lower.tri(valid[, , i])] <- FALSE @@ -57,8 +57,8 @@ make_css_data <- function(y, x, nets, diag, mode) { x[[var]][!valid] <- NA } - vv <- array_to_vector(valid, mode. = mode, diag. = diag) - yv <- array_to_vector(y, mode. = mode, diag. = diag)[vv] + vv <- array_to_vector(valid, directed. = directed, diag. = diag) + yv <- array_to_vector(y, directed. = directed, diag. = diag)[vv] pred <- data.frame(yv = yv, nv = nets) @@ -70,13 +70,13 @@ make_css_data <- function(y, x, nets, diag, mode) { per[, , i] <- i } - pred$sv <- as.factor(array_to_vector(sen, mode. = mode, diag. = diag)[vv]) - pred$rv <- as.factor(array_to_vector(rec, mode. = mode, diag. = diag)[vv]) - pred$pv <- as.factor(array_to_vector(per, mode. = mode, diag. = diag)[vv]) + pred$sv <- as.factor(array_to_vector(sen, directed. = directed, diag. = diag)[vv]) + pred$rv <- as.factor(array_to_vector(rec, directed. = directed, diag. = diag)[vv]) + pred$pv <- as.factor(array_to_vector(per, directed. = directed, diag. = diag)[vv]) for (var in c(1:nx)) { pred[[names(x)[var]]] <- array_to_vector(x[[var]], - mode. = mode, diag. = diag)[vv] + directed. = directed, diag. = diag)[vv] } return(list(pred = pred, valid = valid)) } @@ -85,9 +85,9 @@ make_css_data <- function(y, x, nets, diag, mode) { #' @keywords internal #' @noRd QAPcssPermEst <- function(i, - data., + matlist., perm_var., - mode., + directed., diag., mod., groups., @@ -105,9 +105,9 @@ QAPcssPermEst <- function(i, reference.) { dep <- parsed.$dependent - large <- is.list(data.[[dep]]) + large <- is.list(matlist.[[dep]]) - y_cat <- stats::na.omit(unique(as.vector(unlist(data.[[dep]])))) + y_cat <- stats::na.omit(unique(as.vector(unlist(matlist.[[dep]])))) sufficient_data <- FALSE trial <- 0 @@ -116,7 +116,7 @@ QAPcssPermEst <- function(i, while (!sufficient_data && trial < max_trials) { trial <- trial + 1 - d <- data. + d <- matlist. if (is.null(perm_var.)) { if (!large) { d[[dep]] <- RMPerm(d[[dep]], groups., CSS = TRUE) @@ -137,7 +137,7 @@ QAPcssPermEst <- function(i, names(x_list) <- data_vars. pred <- make_css_data(y = d[[dep]], x = x_list, nets = 1, - diag = diag., mode = mode.)$pred + diag = diag., directed = directed.)$pred } else { pred_list <- vector("list", length(d[[dep]])) for (gr in seq_along(d[[dep]])) { @@ -145,7 +145,7 @@ QAPcssPermEst <- function(i, names(xgr) <- data_vars. pred_list[[gr]] <- make_css_data(y = d[[dep]][[gr]], x = xgr, nets = gr, - diag = diag., mode = mode.)$pred + diag = diag., directed = directed.)$pred } pred <- do.call(rbind, pred_list) } @@ -200,7 +200,7 @@ QAPcssPermEst <- function(i, xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL if (is.null(comp.)) { - # A fit inside the permutation loop runs `reps` times, so a fitter's + # A fit inside the permutation loop runs `times` times, so a fitter's # convergence warning would print once per draw and drown the console. # The count of draws that failed outright is reported by # `aggregate_perm_results()`, which is the number the user needs. @@ -230,7 +230,7 @@ QAPcssPermEst <- function(i, predK <- pred[pred[[dep]] %in% comp.[[k]], ] predK[[dep]] <- ifelse(predK[[dep]] == comp.[[k]][1], 0, 1) - # A fit inside the permutation loop runs `reps` times, so a fitter's + # A fit inside the permutation loop runs `times` times, so a fitter's # convergence warning would print once per draw and drown the console. # The count of draws that failed outright is reported by # `aggregate_perm_results()`, which is the number the user needs. @@ -272,13 +272,13 @@ glm_tab <- function(x, comp) { cmat[, 2] <- format(x$lower[[comp]][2, ]) cmat[, 3] <- format(x$larger[[comp]][2, ]) cmat[, 4] <- format(x$abs[[comp]][2, ]) - if (x$nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (x$permute == "predictor") cmat[1, 2:4] <- "*" colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- names(x$base[[comp]]$coefficients) print.table(cmat) - if (x$nullhyp == "qapspp") - cat("\n* Significance test for the intercept is undefined with qapspp.\n") + if (x$permute == "predictor") + cat("\n* The intercept has no significance test when predictors are permuted.\n") if (!is.null(x$base[[comp]]$base_model)) { cat("\nAIC of base model:", format(stats::AIC(x$base[[comp]]$base_model))) @@ -294,13 +294,13 @@ glm_tab <- function(x, comp) { cmat[, 2] <- format(x$lower[2, ]) cmat[, 3] <- format(x$larger[2, ]) cmat[, 4] <- format(x$abs[2, ]) - if (x$nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (x$permute == "predictor") cmat[1, 2:4] <- "*" colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- names(x$base$coefficients) print.table(cmat) - if (x$nullhyp == "qapspp") - cat("\n* Significance test for the intercept is undefined with qapspp.\n") + if (x$permute == "predictor") + cat("\n* The intercept has no significance test when predictors are permuted.\n") if (!is.null(x$base$base_model)) { cat("\nAIC of base model:", format(stats::AIC(x$base$base_model))) @@ -314,11 +314,11 @@ glm_tab <- function(x, comp) { #' @keywords internal #' @noRd QAPcss <- function(formula, - data, - mode = "directed", + matlist, + directed = TRUE, diag = FALSE, - nullhyp = "qapy", - reps = 1000, + permute = "outcome", + times = 1000, seed = NULL, strategy = "sequential", ncores = NULL, @@ -340,20 +340,20 @@ QAPcss <- function(formula, parsed <- parse_qap_formula(formula, fixest_se_cluster) dep <- parsed$dependent main <- parsed$main - data_vars <- intersect(parsed$all_data_vars, names(data)) + data_vars <- intersect(parsed$all_data_vars, names(matlist)) nx <- length(main) - validate_qap_input(data, parsed, css = TRUE) - large <- is.list(data[[dep]]) + validate_qap_input(matlist, parsed, css = TRUE) + large <- is.list(matlist[[dep]]) if (!large) { - y <- data[[dep]] + y <- matlist[[dep]] if (length(dim(y)) != 3) manynet::snet_abort( "The dependent variable {.val {dep}} must be a 3-dimensional array of sender, receiver, and perceiver.") } else { - for (i in seq_along(data[[dep]])) { - if (length(dim(data[[dep]][[i]])) != 3) + for (i in seq_along(matlist[[dep]])) { + if (length(dim(matlist[[dep]][[i]])) != 3) manynet::snet_abort( "Network {i} of the dependent variable {.val {dep}} must be a 3-dimensional array.") } @@ -390,8 +390,8 @@ QAPcss <- function(formula, "Robust standard errors are not implemented for the multinomial family.") use_robust_errors <- FALSE } - if ((nullhyp == "qapspp") && (nx == 1)) nullhyp <- "qapy" - if (mode == "undirected" && (ris || rir)) { + if ((permute == "predictor") && (nx == 1)) permute <- "outcome" + if (!directed && (ris || rir)) { manynet::snet_warn( c("An undirected network has no senders or receivers.", i = "Setting the sender and receiver random intercepts to {.val FALSE}.")) @@ -408,7 +408,7 @@ QAPcss <- function(formula, if (rir) rand_part <- paste(rand_part, "+ (1|rv)") if (!large) { - n <- dim(data[[dep]])[1] + n <- dim(matlist[[dep]])[1] if (!is.null(groups)) { if (length(groups) != n) manynet::snet_abort( @@ -421,22 +421,22 @@ QAPcss <- function(formula, valid <- NULL; valid_list <- NULL if (!large) { - x_list <- lapply(data_vars, function(v) data[[v]]) + x_list <- lapply(data_vars, function(v) matlist[[v]]) names(x_list) <- data_vars - cssd <- make_css_data(y = data[[dep]], x = x_list, + cssd <- make_css_data(y = matlist[[dep]], x = x_list, nets = 1, - diag = diag, mode = mode) + diag = diag, directed = directed) pred <- cssd$pred valid <- cssd$valid } else { - pred_list <- vector("list", length(data[[dep]])) - valid_list <- vector("list", length(data[[dep]])) - for (gr in seq_along(data[[dep]])) { - xgr <- lapply(data_vars, function(v) data[[v]][[gr]]) + pred_list <- vector("list", length(matlist[[dep]])) + valid_list <- vector("list", length(matlist[[dep]])) + for (gr in seq_along(matlist[[dep]])) { + xgr <- lapply(data_vars, function(v) matlist[[v]][[gr]]) names(xgr) <- data_vars - cssd <- make_css_data(y = data[[dep]][[gr]], x = xgr, + cssd <- make_css_data(y = matlist[[dep]][[gr]], x = xgr, nets = gr, - diag = diag, mode = mode) + diag = diag, directed = directed) pred_list[[gr]] <- cssd$pred valid_list[[gr]] <- cssd$valid } @@ -480,20 +480,20 @@ QAPcss <- function(formula, if (use_gpu && family == "gaussian" && !has_random && !use_fixest && is.null(comparison) && !large) { - if (nullhyp == "qapy") { - gpu_res <- gpu_batch_ols_css(data = data, + if (permute == "outcome") { + gpu_res <- gpu_batch_ols_css(matlist = matlist, parsed = parsed, - mode = mode, + directed = directed, diag = diag, groups = groups, - reps = reps, + times = times, baseline_fit = fit$base, perm_var = NULL) fit$lower <- gpu_res$lower fit$larger <- gpu_res$larger fit$abs <- gpu_res$abs - } else if (nullhyp == "qapspp") { + } else if (permute == "predictor") { n_coefs <- length(fit$base$coefficients) fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) fit$larger <- fit$abs <- fit$lower @@ -501,7 +501,7 @@ QAPcss <- function(formula, colnames(fit$abs) <- names(fit$base$coefficients) for (xi in main) { - test_val <- data[[xi]] + test_val <- matlist[[xi]] if (!is.numeric(test_val)) { manynet::snet_warn( c("Cannot residualise the non-numeric predictor {.val {xi}}.", @@ -511,16 +511,16 @@ QAPcss <- function(formula, xR <- residualise_predictor(xi, pred, main, has_random = has_random, rand_formula = rand_part) - data_resid <- data - data_resid[[xi]] <- residuals_to_array(xR, data[[xi]], valid, pred, + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, large, valid_list) - gpu_res <- gpu_batch_ols_css(data = data_resid, + gpu_res <- gpu_batch_ols_css(matlist = matlist_resid, parsed = parsed, - mode = mode, + directed = directed, diag = diag, groups = groups, - reps = reps, + times = times, baseline_fit = fit$base, perm_var = xi) fit$lower[, xi] <- gpu_res$lower[, xi] @@ -536,12 +536,12 @@ QAPcss <- function(formula, options(future.globals.maxSize = attr(old_plan, "old_maxSize")) }, add = TRUE) - if (nullhyp == "qapy") { + if (permute == "outcome") { res <- run_permutations( - reps, QAPcssPermEst, - data. = data, + times, QAPcssPermEst, + matlist. = matlist, perm_var. = NULL, - mode. = mode, + directed. = directed, diag. = diag, mod. = mod, groups. = groups, @@ -560,7 +560,7 @@ QAPcss <- function(formula, ) if (is.null(comparison)) { - agg <- aggregate_perm_results(res, reps) + agg <- aggregate_perm_results(res, times) fit$lower <- agg$lower fit$larger <- agg$larger fit$abs <- agg$abs @@ -580,7 +580,7 @@ QAPcss <- function(formula, } } - } else if (nullhyp == "qapspp") { + } else if (permute == "predictor") { if (is.null(comparison)) { if (family != "multinom") { n_coefs <- length(fit$base$coefficients) @@ -590,9 +590,9 @@ QAPcss <- function(formula, colnames(fit$abs) <- names(fit$base$coefficients) } else { ncat <- if (large) { - length(stats::na.omit(unique(as.vector(unlist(data[[dep]]))))) + length(stats::na.omit(unique(as.vector(unlist(matlist[[dep]]))))) } else { - length(stats::na.omit(unique(as.vector(data[[dep]])))) + length(stats::na.omit(unique(as.vector(matlist[[dep]])))) } n_coefs <- length(fit$base$coefficients) fit$lower <- matrix(NA, nrow = 2 * (ncat - 1), ncol = n_coefs) @@ -615,7 +615,7 @@ QAPcss <- function(formula, } for (xi in main) { - test_val <- if (!large) data[[xi]] else data[[xi]][[1]] + test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] if (!is.numeric(test_val)) { manynet::snet_warn( c("Cannot residualise the non-numeric predictor {.val {xi}}.", @@ -627,15 +627,15 @@ QAPcss <- function(formula, has_random = has_random, rand_formula = rand_part) - data_resid <- data - data_resid[[xi]] <- residuals_to_array(xR, data[[xi]], valid, pred, + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, large, valid_list) res <- run_permutations( - reps, QAPcssPermEst, - data. = data_resid, + times, QAPcssPermEst, + matlist. = matlist_resid, perm_var. = xi, - mode. = mode, + directed. = directed, diag. = diag, mod. = mod, groups. = groups, @@ -654,7 +654,7 @@ QAPcss <- function(formula, ) if (is.null(comparison)) { - agg <- aggregate_perm_results(res, reps) + agg <- aggregate_perm_results(res, times) fit$lower[, xi] <- agg$lower fit$larger[, xi] <- agg$larger fit$abs[, xi] <- agg$abs @@ -685,12 +685,12 @@ QAPcss <- function(formula, } } - fit$nullhyp <- nullhyp + fit$permute <- permute fit$family <- family fit$groups <- unique(unlist(groups)) fit$diag <- diag - fit$mode <- mode - fit$reps <- reps + fit$directed <- directed + fit$times <- times fit$reference <- reference fit$comp <- comparison fit$random <- c(sender = ris, @@ -742,12 +742,12 @@ print.QAPCSS <- function(x, ...) { if (!is.null(x$groups)) cat("Permutations were performed within groups only.\n") - if (x$nullhyp == "qapy") - cat("The outcome array Y was permuted", format(x$reps), "times.\n") - if (x$nullhyp == "qapspp") { + if (x$permute == "outcome") + cat("The outcome array Y was permuted", format(x$times), "times.\n") + if (x$permute == "predictor") { cat("Significance was estimated using Dekker's\n") cat(" 'semi-partialling plus' procedure with", - format(x$reps), "permutations.\n") + format(x$times), "permutations.\n") } if (x$robust_se) @@ -759,7 +759,8 @@ print.QAPCSS <- function(x, ...) { } else { cat("Diagonal values (loops) were ignored.\n") } - cat("The outcome was treated as", format(paste0(x$mode, ".")), "\n") + cat("The outcome was treated as", + format(paste0(.directed_label(x$directed), ".")), "\n") if (x$family != "multinom") { if (is.null(x$comp)) { @@ -781,15 +782,15 @@ print.QAPCSS <- function(x, ...) { cmat[, 2] <- format(x$lower[row_idx, ]) cmat[, 3] <- format(x$larger[row_idx, ]) cmat[, 4] <- format(x$abs[row_idx, ]) - if (x$nullhyp == "qapspp") cmat[1, 2:4] <- "*" + if (x$permute == "predictor") cmat[1, 2:4] <- "*" colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") rownames(cmat) <- colnames(x$base$coefficients) print.table(cmat) cat("\n\n") } - if (x$nullhyp == "qapspp") - cat("* Significance test for the intercept is undefined with qapspp.\n") + if (x$permute == "predictor") + cat("* The intercept has no significance test when predictors are permuted.\n") cat("\nAIC of base model:", format(stats::AIC(x$base$base_model))) cat("\nBIC of base model:", format(stats::BIC(x$base$base_model))) diff --git a/R/qap_engine.R b/R/qap_engine.R index 8290bf9..ca1a0a5 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -8,13 +8,13 @@ #' @keywords internal #' @noRd QAPglm <- function(formula, - data, + matlist, family = "gaussian", - mode = "directed", + directed = TRUE, diag = FALSE, - nullhyp = "qapspp", + permute = "predictor", estimator = "standard", - reps = 1000, + times = 1000, seed = NULL, groups = NULL, strategy = "sequential", @@ -34,12 +34,11 @@ QAPglm <- function(formula, parsed <- parse_qap_formula(formula, fixest_se_cluster) dep <- parsed$dependent main <- parsed$main - data_vars <- intersect(parsed$all_data_vars, names(data)) + data_vars <- intersect(parsed$all_data_vars, names(matlist)) - validate_qap_input(data, parsed, css = FALSE) - large <- is.list(data[[dep]]) + validate_qap_input(matlist, parsed, css = FALSE) + large <- is.list(matlist[[dep]]) - mode_internal <- if (mode == "directed") "digraph" else "graph" rin <- random_intercept_nets ris <- random_intercept_sender @@ -59,25 +58,25 @@ QAPglm <- function(formula, mod <- stats::as.formula(mod_str) if (!large) { - pred <- make_qap_data(y = data[[dep]], - x = data[data_vars], + pred <- make_qap_data(y = matlist[[dep]], + x = matlist[data_vars], g = groups, diag = diag, - mode = mode_internal, + directed = directed, net = 1, perm = FALSE, xi = NULL) } else { - pred_list <- vector("list", length(data[[dep]])) - for (net in seq_along(data[[dep]])) { - x2 <- lapply(data_vars, function(v) data[[v]][[net]]) + pred_list <- vector("list", length(matlist[[dep]])) + for (net in seq_along(matlist[[dep]])) { + x2 <- lapply(data_vars, function(v) matlist[[v]][[net]]) names(x2) <- data_vars g2 <- if (!is.null(groups)) groups[[net]] else NULL - pred_list[[net]] <- make_qap_data(y = data[[dep]][[net]], + pred_list[[net]] <- make_qap_data(y = matlist[[dep]][[net]], x = x2, g = g2, diag = diag, - mode = mode_internal, + directed = directed, net = net, perm = FALSE, xi = NULL) @@ -128,7 +127,17 @@ QAPglm <- function(formula, } } - if ((nullhyp == "qapspp") && (length(main) == 1)) nullhyp <- "qapy" + # Double semi-partialling residualises a predictor against the others, so + # with one predictor there are none and the scheme reduces to permuting the + # outcome. Say so: the result would otherwise report a scheme nobody chose. + if ((permute == "predictor") && (length(main) == 1)) { + permute <- "outcome" + # `snet_info()` pastes its arguments, so give it separate strings rather + # than a named vector: a named vector loses its bullets and runs together. + manynet::snet_info( + "Permuting {.val outcome}, not {.val predictor}:", + "with one predictor there is nothing to residualise it against.") + } # The GPU path is a shortcut, not a requirement, so an unmet condition falls # back to the CPU permutation loop rather than aborting. `gpu_available()` @@ -144,20 +153,20 @@ QAPglm <- function(formula, if (use_gpu) { - if (nullhyp == "qapy") { - gpu_res <- gpu_batch_ols(data = data, + if (permute == "outcome") { + gpu_res <- gpu_batch_ols(matlist = matlist, parsed = parsed, - mode = mode_internal, + directed = directed, diag = diag, groups = groups, - reps = reps, + times = times, baseline_fit = fit$base, perm_var = NULL) fit$lower <- gpu_res$lower fit$larger <- gpu_res$larger fit$abs <- gpu_res$abs - } else if (nullhyp == "qapspp") { + } else if (permute == "predictor") { n_coefs <- length(fit$base$coefficients) fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, dimnames = list(c("perm_coefs", "perm_t"), @@ -168,15 +177,15 @@ QAPglm <- function(formula, xR <- residualise_predictor(xi, pred, main, has_random = has_random, rand_formula = rand_part) - data_resid <- data - data_resid[[xi]] <- residuals_to_matrix(xR, data[[xi]], pred, large) + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) - gpu_res <- gpu_batch_ols(data = data_resid, + gpu_res <- gpu_batch_ols(matlist = matlist_resid, parsed = parsed, - mode = mode_internal, + directed = directed, diag = diag, groups = groups, - reps = reps, + times = times, baseline_fit = fit$base, perm_var = xi) fit$lower[, xi] <- gpu_res$lower[, xi] @@ -192,12 +201,12 @@ QAPglm <- function(formula, options(future.globals.maxSize = attr(old_plan, "old_maxSize")) }, add = TRUE) - if (nullhyp == "qapy") { + if (permute == "outcome") { res <- run_permutations( - reps, QAPglmPermEst, - data. = data, + times, QAPglmPermEst, + matlist. = matlist, perm_var. = NULL, - mode. = mode_internal, + directed. = directed, diag. = diag, mod. = mod, groups. = groups, @@ -216,7 +225,7 @@ QAPglm <- function(formula, ) if (is.null(comparison)) { - agg <- aggregate_perm_results(res, reps) + agg <- aggregate_perm_results(res, times) fit$lower <- agg$lower fit$larger <- agg$larger fit$abs <- agg$abs @@ -237,7 +246,7 @@ QAPglm <- function(formula, } } - } else if (nullhyp == "qapspp") { + } else if (permute == "predictor") { if (is.null(comparison)) { n_coefs <- length(fit$base$coefficients) fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, @@ -263,14 +272,14 @@ QAPglm <- function(formula, has_random = has_random, rand_formula = rand_part) - data_resid <- data - data_resid[[xi]] <- residuals_to_matrix(xR, data[[xi]], pred, large) + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) res <- run_permutations( - reps, QAPglmPermEst, - data. = data_resid, + times, QAPglmPermEst, + matlist. = matlist_resid, perm_var. = xi, - mode. = mode_internal, + directed. = directed, diag. = diag, mod. = mod, groups. = groups, @@ -289,7 +298,7 @@ QAPglm <- function(formula, ) if (is.null(comparison)) { - agg <- aggregate_perm_results(res, reps) + agg <- aggregate_perm_results(res, times) fit$lower[, xi] <- agg$lower fit$larger[, xi] <- agg$larger fit$abs[, xi] <- agg$abs @@ -340,11 +349,11 @@ QAPglm <- function(formula, } } - fit$nullhyp <- nullhyp + fit$permute <- permute fit$diag <- diag fit$family <- family - fit$mode <- mode - fit$reps <- reps + fit$directed <- directed + fit$times <- times fit$groups <- unique(unlist(groups)) fit$robust_se <- use_robust_errors fit$estimator <- estimator @@ -365,9 +374,9 @@ QAPglm <- function(formula, #' @keywords internal #' @noRd QAPglmPermEst <- function(i, - data., + matlist., perm_var., - mode., + directed., diag., mod., groups., @@ -385,9 +394,9 @@ QAPglmPermEst <- function(i, reference.) { dep <- parsed.$dependent - large <- is.list(data.[[dep]]) + large <- is.list(matlist.[[dep]]) - d <- data. + d <- matlist. if (is.null(perm_var.)) { if (!large) { d[[dep]] <- RMPerm(d[[dep]], groups.) @@ -407,7 +416,7 @@ QAPglmPermEst <- function(i, x = d[data_vars.], g = groups., diag = diag., - mode = mode., + directed = directed., net = 1, perm = FALSE, xi = NULL) @@ -421,7 +430,7 @@ QAPglmPermEst <- function(i, x = x2, g = g2, diag = diag., - mode = mode., + directed = directed., net = net, perm = FALSE, xi = NULL) @@ -434,7 +443,7 @@ QAPglmPermEst <- function(i, xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL if (is.null(comp.)) { - # A fit inside the permutation loop runs `reps` times, so a fitter's + # A fit inside the permutation loop runs `times` times, so a fitter's # convergence warning would print once per draw and drown the console. # The count of draws that failed outright is reported by # `aggregate_perm_results()`, which is the number the user needs. @@ -464,7 +473,7 @@ QAPglmPermEst <- function(i, predK <- pred[pred[[dep]] %in% comp.[[k]], ] predK[[dep]] <- ifelse(predK[[dep]] == comp.[[k]][1], 0, 1) - # A fit inside the permutation loop runs `reps` times, so a fitter's + # A fit inside the permutation loop runs `times` times, so a fitter's # convergence warning would print once per draw and drown the console. # The count of draws that failed outright is reported by # `aggregate_perm_results()`, which is the number the user needs. diff --git a/R/qap_gpu.R b/R/qap_gpu.R index 92cdf97..69a9658 100644 --- a/R/qap_gpu.R +++ b/R/qap_gpu.R @@ -5,7 +5,7 @@ #' @keywords internal #' @noRd -gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, +gpu_batch_ols <- function(matlist, parsed, directed, diag, groups, times, baseline_fit, perm_var = NULL, batch_size = 500, device = "cuda") { @@ -19,11 +19,11 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, dep <- parsed$dependent main <- parsed$main - pred0 <- make_qap_data(y = data[[dep]], - x = data[main], + pred0 <- make_qap_data(y = matlist[[dep]], + x = matlist[main], g = groups, diag = diag, - mode = mode, + directed = directed, net = 1, perm = FALSE, xi = NULL) @@ -51,17 +51,17 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, XtXinv_diag <- torch::torch_diag(XtXinv) reps_done <- 0 - while (reps_done < reps) { - current_batch <- min(batch_size, reps - reps_done) + while (reps_done < times) { + current_batch <- min(batch_size, times - reps_done) Y_batch <- matrix(NA_real_, nrow = n_obs, ncol = current_batch) for (j in seq_len(current_batch)) { - y_perm <- RMPerm(data[[dep]], groups) + y_perm <- RMPerm(matlist[[dep]], groups) perm_pred <- make_qap_data(y = y_perm, - x = data[main], + x = matlist[main], g = groups, diag = diag, - mode = mode, + directed = directed, net = 1, perm = FALSE, xi = NULL) @@ -97,21 +97,21 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, device = device) reps_done <- 0 - while (reps_done < reps) { - current_batch <- min(batch_size, reps - reps_done) + while (reps_done < times) { + current_batch <- min(batch_size, times - reps_done) B_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) T_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) for (j in seq_len(current_batch)) { - d_perm <- data + d_perm <- matlist d_perm[[perm_var]] <- RMPerm(d_perm[[perm_var]], groups) perm_pred <- make_qap_data(y = d_perm[[dep]], x = d_perm[main], g = groups, diag = diag, - mode = mode, + directed = directed, net = 1, perm = FALSE, xi = NULL) @@ -148,11 +148,11 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, } list( - lower = matrix(lower_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + lower = matrix(lower_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))), - larger = matrix(larger_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + larger = matrix(larger_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))), - abs = matrix(abs_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + abs = matrix(abs_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))) ) } @@ -160,7 +160,7 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, #' @keywords internal #' @noRd -gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, +gpu_batch_ols_css <- function(matlist, parsed, directed, diag, groups, times, baseline_fit, perm_var = NULL, batch_size = 500, device = "cuda") { @@ -175,10 +175,10 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, main <- parsed$main data_vars <- parsed$all_data_vars - x_list <- lapply(data_vars, function(v) data[[v]]) + x_list <- lapply(data_vars, function(v) matlist[[v]]) names(x_list) <- data_vars - cssd <- make_css_data(y = data[[dep]], x = x_list, - nets = 1, diag = diag, mode = mode) + cssd <- make_css_data(y = matlist[[dep]], x = x_list, + nets = 1, diag = diag, directed = directed) pred0 <- cssd$pred y_vec <- pred0$yv @@ -199,7 +199,7 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, xl <- lapply(data_vars, function(v) d[[v]]) names(xl) <- data_vars make_css_data(y = d[[dep]], x = xl, - nets = 1, diag = diag, mode = mode)$pred + nets = 1, diag = diag, directed = directed)$pred } if (is.null(perm_var)) { @@ -211,12 +211,12 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, XtXinv_diag <- torch::torch_diag(XtXinv) reps_done <- 0 - while (reps_done < reps) { - current_batch <- min(batch_size, reps - reps_done) + while (reps_done < times) { + current_batch <- min(batch_size, times - reps_done) Y_batch <- matrix(NA_real_, nrow = n_obs, ncol = current_batch) for (j in seq_len(current_batch)) { - d_perm <- data + d_perm <- matlist d_perm[[dep]] <- RMPerm(d_perm[[dep]], groups, CSS = TRUE) perm_pred <- build_css_pred(d_perm) Y_batch[, j] <- perm_pred$yv @@ -251,14 +251,14 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, device = device) reps_done <- 0 - while (reps_done < reps) { - current_batch <- min(batch_size, reps - reps_done) + while (reps_done < times) { + current_batch <- min(batch_size, times - reps_done) B_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) T_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) for (j in seq_len(current_batch)) { - d_perm <- data + d_perm <- matlist d_perm[[perm_var]] <- RMPerm(d_perm[[perm_var]], groups, CSS = TRUE) perm_pred <- build_css_pred(d_perm) @@ -295,11 +295,11 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, } list( - lower = matrix(lower_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + lower = matrix(lower_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))), - larger = matrix(larger_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + larger = matrix(larger_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))), - abs = matrix(abs_sum / reps, nrow = dim_out[1], ncol = dim_out[2], + abs = matrix(abs_sum / times, nrow = dim_out[1], ncol = dim_out[2], dimnames = list(NULL, names(base_coefs))) ) } diff --git a/R/qap_misc.R b/R/qap_misc.R index d362a9b..7a60122 100644 --- a/R/qap_misc.R +++ b/R/qap_misc.R @@ -15,42 +15,42 @@ combine_qap_estimates <- function(res, res2 = NULL) { return_res <- res[[1]] if (is.null(return_res$comp)) { for (i in 1:(n_res - 1)) { - return_res$lower <- (return_res$lower * return_res$reps + - res[[i + 1]]$lower * res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + return_res$lower <- (return_res$lower * return_res$times + + res[[i + 1]]$lower * res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) - return_res$larger <- (return_res$larger * return_res$reps + - res[[i + 1]]$larger * res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + return_res$larger <- (return_res$larger * return_res$times + + res[[i + 1]]$larger * res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) - return_res$abs <- (return_res$abs * return_res$reps + - res[[i + 1]]$abs * res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + return_res$abs <- (return_res$abs * return_res$times + + res[[i + 1]]$abs * res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) - return_res$reps <- return_res$reps + res[[i + 1]]$reps + return_res$times <- return_res$times + res[[i + 1]]$times } } else { for (i in 1:(n_res - 1)) { for (com in names(return_res$comp)) { return_res[[com]]$lower <- (return_res[[com]]$lower * - return_res$reps + + return_res$times + res[[i + 1]][[com]]$lower * - res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) return_res[[com]]$larger <- (return_res[[com]]$larger * - return_res$reps + + return_res$times + res[[i + 1]][[com]]$larger * - res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) return_res[[com]]$abs <- (return_res[[com]]$abs * - return_res$reps + + return_res$times + res[[i + 1]][[com]]$abs * - res[[i + 1]]$reps) / - (return_res$reps + res[[i + 1]]$reps) + res[[i + 1]]$times) / + (return_res$times + res[[i + 1]]$times) } - return_res$reps <- return_res$reps + res[[i + 1]]$reps + return_res$times <- return_res$times + res[[i + 1]]$times } } @@ -64,11 +64,10 @@ df_to_mat <- function(df, sender, receiver, perceiver = NULL, - mode = c("directed", "undirected"), + directed = TRUE, loops = FALSE, multi_mode = FALSE, split_by = NULL) { - mode <- match.arg(mode) var_names <- setdiff(colnames(df), c(sender, receiver, perceiver, split_by)) if (!is.null(split_by)) { @@ -77,7 +76,7 @@ df_to_mat <- function(df, sender = sender, receiver = receiver, perceiver = perceiver, - mode = mode, + directed = directed, loops = loops, multi_mode = multi_mode) return(purrr::transpose(result)) @@ -97,7 +96,7 @@ df_to_mat <- function(df, n_r <- length(nodes_r) n_p <- if (!is.null(perceiver)) length(nodes_p) else NULL - expected <- if (mode == "undirected") { + expected <- if (!directed) { if (loops) n_s * (n_s + 1) / 2 else n_s * (n_s - 1) / 2 } else { if (loops) n_s * n_r else n_s * n_r - min(n_s, n_r) @@ -113,7 +112,7 @@ df_to_mat <- function(df, mat <- matrix(NA_real_, nrow = n_s, ncol = n_r, dimnames = list(nodes_s, nodes_r)) mat[cbind(df[[sender]], df[[receiver]])] <- df[[var]] - if (mode == "undirected") + if (!directed) mat[cbind(df[[receiver]], df[[sender]])] <- df[[var]] if (!loops) diag(mat) <- NA mat @@ -121,7 +120,7 @@ df_to_mat <- function(df, arr <- array(NA_real_, dim = c(n_s, n_r, n_p), dimnames = list(nodes_s, nodes_r, nodes_p)) arr[cbind(df[[sender]], df[[receiver]], df[[perceiver]])] <- df[[var]] - if (mode == "undirected") + if (!directed) arr[cbind(df[[receiver]], df[[sender]], df[[perceiver]])] <- df[[var]] if (!loops && !multi_mode) arr[cbind(nodes_s, nodes_s, rep(nodes_p, each = n_s))] <- NA @@ -131,3 +130,13 @@ df_to_mat <- function(df, stats::setNames(lapply(var_names, make_structure), var_names) } + + +# The fit records directedness as a logical, because that is what +# `manynet::is_directed()` returns and what the engine branches on. Users read +# the word, so the print methods render it here rather than each spelling it. +#' @keywords internal +#' @noRd +.directed_label <- function(directed) { + if (isTRUE(directed)) "directed" else "undirected" +} diff --git a/R/qap_utils.R b/R/qap_utils.R index 2afd0c5..837d06a 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -75,20 +75,20 @@ build_internal_formula <- function(formula, #' @keywords internal #' @noRd -validate_qap_input <- function(data, parsed, css = FALSE) { +validate_qap_input <- function(matlist, parsed, css = FALSE) { dep <- parsed$dependent - if (!(dep %in% names(data))) { + if (!(dep %in% names(matlist))) { manynet::snet_abort("Dependent variable {.val {dep}} not found in the data.") } structural_vars <- c("sv", "rv", "nv", "pv") for (v in parsed$all_data_vars) { if (v %in% structural_vars) next - if (!(v %in% names(data))) { + if (!(v %in% names(matlist))) { manynet::snet_abort("Predictor {.val {v}} not found in the data.") } } - y <- data[[dep]] + y <- matlist[[dep]] large <- is.list(y) if (!css) { @@ -139,9 +139,9 @@ setup_future_plan <- function(strategy = "sequential", ncores = NULL) { #' @keywords internal #' @noRd -run_permutations <- function(reps, FUN, ...) { +run_permutations <- function(times, FUN, ...) { future.apply::future_lapply( - seq_len(reps), + seq_len(times), FUN, ..., future.seed = TRUE @@ -211,7 +211,7 @@ RMPerm <- function(m, groups = NULL, CSS = FALSE) { #' @keywords internal #' @noRd -make_qap_data <- function(y, x, g = NULL, diag = FALSE, mode = "digraph", +make_qap_data <- function(y, x, g = NULL, diag = FALSE, directed = TRUE, net = 1, perm = FALSE, xi = NULL) { nx <- length(x) @@ -235,7 +235,7 @@ make_qap_data <- function(y, x, g = NULL, diag = FALSE, mode = "digraph", # the diagonal. Keeping both halves doubles the sample and shrinks every # standard error, so take the lower triangle only. A two-mode incidence # matrix has no such symmetry, and keeps every cell. - if (identical(mode, "graph") && square) valid[upper.tri(valid)] <- FALSE + if (!directed && square) valid[upper.tri(valid)] <- FALSE for (var in seq_len(nx)) { valid[is.na(x[[var]])] <- FALSE @@ -593,16 +593,16 @@ compare_perm_to_baseline <- function(perm_coefs, perm_t, base_fit, #' @keywords internal #' @noRd -aggregate_perm_results <- function(results, reps) { +aggregate_perm_results <- function(results, times) { results <- Filter(Negate(is.null), results) n_valid <- length(results) if (n_valid == 0) manynet::snet_abort( - c("All {reps} permutations failed to converge.", + c("All {times} permutations failed to converge.", i = "Try a simpler model, another {.arg family}, or fewer predictors.")) - if (n_valid < reps) { + if (n_valid < times) { manynet::snet_warn( - "{reps - n_valid} of {reps} permutation{?s} failed and {?was/were} excluded.") + "{times - n_valid} of {times} permutation{?s} failed and {?was/were} excluded.") } resL <- unlist(results, recursive = FALSE) list( @@ -613,7 +613,7 @@ aggregate_perm_results <- function(results, reps) { } -# ---- residualisation for qapspp --------------------------------------------- +# ---- residualisation for permute = "predictor" --------------------------------------------- #' @keywords internal #' @noRd diff --git a/man/regression.Rd b/man/regression.Rd index 0b5c0bd..7628cec 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -37,16 +37,20 @@ dropped with a warning.} \item{control}{Named list of additional controls; unspecified entries fall back to the defaults below. \itemize{ -\item \code{method}: \code{"qap"} (double semi-partialling plus, default) or \code{"qapy"} -(permute y only). +\item \code{permute}: what the null distribution permutes. \code{"predictor"} (the +default) residualises each main predictor against the others and +permutes that residual, following Dekker et al. (2007). \code{"outcome"} +permutes the dependent matrix and leaves the predictors alone. With one +predictor there is nothing to residualise against, so \code{"predictor"} +reduces to \code{"outcome"} and says so. \item \code{strategy}: future plan, e.g. \code{"sequential"} (default), \code{"multisession"}. \item \code{family}: \code{"auto"} (default; gaussian for weighted networks, binomial for binary), \code{"gaussian"}, \code{"binomial"}, \code{"poisson"}, \code{"negbin"}, \code{"zip"}, or \code{"multinom"}. \item \code{estimator}: \code{"standard"} (default) or \code{"gmm"} (binomial/poisson/ negbin/zip). -\item \code{mode}: \code{"directed"} / \code{"undirected"} (default auto-detected from -\code{.data}). +\item \code{directed}: logical, whether a tie from i to j differs from one from j +to i. Read from \code{.data} unless given, and reported when read. \item \code{diag}: logical, include loops (default auto-detected). \item \code{seed}, \code{groups}, \code{ncores}: passed through to the engine. \item \code{use_robust_errors}: HC3 standard errors. @@ -85,8 +89,8 @@ and handed to a QAP engine (ported from MrQAP) that supports: \itemize{ \item gaussian, binomial, poisson, negbin, zero-inflated Poisson, and multinomial families; -\item \code{"qap"} (Dekker's double semi-partialling plus) and \code{"qapy"} -(permute-y-only) null hypotheses; +\item two permutation schemes: \code{"predictor"} (Dekker's double semi-partialling) +and \code{"outcome"}; \item random intercepts (lme4 / glmmTMB) and fixed effects (fixest); \item robust (HC3) standard errors; \item optional torch-based batch OLS on the GPU; diff --git a/tests/testthat/helper-infernet.R b/tests/testthat/helper-infernet.R index 54f16f1..ad4b002 100644 --- a/tests/testthat/helper-infernet.R +++ b/tests/testthat/helper-infernet.R @@ -58,6 +58,20 @@ qap_net_undirected <- function(n = 24, seed = 105) { manynet::mutate(manynet::as_tidygraph(m, twomode = FALSE), Age = age) } +# A two-mode network with more columns than rows. This shape is what turned the +# square-matrix assumption into an error rather than a wrong number: the +# predictor held more cells than the validity mask, and a logical index longer +# than its target extends that target with NA. See stocnet/infernet#4. +qap_net_twomode_wide <- function(nr = 12, nc = 40, seed = 107) { + set.seed(seed) + m <- matrix(stats::rbinom(nr * nc, size = 4, prob = 0.3), nr, nc) + g <- manynet::as_igraph(m) + manynet::mutate_nodes(g, + GONGO = c(rep(c("GON", "GO"), length.out = nr), rep(NA_character_, nc)), + province = c(rep(LETTERS[1:4], length.out = nr), rep(NA_character_, nc)), + Att = c(stats::runif(nr), rep(NA_real_, nc))) +} + qap_net_twomode <- function(seed = 106) { set.seed(seed) sw <- manynet::ison_southern_women @@ -70,16 +84,14 @@ qap_net_twomode <- function(seed = 106) { # coefficient can be compared against the equivalent standard fit on identical # data. Anything this returns comes from the engine's own internals, so a # comparison against it tests the estimator dispatch, not the vectorisation. -qap_reference_data <- function(formula, .data, mode = NULL, diag = FALSE) { +qap_reference_data <- function(formula, .data, directed = NULL, diag = FALSE) { ml <- convertToMatrixList(formula, .data, advise = FALSE) parsed <- parse_qap_formula(ml$formula) g <- manynet::as_tidygraph(.data) - if (is.null(mode)) { - mode <- if (manynet::is_directed(g)) "digraph" else "graph" - } + if (is.null(directed)) directed <- manynet::is_directed(g) pred <- make_qap_data(y = ml$mydata[[parsed$dependent]], x = ml$mydata[parsed$main], - diag = diag, mode = mode) + diag = diag, directed = directed) names(pred)[names(pred) == "yv"] <- parsed$dependent list(pred = pred, formula = ml$formula, parsed = parsed) } diff --git a/tests/testthat/test-model_tests.R b/tests/testthat/test-model_tests.R index 360dc9f..05ac314 100644 --- a/tests/testthat/test-model_tests.R +++ b/tests/testthat/test-model_tests.R @@ -16,12 +16,12 @@ cugtest2 <- test_random(marvel_friends, test_that("test_random works", { expect_equal(as.numeric(cugtest$testval), -0.85714, tolerance = 0.001) expect_length(cugtest$testdist, 200) # NB: Stochastic - expect_false(cugtest$mode) + expect_false(cugtest$directed) expect_false(cugtest$diag) expect_equal(cugtest$cmode, "edges") expect_type(cugtest$plteobs, "double") expect_type(cugtest$pgteobs, "double") - expect_equal(cugtest$reps, 200) + expect_equal(cugtest$times, 200) expect_s3_class(cugtest, "network_test") expect_equal(as.numeric(cugtest2$testval), 0.2375, tolerance = 0.001) expect_length(cugtest2$testdist, 200) # NB: Stochastic @@ -40,7 +40,7 @@ test_that("test_permutation works", { expect_type(qaptest$plteobs, "double") # NB: Stochastic expect_type(qaptest$pgteobs, "double") # NB: Stochastic expect_length(qaptest$testdist, 200) # NB: Stochastic - expect_equal(qaptest$reps, 200) + expect_equal(qaptest$times, 200) expect_s3_class(qaptest, "network_test") }) diff --git a/tests/testthat/test-net_regression.R b/tests/testthat/test-net_regression.R index bf61f99..9aef90d 100644 --- a/tests/testthat/test-net_regression.R +++ b/tests/testthat/test-net_regression.R @@ -139,12 +139,12 @@ test_that("print.net_regression runs without error for QAPGLM", { # ---- method control -------------------------------------------------------- -test_that("method = 'qapy' runs and flags the nullhyp on the fit", { +test_that("permute = 'outcome' runs and is recorded on the fit", { g <- make_weighted_net() fit <- net_regression(weight ~ ego(Age) + alter(Age), g, times = 10, - control = list(method = "qapy")) - expect_equal(fit$nullhyp, "qapy") + control = list(permute = "outcome")) + expect_equal(fit$permute, "outcome") }) diff --git a/tests/testthat/test-qap_control.R b/tests/testthat/test-qap_control.R index fcb98bd..e76593e 100644 --- a/tests/testthat/test-qap_control.R +++ b/tests/testthat/test-qap_control.R @@ -23,7 +23,7 @@ test_that("an unnamed control entry is rejected", { test_that("an empty control list gives the defaults", { expect_equal(.resolve_control(list()), .resolve_control()) - expect_equal(.resolve_control()$method, "qap") + expect_equal(.resolve_control()$permute, "predictor") expect_equal(.resolve_control()$strategy, "sequential") expect_equal(.resolve_control()$family, "auto") }) @@ -35,41 +35,40 @@ test_that("a named control overrides only that default", { expect_equal(ctrl$estimator, "standard") }) -test_that("method takes only the two spellings it documents", { - expect_equal(.resolve_control(list(method = "qapy"))$method, "qapy") - expect_error(.resolve_control(list(method = "spp"))) +test_that("permute takes only the two spellings it documents", { + expect_equal(.resolve_control(list(permute = "outcome"))$permute, "outcome") + expect_error(.resolve_control(list(permute = "qapy"))) }) -test_that("method = 'qapy' is recorded on the fit and gives a full matrix", { +test_that("permute = 'outcome' is recorded on the fit and gives a full matrix", { g <- qap_net_gaussian(n = 20) fit <- net_regression(FORM, g, times = 10, - control = list(seed = 1, method = "qapy")) - expect_equal(fit$nullhyp, "qapy") + control = list(seed = 1, permute = "outcome")) + expect_equal(fit$permute, "outcome") # Permuting y alone tests every coefficient, the intercept included, whereas # double semi-partialling residualises one predictor at a time. expect_false(anyNA(fit$lower)) spp <- net_regression(FORM, g, times = 10, control = list(seed = 1)) - expect_equal(spp$nullhyp, "qapspp") + expect_equal(spp$permute, "predictor") expect_true(all(is.na(spp$lower[, "(Intercept)"]))) }) -test_that("a single predictor falls back from qapspp to qapy", { +test_that("a single predictor falls back to permuting the outcome", { g <- qap_net_gaussian(n = 20) fit <- net_regression(weight ~ ego(Age), g, times = 10, - control = list(seed = 1, method = "qap")) + control = list(seed = 1, permute = "predictor")) # Double semi-partialling residualises a predictor against the others, and # with one predictor there are none. - expect_equal(fit$nullhyp, "qapy") + expect_equal(fit$permute, "outcome") }) -test_that("mode and diag are read from the network unless set", { +test_that("directed and diag are read from the network unless set", { g <- qap_net_gaussian(n = 15) - expect_equal(net_regression(FORM, g, times = 5, - control = list(seed = 1))$mode, "directed") - expect_equal(net_regression(FORM, g, times = 5, + expect_true(net_regression(FORM, g, times = 5, + control = list(seed = 1))$directed) + expect_false(net_regression(FORM, g, times = 5, control = list(seed = 1, - mode = "undirected"))$mode, - "undirected") + directed = FALSE))$directed) loops <- net_regression(FORM, g, times = 5, control = list(seed = 1, diag = TRUE)) expect_true(loops$diag) diff --git a/tests/testthat/test-qap_reporting.R b/tests/testthat/test-qap_reporting.R new file mode 100644 index 0000000..5af0278 --- /dev/null +++ b/tests/testthat/test-qap_reporting.R @@ -0,0 +1,74 @@ +# Anything the model resolves for itself is reported, so that a user can +# describe the model they actually fitted. +# +# These run with `snet_verbosity = "verbose"`, which is not the default. That +# matters: an informational message is silent in every other test, so a broken +# one stays invisible. A `{cli}` brace expression beginning with a dot is read +# as a style rather than as code, and two of these messages shipped that way +# before this file existed. + +verbosely <- function(expr) { + old <- options(snet_verbosity = "verbose") + on.exit(options(old), add = TRUE) + force(expr) +} + +FORM <- weight ~ ego(Age) + alter(Age) + sim(Age) + +test_that("every reporting message renders", { + g <- qap_net_gaussian(n = 15) + # A message that cli cannot parse aborts, so reaching the end is the test. + expect_no_error(verbosely( + net_regression(FORM, g, times = 5, control = list(seed = 1)))) + expect_no_error(verbosely( + net_regression(weight ~ ego(Age), g, times = 5, control = list(seed = 1)))) + expect_no_error(verbosely( + net_regression(FORM, g, times = 5, + control = list(seed = 1, use_gpu = TRUE)))) +}) + +test_that("a family resolved from the outcome is reported", { + expect_message( + verbosely(net_regression(. ~ ego(Age) + alter(Age), qap_net_binary(n = 15), + times = 5, control = list(seed = 1))), + "binomial") +}) + +test_that("directedness read from the network is reported", { + expect_message( + verbosely(net_regression(weight ~ ego(Age) + alter(Age), + qap_net_undirected(n = 14), + times = 5, control = list(seed = 1))), + "undirected") +}) + +test_that("a stated family and directedness are not reported", { + g <- qap_net_gaussian(n = 15) + expect_no_message( + verbosely(net_regression(FORM, g, times = 5, + control = list(seed = 1, family = "gaussian", + directed = TRUE)))) +}) + +test_that("the fallback to permuting the outcome is reported", { + expect_message( + verbosely(net_regression(weight ~ ego(Age), qap_net_gaussian(n = 15), + times = 5, + control = list(seed = 1, permute = "predictor"))), + "residualise") +}) + +test_that("the GPU falling back to the CPU is reported", { + skip_if(gpu_available(), "a CUDA device is present, so there is no fallback") + expect_message( + verbosely(net_regression(FORM, qap_net_gaussian(n = 15), times = 5, + control = list(seed = 1, use_gpu = TRUE))), + "CPU") +}) + +test_that("the model advice on homophily terms is reported", { + expect_message( + verbosely(net_regression(weight ~ same(Grp), qap_net_gaussian(n = 15), + times = 5, control = list(seed = 1))), + "ego\\(Grp\\)") +}) diff --git a/tests/testthat/test-qap_shapes.R b/tests/testthat/test-qap_shapes.R index f821c93..f208156 100644 --- a/tests/testthat/test-qap_shapes.R +++ b/tests/testthat/test-qap_shapes.R @@ -6,7 +6,7 @@ test_that("a directed network contributes every ordered dyad", { g <- qap_net_gaussian(n = 20) fit <- net_regression(weight ~ ego(Age), g, times = 10, control = list(seed = 1)) - expect_equal(fit$mode, "directed") + expect_true(fit$directed) expect_equal(nrow(fit$pred), 20 * 19) }) @@ -16,7 +16,7 @@ test_that("an undirected network contributes each dyad once", { fit <- net_regression(weight ~ ego(Age), g, times = 10, control = list(seed = 1)) # Both halves of a symmetric matrix hold the same dyad. Keeping both doubles # the sample and shrinks every standard error by about a factor of root two. - expect_equal(fit$mode, "undirected") + expect_false(fit$directed) expect_equal(nrow(fit$pred), 24 * 23 / 2) }) @@ -31,6 +31,28 @@ test_that("a two-mode network contributes every cell of the incidence matrix", { expect_named(fit$coefficients, c("(Intercept)", "ego Att", "alter Att")) }) +test_that("a wide two-mode network fits, and counts its dyads", { + # stocnet/infernet#4: the validity mask was built as nrow-by-nrow, so a + # predictor with more columns than rows extended it with NA, and the dyad + # count came back as NA rather than a number. + g <- qap_net_twomode_wide(nr = 12, nc = 40) + expect_true(manynet::is_twomode(g)) + fit <- net_regression(weight ~ same(GONGO) + same(province), g, times = 10, + control = list(seed = 1)) + expect_equal(nrow(fit$pred), 12 * 40) + expect_false(anyNA(fit$coefficients)) + expect_named(fit$coefficients, + c("(Intercept)", "same GONGO", "same province")) +}) + +test_that("make_qap_data() counts a wide predictor's cells, not the mask's", { + y <- matrix(stats::rnorm(6 * 15), 6, 15) + x <- list(a = matrix(stats::rnorm(6 * 15), 6, 15)) + pred <- make_qap_data(y = y, x = x, diag = FALSE, directed = TRUE) + expect_equal(nrow(pred), 6 * 15) + expect_false(anyNA(pred$a)) +}) + test_that("RMPerm() permutes a rectangular matrix without erroring", { m <- matrix(seq_len(18 * 14), 18, 14) p <- RMPerm(m) From 14b67391fa9bdf6b7fd441bbd882da787500cee8 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:21:48 +0200 Subject: [PATCH 02/12] Branched off the torch GPU path Moved to feature/torch-gpu. Reinstate with `git revert` of this commit. `gpu_batch_ols()` ran only where the family is gaussian and there were no random effects, no fixed effects, and no multinomial comparison. It had no test, no hosted runner has a CUDA device to exercise it, and having {torch} in Suggests broke the CI build: torch installs as an R package before its Lantern backend, so `cuda_is_available()` throws rather than returning FALSE. Removes R/qap_gpu.R, the `use_gpu` control, and the gate in both engines. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 1 - R/model_regression.R | 8 +- R/qap_css.R | 264 +++++++++------------- R/qap_engine.R | 216 +++++++----------- R/qap_gpu.R | 318 --------------------------- man/regression.Rd | 2 - tests/testthat/test-qap_estimators.R | 13 -- tests/testthat/test-qap_reporting.R | 10 - 8 files changed, 186 insertions(+), 646 deletions(-) delete mode 100644 R/qap_gpu.R diff --git a/DESCRIPTION b/DESCRIPTION index 8740e49..1eeb735 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -45,7 +45,6 @@ Suggests: MASS, pscl, testthat (>= 3.0.0), - torch, glmmTMB Config/Needs/build: roxygen2, diff --git a/R/model_regression.R b/R/model_regression.R index 023457f..36293f6 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -16,7 +16,6 @@ #' and `"outcome"`; #' - random intercepts (lme4 / glmmTMB) and fixed effects (fixest); #' - robust (HC3) standard errors; -#' - optional torch-based batch OLS on the GPU; #' - lists of networks, in which graphs that are missing any predictor are #' dropped with a warning and the remaining networks are pooled. #' @@ -61,7 +60,6 @@ #' - `reference`, `comparison`: multinomial / pairwise-comparison options. #' - `random_intercept_nets` / `_sender` / `_receiver`: lme4-style REs. #' - `less_mem`: drop the baseline model object from the return. -#' - `use_gpu`: torch-based batch OLS (gaussian only). #' @return An object of class `net_regression` inheriting from either #' `QAPRegression` (gaussian) or `QAPGLM` (other families). When the #' outcome is binary -- either `family = "binomial"` or `"gaussian"` with @@ -155,8 +153,7 @@ net_regression <- function(formula, random_intercept_sender = ctrl$random_intercept_sender, random_intercept_receiver = ctrl$random_intercept_receiver, use_robust_errors = ctrl$use_robust_errors, - less_mem = ctrl$less_mem, - use_gpu = ctrl$use_gpu + less_mem = ctrl$less_mem ) if (user_requested_gaussian_binary && is.null(ctrl$comparison)) { @@ -236,8 +233,7 @@ net_regression <- function(formula, random_intercept_nets = FALSE, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, - less_mem = FALSE, - use_gpu = FALSE + less_mem = FALSE ) } diff --git a/R/qap_css.R b/R/qap_css.R index c003bee..0abc8b0 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -333,7 +333,7 @@ QAPcss <- function(formula, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, random_intercept_perceiver = FALSE, - use_gpu = FALSE) { + less_mem = FALSE) { if (!is.null(seed)) set.seed(seed) @@ -477,70 +477,111 @@ QAPcss <- function(formula, } } - if (use_gpu && family == "gaussian" && !has_random && !use_fixest && - is.null(comparison) && !large) { - - if (permute == "outcome") { - gpu_res <- gpu_batch_ols_css(matlist = matlist, - parsed = parsed, - directed = directed, - diag = diag, - groups = groups, - times = times, - baseline_fit = fit$base, - perm_var = NULL) - fit$lower <- gpu_res$lower - fit$larger <- gpu_res$larger - fit$abs <- gpu_res$abs - - } else if (permute == "predictor") { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - - for (xi in main) { - test_val <- matlist[[xi]] - if (!is.numeric(test_val)) { - manynet::snet_warn( - c("Cannot residualise the non-numeric predictor {.val {xi}}.", - i = "Skipping double semi-partialling for this predictor.")) - next + old_plan <- setup_future_plan(strategy, ncores) + on.exit({ + future::plan(old_plan) + options(future.globals.maxSize = attr(old_plan, "old_maxSize")) + }, add = TRUE) + + if (permute == "outcome") { + res <- run_permutations( + times, QAPcssPermEst, + matlist. = matlist, + perm_var. = NULL, + directed. = directed, + diag. = diag, + mod. = mod, + groups. = groups, + fit. = if (is.null(comparison)) fit$base else fit$base, + family. = family, + estimator. = estimator, + use_fixest. = use_fixest, + fixest_se_cluster. = fixest_se_cluster, + use_robust_errors. = use_robust_errors, + has_random. = has_random, + main_vars. = main, + data_vars. = data_vars, + parsed. = parsed, + comp. = comparison, + reference. = reference + ) + + if (is.null(comparison)) { + agg <- aggregate_perm_results(res, times) + fit$lower <- agg$lower + fit$larger <- agg$larger + fit$abs <- agg$abs + } else { + res_valid <- Filter(Negate(is.null), res) + n_valid <- length(res_valid) + fit$lower <- fit$larger <- fit$abs <- + vector("list", length(comparison)) + names(fit$lower) <- names(fit$larger) <- + names(fit$abs) <- names(comparison) + resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) + for (k in seq_along(comparison)) { + cn <- names(comparison)[k] + fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid + fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid + fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid + } + } + + } else if (permute == "predictor") { + if (is.null(comparison)) { + if (family != "multinom") { + n_coefs <- length(fit$base$coefficients) + fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) + fit$larger <- fit$abs <- fit$lower + colnames(fit$lower) <- colnames(fit$larger) <- + colnames(fit$abs) <- names(fit$base$coefficients) + } else { + ncat <- if (large) { + length(stats::na.omit(unique(as.vector(unlist(matlist[[dep]]))))) + } else { + length(stats::na.omit(unique(as.vector(matlist[[dep]])))) } - xR <- residualise_predictor(xi, pred, main, - has_random = has_random, - rand_formula = rand_part) - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, - large, valid_list) - - gpu_res <- gpu_batch_ols_css(matlist = matlist_resid, - parsed = parsed, - directed = directed, - diag = diag, - groups = groups, - times = times, - baseline_fit = fit$base, - perm_var = xi) - fit$lower[, xi] <- gpu_res$lower[, xi] - fit$larger[, xi] <- gpu_res$larger[, xi] - fit$abs[, xi] <- gpu_res$abs[, xi] + n_coefs <- length(fit$base$coefficients) + fit$lower <- matrix(NA, nrow = 2 * (ncat - 1), ncol = n_coefs) + fit$larger <- fit$abs <- fit$lower + colnames(fit$lower) <- colnames(fit$larger) <- + colnames(fit$abs) <- names(fit$base$coefficients) + } + } else { + fit$lower <- fit$larger <- fit$abs <- + vector("list", length(comparison)) + names(fit$lower) <- names(fit$larger) <- + names(fit$abs) <- names(comparison) + for (k in seq_along(comparison)) { + n_coefs <- length(fit$base[[k]]$coefficients) + fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs) + fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] + colnames(fit$lower[[k]]) <- colnames(fit$larger[[k]]) <- + colnames(fit$abs[[k]]) <- names(fit$base[[k]]$coefficients) } } - } else { - old_plan <- setup_future_plan(strategy, ncores) - on.exit({ - future::plan(old_plan) - options(future.globals.maxSize = attr(old_plan, "old_maxSize")) - }, add = TRUE) + for (xi in main) { + test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] + if (!is.numeric(test_val)) { + manynet::snet_warn( + c("Cannot residualise the non-numeric predictor {.val {xi}}.", + i = "Skipping double semi-partialling for this predictor.")) + next + } + + xR <- residualise_predictor(xi, pred, main, + has_random = has_random, + rand_formula = rand_part) + + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, + large, valid_list) - if (permute == "outcome") { res <- run_permutations( times, QAPcssPermEst, - matlist. = matlist, - perm_var. = NULL, + matlist. = matlist_resid, + perm_var. = xi, directed. = directed, diag. = diag, mod. = mod, @@ -561,113 +602,18 @@ QAPcss <- function(formula, if (is.null(comparison)) { agg <- aggregate_perm_results(res, times) - fit$lower <- agg$lower - fit$larger <- agg$larger - fit$abs <- agg$abs + fit$lower[, xi] <- agg$lower + fit$larger[, xi] <- agg$larger + fit$abs[, xi] <- agg$abs } else { res_valid <- Filter(Negate(is.null), res) n_valid <- length(res_valid) - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) for (k in seq_along(comparison)) { cn <- names(comparison)[k] - fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } - - } else if (permute == "predictor") { - if (is.null(comparison)) { - if (family != "multinom") { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - } else { - ncat <- if (large) { - length(stats::na.omit(unique(as.vector(unlist(matlist[[dep]]))))) - } else { - length(stats::na.omit(unique(as.vector(matlist[[dep]])))) - } - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2 * (ncat - 1), ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - } - } else { - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) - for (k in seq_along(comparison)) { - n_coefs <- length(fit$base[[k]]$coefficients) - fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] - colnames(fit$lower[[k]]) <- colnames(fit$larger[[k]]) <- - colnames(fit$abs[[k]]) <- names(fit$base[[k]]$coefficients) - } - } - - for (xi in main) { - test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] - if (!is.numeric(test_val)) { - manynet::snet_warn( - c("Cannot residualise the non-numeric predictor {.val {xi}}.", - i = "Skipping double semi-partialling for this predictor.")) - next - } - - xR <- residualise_predictor(xi, pred, main, - has_random = has_random, - rand_formula = rand_part) - - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, - large, valid_list) - - res <- run_permutations( - times, QAPcssPermEst, - matlist. = matlist_resid, - perm_var. = xi, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, - family. = family, - estimator. = estimator, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference - ) - - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower[, xi] <- agg$lower - fit$larger[, xi] <- agg$larger - fit$abs[, xi] <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } + fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid + fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid + fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid } } } diff --git a/R/qap_engine.R b/R/qap_engine.R index ca1a0a5..3dbbc1a 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -26,8 +26,7 @@ QAPglm <- function(formula, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, use_robust_errors = FALSE, - less_mem = FALSE, - use_gpu = FALSE) { + less_mem = FALSE) { if (!is.null(seed)) set.seed(seed) @@ -139,73 +138,90 @@ QAPglm <- function(formula, "with one predictor there is nothing to residualise it against.") } - # The GPU path is a shortcut, not a requirement, so an unmet condition falls - # back to the CPU permutation loop rather than aborting. `gpu_available()` - # covers the two conditions the user cannot see from the call: whether - # {torch} is installed, and whether CUDA is reachable. - use_gpu <- use_gpu && family == "gaussian" && !has_random && !use_fixest && - is.null(comparison) && !large - if (use_gpu && !gpu_available()) { - manynet::snet_info( - "No CUDA device is reachable, so using the CPU permutation path.") - use_gpu <- FALSE - } + old_plan <- setup_future_plan(strategy, ncores) + on.exit({ + future::plan(old_plan) + options(future.globals.maxSize = attr(old_plan, "old_maxSize")) + }, add = TRUE) + + if (permute == "outcome") { + res <- run_permutations( + times, QAPglmPermEst, + matlist. = matlist, + perm_var. = NULL, + directed. = directed, + diag. = diag, + mod. = mod, + groups. = groups, + fit. = if (is.null(comparison)) fit$base else fit$base, + family. = family, + estimator. = estimator, + use_fixest. = use_fixest, + fixest_se_cluster. = fixest_se_cluster, + use_robust_errors. = use_robust_errors, + has_random. = has_random, + main_vars. = main, + data_vars. = data_vars, + parsed. = parsed, + comp. = comparison, + reference. = reference + ) + + if (is.null(comparison)) { + agg <- aggregate_perm_results(res, times) + fit$lower <- agg$lower + fit$larger <- agg$larger + fit$abs <- agg$abs + } else { + res_valid <- Filter(Negate(is.null), res) + n_valid <- length(res_valid) + fit$lower <- fit$larger <- fit$abs <- + vector("list", length(comparison)) + names(fit$lower) <- names(comparison) + names(fit$larger) <- names(comparison) + names(fit$abs) <- names(comparison) + resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) + for (k in seq_along(comparison)) { + cn <- names(comparison)[k] + fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid + fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid + fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid + } + } - if (use_gpu) { - - if (permute == "outcome") { - gpu_res <- gpu_batch_ols(matlist = matlist, - parsed = parsed, - directed = directed, - diag = diag, - groups = groups, - times = times, - baseline_fit = fit$base, - perm_var = NULL) - fit$lower <- gpu_res$lower - fit$larger <- gpu_res$larger - fit$abs <- gpu_res$abs - - } else if (permute == "predictor") { + } else if (permute == "predictor") { + if (is.null(comparison)) { n_coefs <- length(fit$base$coefficients) fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, dimnames = list(c("perm_coefs", "perm_t"), names(fit$base$coefficients))) fit$larger <- fit$abs <- fit$lower - - for (xi in main) { - xR <- residualise_predictor(xi, pred, main, - has_random = has_random, - rand_formula = rand_part) - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) - - gpu_res <- gpu_batch_ols(matlist = matlist_resid, - parsed = parsed, - directed = directed, - diag = diag, - groups = groups, - times = times, - baseline_fit = fit$base, - perm_var = xi) - fit$lower[, xi] <- gpu_res$lower[, xi] - fit$larger[, xi] <- gpu_res$larger[, xi] - fit$abs[, xi] <- gpu_res$abs[, xi] + } else { + fit$lower <- fit$larger <- fit$abs <- + vector("list", length(comparison)) + names(fit$lower) <- names(fit$larger) <- + names(fit$abs) <- names(comparison) + for (k in seq_along(comparison)) { + n_coefs <- length(fit$base[[k]]$coefficients) + fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs, + dimnames = list(c("perm_coefs", "perm_t"), + names(fit$base[[k]]$coefficients))) + fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] } } - } else { - old_plan <- setup_future_plan(strategy, ncores) - on.exit({ - future::plan(old_plan) - options(future.globals.maxSize = attr(old_plan, "old_maxSize")) - }, add = TRUE) + for (xi in main) { + xR <- residualise_predictor(xi, pred, main, + has_random = has_random, + rand_formula = rand_part) + + matlist_resid <- matlist + matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) - if (permute == "outcome") { res <- run_permutations( times, QAPglmPermEst, - matlist. = matlist, - perm_var. = NULL, + matlist. = matlist_resid, + perm_var. = xi, directed. = directed, diag. = diag, mod. = mod, @@ -226,92 +242,18 @@ QAPglm <- function(formula, if (is.null(comparison)) { agg <- aggregate_perm_results(res, times) - fit$lower <- agg$lower - fit$larger <- agg$larger - fit$abs <- agg$abs + fit$lower[, xi] <- agg$lower + fit$larger[, xi] <- agg$larger + fit$abs[, xi] <- agg$abs } else { res_valid <- Filter(Negate(is.null), res) n_valid <- length(res_valid) - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(comparison) - names(fit$larger) <- names(comparison) - names(fit$abs) <- names(comparison) resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) for (k in seq_along(comparison)) { cn <- names(comparison)[k] - fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } - - } else if (permute == "predictor") { - if (is.null(comparison)) { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, - dimnames = list(c("perm_coefs", "perm_t"), - names(fit$base$coefficients))) - fit$larger <- fit$abs <- fit$lower - } else { - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) - for (k in seq_along(comparison)) { - n_coefs <- length(fit$base[[k]]$coefficients) - fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs, - dimnames = list(c("perm_coefs", "perm_t"), - names(fit$base[[k]]$coefficients))) - fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] - } - } - - for (xi in main) { - xR <- residualise_predictor(xi, pred, main, - has_random = has_random, - rand_formula = rand_part) - - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) - - res <- run_permutations( - times, QAPglmPermEst, - matlist. = matlist_resid, - perm_var. = xi, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, - family. = family, - estimator. = estimator, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference - ) - - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower[, xi] <- agg$lower - fit$larger[, xi] <- agg$larger - fit$abs[, xi] <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } + fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid + fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid + fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid } } } diff --git a/R/qap_gpu.R b/R/qap_gpu.R deleted file mode 100644 index 69a9658..0000000 --- a/R/qap_gpu.R +++ /dev/null @@ -1,318 +0,0 @@ -# GPU batch OLS -------------------------------------------------------------- -# -# Optional torch-based batch OLS for gaussian QAP permutations. Internal; -# triggered by control$use_gpu = TRUE in net_regression(). - -#' @keywords internal -#' @noRd -gpu_batch_ols <- function(matlist, parsed, directed, diag, groups, times, - baseline_fit, perm_var = NULL, - batch_size = 500, device = "cuda") { - - thisRequires("torch", "for GPU acceleration") - - if (device == "cuda" && !torch::cuda_is_available()) { - manynet::snet_info("CUDA is not available, so falling back to CPU {.pkg torch}.") - device <- "cpu" - } - - dep <- parsed$dependent - main <- parsed$main - - pred0 <- make_qap_data(y = matlist[[dep]], - x = matlist[main], - g = groups, - diag = diag, - directed = directed, - net = 1, - perm = FALSE, - xi = NULL) - - y_vec <- pred0$yv - X_mat <- cbind(1, as.matrix(pred0[, main, drop = FALSE])) - n_obs <- nrow(X_mat) - p <- ncol(X_mat) - - base_coefs <- baseline_fit$coefficients - base_t <- baseline_fit$t - - lower_sum <- rep(0, length(base_coefs) * 2) - larger_sum <- rep(0, length(base_coefs) * 2) - abs_sum <- rep(0, length(base_coefs) * 2) - - dim_out <- c(2, length(base_coefs)) - - if (is.null(perm_var)) { - X_t <- torch::torch_tensor(X_mat, dtype = torch::torch_float64(), - device = device) - XtX <- torch::torch_mm(torch::torch_t(X_t), X_t) - XtXinv <- torch::torch_inverse(XtX) - M <- torch::torch_mm(XtXinv, torch::torch_t(X_t)) - XtXinv_diag <- torch::torch_diag(XtXinv) - - reps_done <- 0 - while (reps_done < times) { - current_batch <- min(batch_size, times - reps_done) - - Y_batch <- matrix(NA_real_, nrow = n_obs, ncol = current_batch) - for (j in seq_len(current_batch)) { - y_perm <- RMPerm(matlist[[dep]], groups) - perm_pred <- make_qap_data(y = y_perm, - x = matlist[main], - g = groups, - diag = diag, - directed = directed, - net = 1, - perm = FALSE, - xi = NULL) - Y_batch[, j] <- perm_pred$yv - } - - Y_t <- torch::torch_tensor(Y_batch, dtype = torch::torch_float64(), - device = device) - - B <- torch::torch_mm(M, Y_t) - E <- Y_t - torch::torch_mm(X_t, B) - MSE <- torch::torch_sum(E^2, dim = 1) / (n_obs - p) - SE <- torch::torch_sqrt(torch::torch_ger(XtXinv_diag, MSE)) - T_vals <- B / SE - - B_cpu <- as.matrix(B$cpu()) - T_cpu <- as.matrix(T_vals$cpu()) - - for (j in seq_len(current_batch)) { - pres <- rbind(B_cpu[, j], T_cpu[, j]) - bres <- rbind(base_coefs, base_t) - lower_sum <- lower_sum + as.vector(pres <= bres) - larger_sum <- larger_sum + as.vector(pres >= bres) - abs_sum <- abs_sum + as.vector(abs(pres) >= abs(bres)) - } - - reps_done <- reps_done + current_batch - } - - } else { - y_t <- torch::torch_tensor(matrix(y_vec, ncol = 1), - dtype = torch::torch_float64(), - device = device) - - reps_done <- 0 - while (reps_done < times) { - current_batch <- min(batch_size, times - reps_done) - - B_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) - T_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) - - for (j in seq_len(current_batch)) { - d_perm <- matlist - d_perm[[perm_var]] <- RMPerm(d_perm[[perm_var]], groups) - - perm_pred <- make_qap_data(y = d_perm[[dep]], - x = d_perm[main], - g = groups, - diag = diag, - directed = directed, - net = 1, - perm = FALSE, - xi = NULL) - X_perm <- cbind(1, as.matrix(perm_pred[, main, drop = FALSE])) - Xp_t <- torch::torch_tensor(X_perm, dtype = torch::torch_float64(), - device = device) - - XpXp <- torch::torch_mm(torch::torch_t(Xp_t), Xp_t) - XpXp_inv <- torch::torch_inverse(XpXp) - b_perm <- torch::torch_mm(XpXp_inv, - torch::torch_mm(torch::torch_t(Xp_t), y_t)) - e_perm <- y_t - torch::torch_mm(Xp_t, b_perm) - mse_perm <- (torch::torch_sum(e_perm^2) / (n_obs - p))$item() - XpXp_diag <- as.numeric(torch::torch_diag(XpXp_inv)$cpu()) - se_perm <- sqrt(XpXp_diag * mse_perm) - - b_cpu <- as.numeric(b_perm$cpu()) - t_cpu <- b_cpu / se_perm - - B_batch[, j] <- b_cpu - T_batch[, j] <- t_cpu - } - - for (j in seq_len(current_batch)) { - pres <- rbind(B_batch[, j], T_batch[, j]) - bres <- rbind(base_coefs, base_t) - lower_sum <- lower_sum + as.vector(pres <= bres) - larger_sum <- larger_sum + as.vector(pres >= bres) - abs_sum <- abs_sum + as.vector(abs(pres) >= abs(bres)) - } - - reps_done <- reps_done + current_batch - } - } - - list( - lower = matrix(lower_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))), - larger = matrix(larger_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))), - abs = matrix(abs_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))) - ) -} - - -#' @keywords internal -#' @noRd -gpu_batch_ols_css <- function(matlist, parsed, directed, diag, groups, times, - baseline_fit, perm_var = NULL, - batch_size = 500, device = "cuda") { - - thisRequires("torch", "for GPU acceleration") - - if (device == "cuda" && !torch::cuda_is_available()) { - manynet::snet_info("CUDA is not available, so falling back to CPU {.pkg torch}.") - device <- "cpu" - } - - dep <- parsed$dependent - main <- parsed$main - data_vars <- parsed$all_data_vars - - x_list <- lapply(data_vars, function(v) matlist[[v]]) - names(x_list) <- data_vars - cssd <- make_css_data(y = matlist[[dep]], x = x_list, - nets = 1, diag = diag, directed = directed) - pred0 <- cssd$pred - - y_vec <- pred0$yv - X_mat <- cbind(1, as.matrix(pred0[, main, drop = FALSE])) - n_obs <- nrow(X_mat) - p <- ncol(X_mat) - - base_coefs <- baseline_fit$coefficients - base_t <- baseline_fit$t - - lower_sum <- rep(0, length(base_coefs) * 2) - larger_sum <- rep(0, length(base_coefs) * 2) - abs_sum <- rep(0, length(base_coefs) * 2) - - dim_out <- c(2, length(base_coefs)) - - build_css_pred <- function(d) { - xl <- lapply(data_vars, function(v) d[[v]]) - names(xl) <- data_vars - make_css_data(y = d[[dep]], x = xl, - nets = 1, diag = diag, directed = directed)$pred - } - - if (is.null(perm_var)) { - X_t <- torch::torch_tensor(X_mat, dtype = torch::torch_float64(), - device = device) - XtX <- torch::torch_mm(torch::torch_t(X_t), X_t) - XtXinv <- torch::torch_inverse(XtX) - M <- torch::torch_mm(XtXinv, torch::torch_t(X_t)) - XtXinv_diag <- torch::torch_diag(XtXinv) - - reps_done <- 0 - while (reps_done < times) { - current_batch <- min(batch_size, times - reps_done) - - Y_batch <- matrix(NA_real_, nrow = n_obs, ncol = current_batch) - for (j in seq_len(current_batch)) { - d_perm <- matlist - d_perm[[dep]] <- RMPerm(d_perm[[dep]], groups, CSS = TRUE) - perm_pred <- build_css_pred(d_perm) - Y_batch[, j] <- perm_pred$yv - } - - Y_t <- torch::torch_tensor(Y_batch, dtype = torch::torch_float64(), - device = device) - - B <- torch::torch_mm(M, Y_t) - E <- Y_t - torch::torch_mm(X_t, B) - MSE <- torch::torch_sum(E^2, dim = 1) / (n_obs - p) - SE <- torch::torch_sqrt(torch::torch_ger(XtXinv_diag, MSE)) - T_vals <- B / SE - - B_cpu <- as.matrix(B$cpu()) - T_cpu <- as.matrix(T_vals$cpu()) - - for (j in seq_len(current_batch)) { - pres <- rbind(B_cpu[, j], T_cpu[, j]) - bres <- rbind(base_coefs, base_t) - lower_sum <- lower_sum + as.vector(pres <= bres) - larger_sum <- larger_sum + as.vector(pres >= bres) - abs_sum <- abs_sum + as.vector(abs(pres) >= abs(bres)) - } - - reps_done <- reps_done + current_batch - } - - } else { - y_t <- torch::torch_tensor(matrix(y_vec, ncol = 1), - dtype = torch::torch_float64(), - device = device) - - reps_done <- 0 - while (reps_done < times) { - current_batch <- min(batch_size, times - reps_done) - - B_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) - T_batch <- matrix(NA_real_, nrow = p, ncol = current_batch) - - for (j in seq_len(current_batch)) { - d_perm <- matlist - d_perm[[perm_var]] <- RMPerm(d_perm[[perm_var]], groups, CSS = TRUE) - - perm_pred <- build_css_pred(d_perm) - X_perm <- cbind(1, as.matrix(perm_pred[, main, drop = FALSE])) - Xp_t <- torch::torch_tensor(X_perm, dtype = torch::torch_float64(), - device = device) - - XpXp <- torch::torch_mm(torch::torch_t(Xp_t), Xp_t) - XpXp_inv <- torch::torch_inverse(XpXp) - b_perm <- torch::torch_mm(XpXp_inv, - torch::torch_mm(torch::torch_t(Xp_t), y_t)) - e_perm <- y_t - torch::torch_mm(Xp_t, b_perm) - mse_perm <- (torch::torch_sum(e_perm^2) / (n_obs - p))$item() - XpXp_diag <- as.numeric(torch::torch_diag(XpXp_inv)$cpu()) - se_perm <- sqrt(XpXp_diag * mse_perm) - - b_cpu <- as.numeric(b_perm$cpu()) - t_cpu <- b_cpu / se_perm - - B_batch[, j] <- b_cpu - T_batch[, j] <- t_cpu - } - - for (j in seq_len(current_batch)) { - pres <- rbind(B_batch[, j], T_batch[, j]) - bres <- rbind(base_coefs, base_t) - lower_sum <- lower_sum + as.vector(pres <= bres) - larger_sum <- larger_sum + as.vector(pres >= bres) - abs_sum <- abs_sum + as.vector(abs(pres) >= abs(bres)) - } - - reps_done <- reps_done + current_batch - } - } - - list( - lower = matrix(lower_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))), - larger = matrix(larger_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))), - abs = matrix(abs_sum / times, nrow = dim_out[1], ncol = dim_out[2], - dimnames = list(NULL, names(base_coefs))) - ) -} - - -#' @keywords internal -#' @noRd -gpu_available <- function() { - if (!requireNamespace("torch", quietly = TRUE)) return(FALSE) - # {torch} installs as an R package before its Lantern backend is downloaded, - # so `cuda_is_available()` throws rather than returning FALSE on a machine - # that has the package but not the runtime. That is the state of a CI runner - # that installed Suggests, and it must read as "no GPU", not as an error. - isTRUE(tryCatch(torch::cuda_is_available(), - error = function(e) FALSE, warning = function(w) FALSE)) -} diff --git a/man/regression.Rd b/man/regression.Rd index 7628cec..fd23bf5 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -58,7 +58,6 @@ to i. Read from \code{.data} unless given, and reported when read. \item \code{reference}, \code{comparison}: multinomial / pairwise-comparison options. \item \code{random_intercept_nets} / \verb{_sender} / \verb{_receiver}: lme4-style REs. \item \code{less_mem}: drop the baseline model object from the return. -\item \code{use_gpu}: torch-based batch OLS (gaussian only). }} \item{x}{A \code{net_regression} object.} @@ -93,7 +92,6 @@ multinomial families; and \code{"outcome"}; \item random intercepts (lme4 / glmmTMB) and fixed effects (fixest); \item robust (HC3) standard errors; -\item optional torch-based batch OLS on the GPU; \item lists of networks, in which graphs that are missing any predictor are dropped with a warning and the remaining networks are pooled. } diff --git a/tests/testthat/test-qap_estimators.R b/tests/testthat/test-qap_estimators.R index fe9d8c4..ee27f10 100644 --- a/tests/testthat/test-qap_estimators.R +++ b/tests/testthat/test-qap_estimators.R @@ -239,16 +239,3 @@ test_that("combining fixed and random effects warns", { "random effects") }) - -# ---- GPU ------------------------------------------------------------------- - -test_that("use_gpu falls back to the CPU rather than aborting", { - # The GPU path is a shortcut, so an unmet condition must not stop the run. - gpu <- suppressMessages( - net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1, use_gpu = TRUE))) - cpu <- net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1)) - expect_qap_shape(gpu, COEFS3) - if (!gpu_available()) expect_equal(gpu$lower, cpu$lower) -}) diff --git a/tests/testthat/test-qap_reporting.R b/tests/testthat/test-qap_reporting.R index 5af0278..2e2fb8b 100644 --- a/tests/testthat/test-qap_reporting.R +++ b/tests/testthat/test-qap_reporting.R @@ -22,9 +22,6 @@ test_that("every reporting message renders", { net_regression(FORM, g, times = 5, control = list(seed = 1)))) expect_no_error(verbosely( net_regression(weight ~ ego(Age), g, times = 5, control = list(seed = 1)))) - expect_no_error(verbosely( - net_regression(FORM, g, times = 5, - control = list(seed = 1, use_gpu = TRUE)))) }) test_that("a family resolved from the outcome is reported", { @@ -58,13 +55,6 @@ test_that("the fallback to permuting the outcome is reported", { "residualise") }) -test_that("the GPU falling back to the CPU is reported", { - skip_if(gpu_available(), "a CUDA device is present, so there is no fallback") - expect_message( - verbosely(net_regression(FORM, qap_net_gaussian(n = 15), times = 5, - control = list(seed = 1, use_gpu = TRUE))), - "CPU") -}) test_that("the model advice on homophily terms is reported", { expect_message( From 1113795d2d6554396ab7296e23c61abee589cbc9 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:24:32 +0200 Subject: [PATCH 03/12] Branched off the GMM estimator Moved to feature/gmm-estimator. Reinstate with `git revert` of this commit. The GMM path warned "the covariance matrix of the coefficients is singular" on every family it declares, on well-conditioned data with real signal, so its standard errors are not what they claim to be. It returns coefficients, so a test that only checks the shape passes. Removes R/qap_gmm.R, the 68-line branch in fit_qap_model(), and the `estimator` control, which had one remaining value. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 1 - R/model_regression.R | 9 +-- R/qap_css.R | 24 ++---- R/qap_engine.R | 23 ++---- R/qap_gmm.R | 111 --------------------------- R/qap_utils.R | 73 +----------------- man/regression.Rd | 2 - tests/testthat/test-qap_control.R | 2 +- tests/testthat/test-qap_estimators.R | 35 +-------- 9 files changed, 17 insertions(+), 263 deletions(-) delete mode 100644 R/qap_gmm.R diff --git a/DESCRIPTION b/DESCRIPTION index 1eeb735..fe01816 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,7 +41,6 @@ Suggests: lme4, nnet, fixest, - gmm, MASS, pscl, testthat (>= 3.0.0), diff --git a/R/model_regression.R b/R/model_regression.R index 36293f6..16ab4da 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -49,8 +49,6 @@ #' - `family`: `"auto"` (default; gaussian for weighted networks, binomial #' for binary), `"gaussian"`, `"binomial"`, `"poisson"`, `"negbin"`, #' `"zip"`, or `"multinom"`. -#' - `estimator`: `"standard"` (default) or `"gmm"` (binomial/poisson/ -#' negbin/zip). #' - `directed`: logical, whether a tie from i to j differs from one from j #' to i. Read from `.data` unless given, and reported when read. #' - `diag`: logical, include loops (default auto-detected). @@ -140,7 +138,6 @@ net_regression <- function(formula, directed = ctrl$directed, diag = ctrl$diag, permute = ctrl$permute, - estimator = ctrl$estimator, times = times, seed = ctrl$seed, groups = ctrl$groups, @@ -220,7 +217,6 @@ net_regression <- function(formula, permute = .permute_schemes(), strategy = "sequential", family = "auto", - estimator = "standard", directed = NULL, diag = NULL, seed = NULL, @@ -456,8 +452,6 @@ print.net_regression <- function(x, ..., } else { cat("\nGeneralized Linear Mixed Network Model fit by REML\n") } - if (!is.null(x$estimator) && x$estimator == "gmm") - cat("\nEstimator: Generalized Method-of-Moments.") if (!is.null(x$theta)) cat("\nNegative binomial dispersion (theta):", format(round(x$theta, 4))) if (!is.null(x$zi_coefficients)) { @@ -534,8 +528,7 @@ print.net_regression <- function(x, ..., cat("--------------\n") } - if (!is.null(x$simple_fit) && !is.null(x$estimator) && - x$estimator != "gmm") { + if (!is.null(x$simple_fit)) { cat("\nAIC:", format(stats::AIC(x$simple_fit))) cat("\nBIC:", format(stats::BIC(x$simple_fit))) } diff --git a/R/qap_css.R b/R/qap_css.R index 0abc8b0..024b4d6 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -93,7 +93,6 @@ QAPcssPermEst <- function(i, groups., fit., family., - estimator., use_fixest., fixest_se_cluster., use_robust_errors., @@ -208,7 +207,6 @@ QAPcssPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., - estimator = estimator., use_fixest = use_fixest., fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., @@ -238,7 +236,6 @@ QAPcssPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = predK, family = family., - estimator = estimator., use_fixest = use_fixest., fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., @@ -323,7 +320,6 @@ QAPcss <- function(formula, strategy = "sequential", ncores = NULL, family = "gaussian", - estimator = "standard", groups = NULL, fixest_se_cluster = NULL, reference = NULL, @@ -451,7 +447,6 @@ QAPcss <- function(formula, fit$base <- fit_qap_model(mod = mod, pred = pred, family = family, - estimator = estimator, use_fixest = use_fixest, fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, @@ -467,7 +462,6 @@ QAPcss <- function(formula, fit$base[[k]] <- fit_qap_model(mod = mod, pred = predK, family = family, - estimator = estimator, use_fixest = use_fixest, fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, @@ -494,7 +488,6 @@ QAPcss <- function(formula, groups. = groups, fit. = if (is.null(comparison)) fit$base else fit$base, family. = family, - estimator. = estimator, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, @@ -588,7 +581,6 @@ QAPcss <- function(formula, groups. = groups, fit. = if (is.null(comparison)) fit$base else fit$base, family. = family, - estimator. = estimator, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, @@ -620,15 +612,11 @@ QAPcss <- function(formula, } if (family == "binomial" && is.null(comparison)) { - bm <- fit$base$base_model - if (!inherits(bm, "gmm")) { - predicted <- stats::fitted(bm) - actual <- pred[[dep]] - fit$confusion_matrix <- probabilistic_confusion_matrix( - actual = actual, predicted_prob = predicted, - n_draws = 1000, seed = seed - ) - } + fit$confusion_matrix <- probabilistic_confusion_matrix( + actual = pred[[dep]], + predicted_prob = stats::fitted(fit$base$base_model), + n_draws = 1000, seed = seed + ) } fit$permute <- permute @@ -644,7 +632,6 @@ QAPcss <- function(formula, perceiver = rip, nets = rin) fit$robust_se <- use_robust_errors - fit$estimator <- estimator if (is.null(comparison) && !is.null(fit$base$theta)) fit$theta <- fit$base$theta @@ -674,7 +661,6 @@ print.QAPCSS <- function(x, ...) { cat("The reference group was", format(paste0(x$reference, ".")), "\n") } - if (!is.null(x$estimator) && x$estimator == "gmm") cat("Estimator: Generalized Method-of-Moments.\n") if (!is.null(x$theta)) cat("Negative binomial dispersion (theta):", format(round(x$theta, 4)), "\n") diff --git a/R/qap_engine.R b/R/qap_engine.R index 3dbbc1a..1e5c935 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -13,7 +13,6 @@ QAPglm <- function(formula, directed = TRUE, diag = FALSE, permute = "predictor", - estimator = "standard", times = 1000, seed = NULL, groups = NULL, @@ -100,7 +99,6 @@ QAPglm <- function(formula, fit$base <- fit_qap_model(mod = mod, pred = pred, family = family, - estimator = estimator, use_fixest = use_fixest, fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, @@ -116,7 +114,6 @@ QAPglm <- function(formula, fit$base[[k]] <- fit_qap_model(mod = mod, pred = predK, family = family, - estimator = estimator, use_fixest = use_fixest, fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, @@ -155,7 +152,6 @@ QAPglm <- function(formula, groups. = groups, fit. = if (is.null(comparison)) fit$base else fit$base, family. = family, - estimator. = estimator, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, @@ -228,7 +224,6 @@ QAPglm <- function(formula, groups. = groups, fit. = if (is.null(comparison)) fit$base else fit$base, family. = family, - estimator. = estimator, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, @@ -280,15 +275,11 @@ QAPglm <- function(formula, } if (family == "binomial" && is.null(comparison)) { - bm <- fit$base$base_model - if (!inherits(bm, "gmm")) { - predicted <- stats::fitted(bm) - actual <- pred[[dep]] - fit$confusion_matrix <- probabilistic_confusion_matrix( - actual = actual, predicted_prob = predicted, - n_draws = 1000, seed = seed - ) - } + fit$confusion_matrix <- probabilistic_confusion_matrix( + actual = pred[[dep]], + predicted_prob = stats::fitted(fit$base$base_model), + n_draws = 1000, seed = seed + ) } fit$permute <- permute @@ -298,7 +289,6 @@ QAPglm <- function(formula, fit$times <- times fit$groups <- unique(unlist(groups)) fit$robust_se <- use_robust_errors - fit$estimator <- estimator fit$comp <- comparison fit$reference <- reference fit$pred <- pred @@ -324,7 +314,6 @@ QAPglmPermEst <- function(i, groups., fit., family., - estimator., use_fixest., fixest_se_cluster., use_robust_errors., @@ -393,7 +382,6 @@ QAPglmPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., - estimator = estimator., use_fixest = use_fixest., fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., @@ -423,7 +411,6 @@ QAPglmPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = predK, family = family., - estimator = estimator., use_fixest = use_fixest., fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., diff --git a/R/qap_gmm.R b/R/qap_gmm.R deleted file mode 100644 index 6863472..0000000 --- a/R/qap_gmm.R +++ /dev/null @@ -1,111 +0,0 @@ -# GMM moment conditions and residuals ---------------------------------------- -# -# Auxiliary estimators for the `estimator = "gmm"` path in fit_qap_model(). -# Internal; ported from MrQAP. - -#' @keywords internal -#' @noRd -poisson_moments <- function(theta, data) { - Y <- as.numeric(data$y) - X <- data.matrix(data$x) - lambda_hat <- exp(X %*% theta) - residuals <- as.vector(Y - lambda_hat) - g <- residuals * X - return(g) -} - -#' @keywords internal -#' @noRd -logit_moments <- function(theta, data) { - Y <- data$y - X <- data$x - prob <- 1 / (1 + exp(-1 * (X %*% theta))) - residuals <- as.vector(Y - prob) - g <- residuals * X - return(g) -} - -#' @keywords internal -#' @noRd -logit_resid <- function(gmmo) { - Y <- gmmo$dat$y - X <- gmmo$dat$x - prob <- 1 / (1 + exp(-1 * (X %*% gmmo$coefficients))) - residuals <- as.vector(Y - prob) - return(residuals) -} - -#' @keywords internal -#' @noRd -poisson_resid <- function(gmmo) { - Y <- gmmo$dat$y - X <- gmmo$dat$x - lambda_hat <- exp(X %*% gmmo$coefficients) - residuals <- as.vector(Y - lambda_hat) - return(residuals) -} - -#' @keywords internal -#' @noRd -negbin_moments <- function(theta, data) { - Y <- as.numeric(data$y) - X <- data.matrix(data$x) - p <- ncol(X) - beta <- theta[1:p] - alpha <- exp(theta[p + 1]) - - mu <- as.vector(exp(X %*% beta)) - resid <- Y - mu - V <- mu + alpha * mu^2 - - g1 <- (resid / V) * X - g2 <- (resid^2 / V) - 1 - - cbind(g1, g2) -} - -#' @keywords internal -#' @noRd -negbin_resid <- function(gmmo) { - Y <- as.numeric(gmmo$dat$y) - X <- data.matrix(gmmo$dat$x) - p <- ncol(X) - beta <- gmmo$coefficients[1:p] - mu <- as.vector(exp(X %*% beta)) - as.vector(Y - mu) -} - -#' @keywords internal -#' @noRd -zip_moments <- function(theta, data) { - Y <- as.numeric(data$y) - X <- data.matrix(data$x) - p <- ncol(X) - beta <- theta[1:p] - pi_z <- 1 / (1 + exp(-theta[p + 1])) - - lambda <- as.vector(exp(X %*% beta)) - mu <- (1 - pi_z) * lambda - resid <- Y - mu - - g1 <- resid * X - - p0 <- pi_z + (1 - pi_z) * exp(-lambda) - is_zero <- as.numeric(Y == 0) - g2 <- is_zero - p0 - - cbind(g1, g2) -} - -#' @keywords internal -#' @noRd -zip_resid <- function(gmmo) { - Y <- as.numeric(gmmo$dat$y) - X <- data.matrix(gmmo$dat$x) - p <- ncol(X) - beta <- gmmo$coefficients[1:p] - pi_z <- 1 / (1 + exp(-gmmo$coefficients[p + 1])) - lambda <- as.vector(exp(X %*% beta)) - mu <- (1 - pi_z) * lambda - as.vector(Y - mu) -} diff --git a/R/qap_utils.R b/R/qap_utils.R index 837d06a..aeb6e02 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -302,7 +302,7 @@ HC3 <- function(X, e) { # The predictor names carry spaces ("ego Age"), so the model formula quotes them # and several fitters hand the backticks back in the coefficient names. Double # semi-partialling then looks a column up by the unquoted name and fails with a -# subscript error. The inner function has many early returns, one per estimator, +# subscript error. The inner function has many early returns, one per family, # so the names are cleaned here, where every path passes through exactly once. #' @keywords internal #' @noRd @@ -319,7 +319,6 @@ fit_qap_model <- function(...) { #' @keywords internal #' @noRd .fit_qap_model <- function(mod, pred, family, - estimator = "standard", use_fixest = FALSE, fixest_se_cluster = NULL, use_robust_errors = FALSE, @@ -344,76 +343,8 @@ fit_qap_model <- function(...) { return(fit) } - if (estimator == "gmm") { - thisRequires("gmm", "for GMM estimation") - y_vec <- pred[[dep_var]] - x_mat <- cbind(1, as.matrix(pred[, main_vars, drop = FALSE])) - - gmm_args <- list( - x = list(y = y_vec, x = x_mat), - t0 = stats::rnorm(nx + 1), - wmatrix = "optimal", vcov = "MDS", - optfct = "nlminb", - control = list(eval.max = 10000) - ) - - has_extra_param <- FALSE - - if (family == "binomial") { - gmm_args$g <- logit_moments - base_model <- do.call(gmm::gmm, gmm_args) - resid <- logit_resid(base_model) - } else if (family == "poisson") { - gmm_args$g <- poisson_moments - base_model <- do.call(gmm::gmm, gmm_args) - resid <- poisson_resid(base_model) - } else if (family == "negbin") { - gmm_args$g <- negbin_moments - gmm_args$t0 <- stats::rnorm(nx + 2) - base_model <- do.call(gmm::gmm, gmm_args) - resid <- negbin_resid(base_model) - has_extra_param <- TRUE - } else if (family == "zip") { - gmm_args$g <- zip_moments - gmm_args$t0 <- stats::rnorm(nx + 2) - base_model <- do.call(gmm::gmm, gmm_args) - resid <- zip_resid(base_model) - has_extra_param <- TRUE - } else { - manynet::snet_abort( - c("The GMM estimator is not available for the {.val {family}} family.", - i = "It is available for the binomial, poisson, negbin, and zip families.")) - } - - all_coefs <- base_model$coefficients - if (!use_robust_errors) { - all_t <- summary(base_model)$coefficients[, 3] - } - - if (has_extra_param) { - fit$coefficients <- all_coefs[1:(nx + 1)] - } else { - fit$coefficients <- all_coefs - } - names(fit$coefficients) <- c("(Intercept)", main_vars) - - if (use_robust_errors) { - xv <- as.matrix(pred[, main_vars, drop = FALSE]) - fit$t <- fit$coefficients / HC3(xv, resid) - } else { - if (has_extra_param) { - fit$t <- all_t[1:(nx + 1)] - } else { - fit$t <- all_t - } - } - names(fit$t) <- names(fit$coefficients) - fit$base_model <- base_model - fit$estimator <- "gmm" - return(fit) - } - if (family == "zip" && estimator == "standard") { + if (family == "zip") { if (has_random) { thisRequires("glmmTMB", "for mixed zero-inflated Poisson models") base_model <- glmmTMB::glmmTMB(mod, data = pred, diff --git a/man/regression.Rd b/man/regression.Rd index fd23bf5..9f5aa0f 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -47,8 +47,6 @@ reduces to \code{"outcome"} and says so. \item \code{family}: \code{"auto"} (default; gaussian for weighted networks, binomial for binary), \code{"gaussian"}, \code{"binomial"}, \code{"poisson"}, \code{"negbin"}, \code{"zip"}, or \code{"multinom"}. -\item \code{estimator}: \code{"standard"} (default) or \code{"gmm"} (binomial/poisson/ -negbin/zip). \item \code{directed}: logical, whether a tie from i to j differs from one from j to i. Read from \code{.data} unless given, and reported when read. \item \code{diag}: logical, include loops (default auto-detected). diff --git a/tests/testthat/test-qap_control.R b/tests/testthat/test-qap_control.R index e76593e..2e4f30c 100644 --- a/tests/testthat/test-qap_control.R +++ b/tests/testthat/test-qap_control.R @@ -32,7 +32,7 @@ test_that("a named control overrides only that default", { ctrl <- .resolve_control(list(family = "poisson")) expect_equal(ctrl$family, "poisson") expect_equal(ctrl$strategy, "sequential") - expect_equal(ctrl$estimator, "standard") + expect_null(ctrl$directed) }) test_that("permute takes only the two spellings it documents", { diff --git a/tests/testthat/test-qap_estimators.R b/tests/testthat/test-qap_estimators.R index ee27f10..b1a01b5 100644 --- a/tests/testthat/test-qap_estimators.R +++ b/tests/testthat/test-qap_estimators.R @@ -1,7 +1,7 @@ # Every estimator path in fit_qap_model() is selected by a combination of -# `family`, `estimator`, and the random and fixed effect flags. This file names -# each combination, so that a path with no test fails the build rather than -# going unnoticed. +# `family` and the random and fixed effect flags. This file names each +# combination, so that a path with no test fails the build rather than going +# unnoticed. # # Two things are asserted for each. First, the baseline coefficients equal those # of the equivalent standard fit on the same dyad-level data: the permutation @@ -110,35 +110,6 @@ test_that("zip baseline matches pscl::zeroinfl() and names its coefficients", { }) -# ---- GMM ------------------------------------------------------------------- - -test_that("the GMM estimator runs for each family it declares", { - skip_if_not_installed("gmm") - cases <- list( - list(form = FORM_B, net = qap_net_binary(), family = "binomial"), - list(form = FORM, net = qap_net_count(), family = "poisson"), - list(form = FORM, net = qap_net_count(), family = "negbin"), - list(form = FORM, net = qap_net_zip(), family = "zip") - ) - for (case in cases) { - fit <- suppressWarnings( - net_regression(case$form, case$net, times = 10, - control = list(seed = 1, family = case$family, - estimator = "gmm"))) - expect_qap_shape(fit, COEFS3) - expect_equal(fit$estimator, "gmm", info = case$family) - } -}) - -test_that("the GMM estimator rejects a family it cannot fit", { - skip_if_not_installed("gmm") - expect_error( - net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(family = "gaussian", estimator = "gmm")), - "binomial") -}) - - # ---- random effects -------------------------------------------------------- test_that("gaussian random intercepts match lme4::lmer() on the same dyads", { From ebb3d1420e839460f0d67780e82ecf67ac8a3a98 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:25:48 +0200 Subject: [PATCH 04/12] Branched off the glmmTMB mixed models Moved to feature/glmmtmb-mixed. Reinstate with `git revert` of this commit. Mixed negative binomial and mixed zero-inflated Poisson were the only two paths needing {glmmTMB}, which carries 62 recursive dependencies, an order of magnitude more than anything else in Suggests, and pulls {lme4} anyway. Its build must match {TMB}; the pair already fell out of step locally, so both paths went untested. Each combination now aborts and names the alternative. The standard zip path through {pscl} and the standard negbin path through {MASS} are unaffected. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 3 +-- R/model_regression.R | 2 +- R/qap_utils.R | 42 ++++++++---------------------------------- man/regression.Rd | 2 +- 4 files changed, 11 insertions(+), 38 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fe01816..7f58c85 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -43,8 +43,7 @@ Suggests: fixest, MASS, pscl, - testthat (>= 3.0.0), - glmmTMB + testthat (>= 3.0.0) Config/Needs/build: roxygen2, devtools diff --git a/R/model_regression.R b/R/model_regression.R index 16ab4da..b6b17e0 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -14,7 +14,7 @@ #' multinomial families; #' - two permutation schemes: `"predictor"` (Dekker's double semi-partialling) #' and `"outcome"`; -#' - random intercepts (lme4 / glmmTMB) and fixed effects (fixest); +#' - random intercepts (lme4) and fixed effects (fixest); #' - robust (HC3) standard errors; #' - lists of networks, in which graphs that are missing any predictor are #' dropped with a warning and the remaining networks are pooled. diff --git a/R/qap_utils.R b/R/qap_utils.R index aeb6e02..c2ee58c 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -346,20 +346,10 @@ fit_qap_model <- function(...) { if (family == "zip") { if (has_random) { - thisRequires("glmmTMB", "for mixed zero-inflated Poisson models") - base_model <- glmmTMB::glmmTMB(mod, data = pred, - family = stats::poisson(), - ziformula = ~1) - fit$coefficients <- glmmTMB::fixef(base_model)$cond - resid <- stats::residuals(base_model, type = "response") - fit$t <- summary(base_model)$coefficients$cond[, 3] - names(fit$t) <- names(fit$coefficients) - fit$zi_coefficients <- glmmTMB::fixef(base_model)$zi - fit$random.intercepts <- list() - re <- glmmTMB::ranef(base_model)$cond - for (rV in names(re)) { - fit$random.intercepts[[rV]] <- re[[rV]][, 1] - } + # The mixed variant needs {glmmTMB}, which is on feature/glmmtmb-mixed. + manynet::snet_abort( + c("Random intercepts are not available for the {.val zip} family.", + i = "Drop the random intercepts, or use {.val poisson}.")) } else { thisRequires("pscl", "for zero-inflated Poisson models") base_model <- pscl::zeroinfl(mod, data = pred, dist = "poisson") @@ -448,26 +438,10 @@ fit_qap_model <- function(...) { thisRequires("lme4", "for random effects") base_model <- lme4::lmer(mod, data = pred) } else if (family == "negbin") { - thisRequires("glmmTMB", "for mixed negative binomial models") - base_model <- glmmTMB::glmmTMB(mod, data = pred, - family = glmmTMB::nbinom2()) - fit$coefficients <- glmmTMB::fixef(base_model)$cond - resid <- stats::residuals(base_model, type = "response") - if (use_robust_errors) { - xv <- as.matrix(pred[, main_vars, drop = FALSE]) - fit$t <- fit$coefficients / HC3(xv, resid) - } else { - fit$t <- summary(base_model)$coefficients$cond[, 3] - names(fit$t) <- names(fit$coefficients) - } - fit$theta <- glmmTMB::sigma(base_model) - fit$random.intercepts <- list() - re <- glmmTMB::ranef(base_model)$cond - for (rV in names(re)) { - fit$random.intercepts[[rV]] <- re[[rV]][, 1] - } - fit$base_model <- base_model - return(fit) + # The mixed variant needs {glmmTMB}, which is on feature/glmmtmb-mixed. + manynet::snet_abort( + c("Random intercepts are not available for the {.val negbin} family.", + i = "Drop the random intercepts, or use {.val poisson}.")) } else { thisRequires("lme4", "for random effects") base_model <- lme4::glmer(mod, data = pred, family = family, diff --git a/man/regression.Rd b/man/regression.Rd index 9f5aa0f..4387e9f 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -88,7 +88,7 @@ and handed to a QAP engine (ported from MrQAP) that supports: multinomial families; \item two permutation schemes: \code{"predictor"} (Dekker's double semi-partialling) and \code{"outcome"}; -\item random intercepts (lme4 / glmmTMB) and fixed effects (fixest); +\item random intercepts (lme4) and fixed effects (fixest); \item robust (HC3) standard errors; \item lists of networks, in which graphs that are missing any predictor are dropped with a warning and the remaining networks are pooled. From 9e9c56618471f5062cecfdd3a8e998c91d8ecfe3 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:32:00 +0200 Subject: [PATCH 05/12] Branched off multinomial models and the comparison machinery Moved to feature/multinomial-comparison. Reinstate with `git revert` of this commit. `family = "multinom"` coerced the outcome to a factor, so a numeric tie weight produced one level per distinct weight and {nnet} refused. No call through `net_regression()` reached it. The `comparison` and `reference` controls belonged to it: they ran a pairwise branch returning a list of p-value matrices instead of one, so the shape of the result changed. That fork appeared at 21 points across the two engines, four of them reading `if (is.null(comparison)) fit$base else fit$base` -- the same value on both arms. Neither control was documented beyond a line, and neither was tested. This is the largest obstacle to merging the two engines, so it goes first. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 1 - R/model_regression.R | 58 +++---- R/qap_css.R | 373 +++++++++---------------------------------- R/qap_engine.R | 236 +++++++-------------------- R/qap_utils.R | 18 +-- man/regression.Rd | 7 +- 6 files changed, 160 insertions(+), 533 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7f58c85..af7480d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -39,7 +39,6 @@ Imports: reformulas Suggests: lme4, - nnet, fixest, MASS, pscl, diff --git a/R/model_regression.R b/R/model_regression.R index b6b17e0..f88bca9 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -10,8 +10,8 @@ #' object. Internally the response and predictors are packed into matrices #' and handed to a QAP engine (ported from MrQAP) that supports: #' -#' - gaussian, binomial, poisson, negbin, zero-inflated Poisson, and -#' multinomial families; +#' - gaussian, binomial, poisson, negative binomial, and zero-inflated +#' Poisson families; #' - two permutation schemes: `"predictor"` (Dekker's double semi-partialling) #' and `"outcome"`; #' - random intercepts (lme4) and fixed effects (fixest); @@ -48,14 +48,13 @@ #' - `strategy`: future plan, e.g. `"sequential"` (default), `"multisession"`. #' - `family`: `"auto"` (default; gaussian for weighted networks, binomial #' for binary), `"gaussian"`, `"binomial"`, `"poisson"`, `"negbin"`, -#' `"zip"`, or `"multinom"`. +#' or `"zip"`. #' - `directed`: logical, whether a tie from i to j differs from one from j #' to i. Read from `.data` unless given, and reported when read. #' - `diag`: logical, include loops (default auto-detected). #' - `seed`, `groups`, `ncores`: passed through to the engine. #' - `use_robust_errors`: HC3 standard errors. #' - `fixest_se_cluster`: cluster variable for fixest. -#' - `reference`, `comparison`: multinomial / pairwise-comparison options. #' - `random_intercept_nets` / `_sender` / `_receiver`: lme4-style REs. #' - `less_mem`: drop the baseline model object from the return. #' @return An object of class `net_regression` inheriting from either @@ -144,8 +143,6 @@ net_regression <- function(formula, strategy = ctrl$strategy, ncores = ctrl$ncores, fixest_se_cluster = ctrl$fixest_se_cluster, - comparison = ctrl$comparison, - reference = ctrl$reference, random_intercept_nets = ctrl$random_intercept_nets, random_intercept_sender = ctrl$random_intercept_sender, random_intercept_receiver = ctrl$random_intercept_receiver, @@ -153,7 +150,7 @@ net_regression <- function(formula, less_mem = ctrl$less_mem ) - if (user_requested_gaussian_binary && is.null(ctrl$comparison)) { + if (user_requested_gaussian_binary) { fit$confusion_matrix <- .lpm_confusion_matrix(fit, ctrl$seed) } @@ -224,8 +221,6 @@ net_regression <- function(formula, ncores = NULL, use_robust_errors = FALSE, fixest_se_cluster = NULL, - reference = NULL, - comparison = NULL, random_intercept_nets = FALSE, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, @@ -480,41 +475,32 @@ print.net_regression <- function(x, ..., format(paste0(.directed_label(x$directed), "."))) cat("\nModel family:", format(x$family)) - if (!is.null(x$comp)) { - for (k in seq_along(x$comp)) { - cat("\n\n--- Comparison:", names(x$comp)[k], "---") - cat("\n ", x$comp[[k]][1], "vs", x$comp[[k]][2]) - .print_glm_table(x$base[[k]], x$lower[[k]], x$larger[[k]], x$abs[[k]], - x$permute, print_b) - } - } else { - cat("\n\nCoefficients:\n") - if (print_b) { - cmat <- matrix(NA, nrow = length(x$coefficients), ncol = 5) - cmat[, 1] <- format(as.numeric(x$coefficients)) - cmat[, 2] <- format(exp(as.numeric(x$coefficients))) - cmat[, 3] <- format(x$lower[1, ]) - cmat[, 4] <- format(x$larger[1, ]) - cmat[, 5] <- format(x$abs[1, ]) - if (x$permute == "predictor") cmat[1, 3:5] <- "*" - colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=b)", "Pr(>=b)", "Pr(>=|b|)") - rownames(cmat) <- names(x$coefficients) - print.table(cmat) - cat("--------------\n") - } - + cat("\n\nCoefficients:\n") + if (print_b) { cmat <- matrix(NA, nrow = length(x$coefficients), ncol = 5) cmat[, 1] <- format(as.numeric(x$coefficients)) cmat[, 2] <- format(exp(as.numeric(x$coefficients))) - cmat[, 3] <- format(x$lower[2, ]) - cmat[, 4] <- format(x$larger[2, ]) - cmat[, 5] <- format(x$abs[2, ]) + cmat[, 3] <- format(x$lower[1, ]) + cmat[, 4] <- format(x$larger[1, ]) + cmat[, 5] <- format(x$abs[1, ]) if (x$permute == "predictor") cmat[1, 3:5] <- "*" - colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") + colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=b)", "Pr(>=b)", "Pr(>=|b|)") rownames(cmat) <- names(x$coefficients) print.table(cmat) + cat("--------------\n") } + cmat <- matrix(NA, nrow = length(x$coefficients), ncol = 5) + cmat[, 1] <- format(as.numeric(x$coefficients)) + cmat[, 2] <- format(exp(as.numeric(x$coefficients))) + cmat[, 3] <- format(x$lower[2, ]) + cmat[, 4] <- format(x$larger[2, ]) + cmat[, 5] <- format(x$abs[2, ]) + if (x$permute == "predictor") cmat[1, 3:5] <- "*" + colnames(cmat) <- c("Estimate", "Exp(b)", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") + rownames(cmat) <- names(x$coefficients) + print.table(cmat) + if (x$permute == "predictor") cat("\n* The intercept has no significance test when predictors are permuted.\n") diff --git a/R/qap_css.R b/R/qap_css.R index 024b4d6..a51c85d 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -99,9 +99,7 @@ QAPcssPermEst <- function(i, has_random., main_vars., data_vars., - parsed., - comp., - reference.) { + parsed.) { dep <- parsed.$dependent large <- is.list(matlist.[[dep]]) @@ -151,14 +149,7 @@ QAPcssPermEst <- function(i, names(pred)[names(pred) == "yv"] <- dep - if (family. != "multinom" && is.null(comp.)) { - y_ok <- length(stats::na.omit(unique(pred[[dep]]))) > 1 - } else { - y2_cat <- stats::na.omit(unique(pred[[dep]])) - y_present <- all(y_cat %in% y2_cat) - y_mult <- all(table(pred[[dep]]) > 2) - y_ok <- y_present && y_mult - } + y_ok <- length(stats::na.omit(unique(pred[[dep]]))) > 1 x_ok <- TRUE num_preds <- pred[, data_vars.[data_vars. %in% names(pred)], drop = FALSE] @@ -167,26 +158,6 @@ QAPcssPermEst <- function(i, x_ok <- all(sapply(num_preds, function(col) length(unique(col)) > 1)) } - if (nrow(pred) != 0 && x_ok && y_ok && !is.null(comp.)) { - for (k in seq_along(comp.)) { - pred2 <- pred[pred[[dep]] %in% comp.[[k]], ] - pred2[[dep]] <- ifelse(pred2[[dep]] == comp.[[k]][1], 0, 1) - check_cols <- c(dep, intersect(main_vars., names(pred2))) - if (length(check_cols) > 1) { - cors <- tryCatch( - stats::cor(pred2[, check_cols, drop = FALSE], use = "complete.obs"), - error = function(e) NULL - ) - if (is.null(cors) || any(is.na(cors))) { - y_ok <- x_ok <- FALSE - } else { - diag(cors) <- 0 - if (any(abs(cors) > 0.9999)) y_ok <- x_ok <- FALSE - } - } - } - } - sufficient_data <- y_ok && x_ok } @@ -198,113 +169,53 @@ QAPcssPermEst <- function(i, xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL - if (is.null(comp.)) { - # A fit inside the permutation loop runs `times` times, so a fitter's - # convergence warning would print once per draw and drown the console. - # The count of draws that failed outright is reported by - # `aggregate_perm_results()`, which is the number the user needs. - perm_fit <- tryCatch( - suppressWarnings(fit_qap_model(mod = mod., - pred = pred, - family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random., - reference = reference.)), - error = function(e) NULL - ) - if (is.null(perm_fit)) return(NULL) - - return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit., xi = xi_arg)) - } - - xresL <- vector("list", length(comp.)) - names(xresL) <- names(comp.) - - for (k in seq_along(comp.)) { - predK <- pred[pred[[dep]] %in% comp.[[k]], ] - predK[[dep]] <- ifelse(predK[[dep]] == comp.[[k]][1], 0, 1) - - # A fit inside the permutation loop runs `times` times, so a fitter's - # convergence warning would print once per draw and drown the console. - # The count of draws that failed outright is reported by - # `aggregate_perm_results()`, which is the number the user needs. - perm_fit <- tryCatch( - suppressWarnings(fit_qap_model(mod = mod., - pred = predK, - family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random., - reference = reference.)), - error = function(e) NULL - ) - if (is.null(perm_fit)) return(NULL) - - xresL[[k]] <- compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit.[[k]], xi = xi_arg) - } - - return(xresL) + # A fit inside the permutation loop runs `times` times, so a fitter's + # convergence warning would print once per draw and drown the console. + # The count of draws that failed outright is reported by + # `aggregate_perm_results()`, which is the number the user needs. + perm_fit <- tryCatch( + suppressWarnings(fit_qap_model(mod = mod., + pred = pred, + family = family., + use_fixest = use_fixest., + fixest_se_cluster = fixest_se_cluster., + use_robust_errors = use_robust_errors., + main_vars = main_vars., + has_random = has_random.)), + error = function(e) NULL + ) + if (is.null(perm_fit)) return(NULL) + + return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, + fit., xi = xi_arg)) } # Coefficient-table helper used by print.QAPCSS #' @keywords internal #' @noRd -glm_tab <- function(x, comp) { - if (!is.null(comp)) { - cat("\n\nComparison between", - x$comp[[comp]][1], "and", x$comp[[comp]][2]) - cat("\n\nCoefficients:\n") - - nc <- length(x$base[[comp]]$coefficients) - cmat <- matrix(NA, nrow = nc, ncol = 4) - cmat[, 1] <- format(round(as.numeric(x$base[[comp]]$coefficients), 3)) - cmat[, 2] <- format(x$lower[[comp]][2, ]) - cmat[, 3] <- format(x$larger[[comp]][2, ]) - cmat[, 4] <- format(x$abs[[comp]][2, ]) - if (x$permute == "predictor") cmat[1, 2:4] <- "*" - colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") - rownames(cmat) <- names(x$base[[comp]]$coefficients) - print.table(cmat) - - if (x$permute == "predictor") - cat("\n* The intercept has no significance test when predictors are permuted.\n") - - if (!is.null(x$base[[comp]]$base_model)) { - cat("\nAIC of base model:", format(stats::AIC(x$base[[comp]]$base_model))) - cat("\nBIC of base model:", format(stats::BIC(x$base[[comp]]$base_model))) - } - cat("\n") - } else { - cat("\n\nCoefficients:\n") - - nc <- length(x$base$coefficients) - cmat <- matrix(NA, nrow = nc, ncol = 4) - cmat[, 1] <- format(round(as.numeric(x$base$coefficients), 3)) - cmat[, 2] <- format(x$lower[2, ]) - cmat[, 3] <- format(x$larger[2, ]) - cmat[, 4] <- format(x$abs[2, ]) - if (x$permute == "predictor") cmat[1, 2:4] <- "*" - colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") - rownames(cmat) <- names(x$base$coefficients) - print.table(cmat) - - if (x$permute == "predictor") - cat("\n* The intercept has no significance test when predictors are permuted.\n") - - if (!is.null(x$base$base_model)) { - cat("\nAIC of base model:", format(stats::AIC(x$base$base_model))) - cat("\nBIC of base model:", format(stats::BIC(x$base$base_model))) - } - cat("\n") +glm_tab <- function(x) { + cat("\n\nCoefficients:\n") + + nc <- length(x$base$coefficients) + cmat <- matrix(NA, nrow = nc, ncol = 4) + cmat[, 1] <- format(round(as.numeric(x$base$coefficients), 3)) + cmat[, 2] <- format(x$lower[2, ]) + cmat[, 3] <- format(x$larger[2, ]) + cmat[, 4] <- format(x$abs[2, ]) + if (x$permute == "predictor") cmat[1, 2:4] <- "*" + colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") + rownames(cmat) <- names(x$base$coefficients) + print.table(cmat) + + if (x$permute == "predictor") + cat("\n* The intercept has no significance test when predictors are permuted.\n") + + if (!is.null(x$base$base_model)) { + cat("\nAIC of base model:", format(stats::AIC(x$base$base_model))) + cat("\nBIC of base model:", format(stats::BIC(x$base$base_model))) } + cat("\n") } @@ -322,8 +233,6 @@ QAPcss <- function(formula, family = "gaussian", groups = NULL, fixest_se_cluster = NULL, - reference = NULL, - comparison = NULL, use_robust_errors = FALSE, random_intercept_nets = FALSE, random_intercept_sender = FALSE, @@ -373,19 +282,6 @@ QAPcss <- function(formula, } mod <- stats::as.formula(mod_str) - if (has_random && family == "multinom") { - manynet::snet_warn( - c("Random intercepts are not implemented for the multinomial family.", - i = "Using {.fn nnet::multinom} instead.")) - has_random <- FALSE - } - if (!is.null(reference) && !is.character(reference) && family == "multinom") - reference <- as.character(reference) - if (use_robust_errors && family == "multinom") { - manynet::snet_warn( - "Robust standard errors are not implemented for the multinomial family.") - use_robust_errors <- FALSE - } if ((permute == "predictor") && (nx == 1)) permute <- "outcome" if (!directed && (ris || rir)) { manynet::snet_warn( @@ -443,33 +339,14 @@ QAPcss <- function(formula, fit <- list() - if (is.null(comparison)) { - fit$base <- fit_qap_model(mod = mod, - pred = pred, - family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, - use_robust_errors = use_robust_errors, - main_vars = main, - has_random = has_random, - reference = reference) - } else { - fit$base <- vector("list", length(comparison)) - names(fit$base) <- names(comparison) - for (k in seq_along(comparison)) { - predK <- pred[pred[[dep]] %in% comparison[[k]], ] - predK[[dep]] <- ifelse(predK[[dep]] == comparison[[k]][1], 0, 1) - fit$base[[k]] <- fit_qap_model(mod = mod, - pred = predK, - family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, - use_robust_errors = use_robust_errors, - main_vars = main, - has_random = has_random, - reference = reference) - } - } + fit$base <- fit_qap_model(mod = mod, + pred = pred, + family = family, + use_fixest = use_fixest, + fixest_se_cluster = fixest_se_cluster, + use_robust_errors = use_robust_errors, + main_vars = main, + has_random = has_random) old_plan <- setup_future_plan(strategy, ncores) on.exit({ @@ -486,7 +363,7 @@ QAPcss <- function(formula, diag. = diag, mod. = mod, groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, + fit. = fit$base, family. = family, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, @@ -494,65 +371,20 @@ QAPcss <- function(formula, has_random. = has_random, main_vars. = main, data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference + parsed. = parsed ) - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower <- agg$lower - fit$larger <- agg$larger - fit$abs <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } + agg <- aggregate_perm_results(res, times) + fit$lower <- agg$lower + fit$larger <- agg$larger + fit$abs <- agg$abs } else if (permute == "predictor") { - if (is.null(comparison)) { - if (family != "multinom") { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - } else { - ncat <- if (large) { - length(stats::na.omit(unique(as.vector(unlist(matlist[[dep]]))))) - } else { - length(stats::na.omit(unique(as.vector(matlist[[dep]])))) - } - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2 * (ncat - 1), ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - } - } else { - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) - for (k in seq_along(comparison)) { - n_coefs <- length(fit$base[[k]]$coefficients) - fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] - colnames(fit$lower[[k]]) <- colnames(fit$larger[[k]]) <- - colnames(fit$abs[[k]]) <- names(fit$base[[k]]$coefficients) - } - } + n_coefs <- length(fit$base$coefficients) + fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) + fit$larger <- fit$abs <- fit$lower + colnames(fit$lower) <- colnames(fit$larger) <- + colnames(fit$abs) <- names(fit$base$coefficients) for (xi in main) { test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] @@ -579,7 +411,7 @@ QAPcss <- function(formula, diag. = diag, mod. = mod, groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, + fit. = fit$base, family. = family, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, @@ -587,31 +419,17 @@ QAPcss <- function(formula, has_random. = has_random, main_vars. = main, data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference + parsed. = parsed ) - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower[, xi] <- agg$lower - fit$larger[, xi] <- agg$larger - fit$abs[, xi] <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } + agg <- aggregate_perm_results(res, times) + fit$lower[, xi] <- agg$lower + fit$larger[, xi] <- agg$larger + fit$abs[, xi] <- agg$abs } } - if (family == "binomial" && is.null(comparison)) { + if (family == "binomial") { fit$confusion_matrix <- probabilistic_confusion_matrix( actual = pred[[dep]], predicted_prob = stats::fitted(fit$base$base_model), @@ -625,22 +443,17 @@ QAPcss <- function(formula, fit$diag <- diag fit$directed <- directed fit$times <- times - fit$reference <- reference - fit$comp <- comparison fit$random <- c(sender = ris, receiver = rir, perceiver = rip, nets = rin) fit$robust_se <- use_robust_errors - if (is.null(comparison) && !is.null(fit$base$theta)) + if (!is.null(fit$base$theta)) fit$theta <- fit$base$theta - if (is.null(comparison) && !is.null(fit$base$zi_coefficients)) + if (!is.null(fit$base$zi_coefficients)) fit$zi_coefficients <- fit$base$zi_coefficients - if (family == "multinom") - names(fit)[names(fit) == "t"] <- "z" - class(fit) <- "QAPCSS" return(fit) } @@ -650,15 +463,10 @@ QAPcss <- function(formula, #' @noRd print.QAPCSS <- function(x, ...) { - if (x$family != "multinom") { - if (!any(x$random)) { - cat("\nGeneralized Linear Network Model for CSS\n\n") - } else { - cat("\nGeneralized Linear Mixed Network Model for CSS fit by REML\n\n") - } + if (!any(x$random)) { + cat("\nGeneralized Linear Network Model for CSS\n\n") } else { - cat("\nMultinomial Choice Network Model for CSS\n\n") - cat("The reference group was", format(paste0(x$reference, ".")), "\n") + cat("\nGeneralized Linear Mixed Network Model for CSS fit by REML\n\n") } cat("Estimator: Generalized Method-of-Moments.\n") @@ -694,40 +502,7 @@ print.QAPCSS <- function(x, ...) { cat("The outcome was treated as", format(paste0(.directed_label(x$directed), ".")), "\n") - if (x$family != "multinom") { - if (is.null(x$comp)) { - glm_tab(x, comp = x$comp) - } else { - for (mod in seq_along(x$comp)) { - glm_tab(x, comp = names(x$comp)[[mod]]) - } - } - } else { - cat("\nCoefficients:\n\n") - for (option in seq_len(nrow(x$base$coefficients))) { - cat(format(paste0("-- ", rownames(x$base$coefficients)[option], "\n"))) - - nc <- ncol(x$base$coefficients) - cmat <- matrix(NA, nrow = nc, ncol = 4) - row_idx <- option + nrow(x$base$coefficients) - cmat[, 1] <- format(as.numeric(x$base$coefficients[option, ])) - cmat[, 2] <- format(x$lower[row_idx, ]) - cmat[, 3] <- format(x$larger[row_idx, ]) - cmat[, 4] <- format(x$abs[row_idx, ]) - if (x$permute == "predictor") cmat[1, 2:4] <- "*" - colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") - rownames(cmat) <- colnames(x$base$coefficients) - print.table(cmat) - cat("\n\n") - } - - if (x$permute == "predictor") - cat("* The intercept has no significance test when predictors are permuted.\n") - - cat("\nAIC of base model:", format(stats::AIC(x$base$base_model))) - cat("\nBIC of base model:", format(stats::BIC(x$base$base_model))) - cat("\n") - } + glm_tab(x) if (!is.null(x$confusion_matrix)) { cat("\n") diff --git a/R/qap_engine.R b/R/qap_engine.R index 1e5c935..223e13a 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -19,8 +19,6 @@ QAPglm <- function(formula, strategy = "sequential", ncores = NULL, fixest_se_cluster = NULL, - comparison = NULL, - reference = NULL, random_intercept_nets = FALSE, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, @@ -84,10 +82,6 @@ QAPglm <- function(formula, names(pred)[names(pred) == "yv"] <- dep - if (!is.null(comparison) && is.null(reference)) { - reference <- NULL - } - fit <- list() rand_part <- "" @@ -95,33 +89,14 @@ QAPglm <- function(formula, if (ris) rand_part <- paste(rand_part, "+ (1|sv)") if (rir) rand_part <- paste(rand_part, "+ (1|rv)") - if (is.null(comparison)) { - fit$base <- fit_qap_model(mod = mod, - pred = pred, - family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, - use_robust_errors = use_robust_errors, - main_vars = main, - has_random = has_random, - reference = reference) - } else { - fit$base <- vector("list", length(comparison)) - names(fit$base) <- names(comparison) - for (k in seq_along(comparison)) { - predK <- pred[pred[[dep]] %in% comparison[[k]], ] - predK[[dep]] <- ifelse(predK[[dep]] == comparison[[k]][1], 0, 1) - fit$base[[k]] <- fit_qap_model(mod = mod, - pred = predK, - family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, - use_robust_errors = use_robust_errors, - main_vars = main, - has_random = has_random, - reference = reference) - } - } + fit$base <- fit_qap_model(mod = mod, + pred = pred, + family = family, + use_fixest = use_fixest, + fixest_se_cluster = fixest_se_cluster, + use_robust_errors = use_robust_errors, + main_vars = main, + has_random = has_random) # Double semi-partialling residualises a predictor against the others, so # with one predictor there are none and the scheme reduces to permuting the @@ -150,7 +125,7 @@ QAPglm <- function(formula, diag. = diag, mod. = mod, groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, + fit. = fit$base, family. = family, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, @@ -158,53 +133,20 @@ QAPglm <- function(formula, has_random. = has_random, main_vars. = main, data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference + parsed. = parsed ) - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower <- agg$lower - fit$larger <- agg$larger - fit$abs <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(comparison) - names(fit$larger) <- names(comparison) - names(fit$abs) <- names(comparison) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } + agg <- aggregate_perm_results(res, times) + fit$lower <- agg$lower + fit$larger <- agg$larger + fit$abs <- agg$abs } else if (permute == "predictor") { - if (is.null(comparison)) { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, - dimnames = list(c("perm_coefs", "perm_t"), - names(fit$base$coefficients))) - fit$larger <- fit$abs <- fit$lower - } else { - fit$lower <- fit$larger <- fit$abs <- - vector("list", length(comparison)) - names(fit$lower) <- names(fit$larger) <- - names(fit$abs) <- names(comparison) - for (k in seq_along(comparison)) { - n_coefs <- length(fit$base[[k]]$coefficients) - fit$lower[[k]] <- matrix(NA, nrow = 2, ncol = n_coefs, - dimnames = list(c("perm_coefs", "perm_t"), - names(fit$base[[k]]$coefficients))) - fit$larger[[k]] <- fit$abs[[k]] <- fit$lower[[k]] - } - } + n_coefs <- length(fit$base$coefficients) + fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs, + dimnames = list(c("perm_coefs", "perm_t"), + names(fit$base$coefficients))) + fit$larger <- fit$abs <- fit$lower for (xi in main) { xR <- residualise_predictor(xi, pred, main, @@ -222,7 +164,7 @@ QAPglm <- function(formula, diag. = diag, mod. = mod, groups. = groups, - fit. = if (is.null(comparison)) fit$base else fit$base, + fit. = fit$base, family. = family, use_fixest. = use_fixest, fixest_se_cluster. = fixest_se_cluster, @@ -230,51 +172,31 @@ QAPglm <- function(formula, has_random. = has_random, main_vars. = main, data_vars. = data_vars, - parsed. = parsed, - comp. = comparison, - reference. = reference + parsed. = parsed ) - if (is.null(comparison)) { - agg <- aggregate_perm_results(res, times) - fit$lower[, xi] <- agg$lower - fit$larger[, xi] <- agg$larger - fit$abs[, xi] <- agg$abs - } else { - res_valid <- Filter(Negate(is.null), res) - n_valid <- length(res_valid) - resL <- unlist(unlist(res_valid, recursive = FALSE), recursive = FALSE) - for (k in seq_along(comparison)) { - cn <- names(comparison)[k] - fit$lower[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".lower")], 0) / n_valid - fit$larger[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".larger")], 0) / n_valid - fit$abs[[k]][, xi] <- Reduce("+", resL[names(resL) == paste0(cn, ".abs")], 0) / n_valid - } - } + agg <- aggregate_perm_results(res, times) + fit$lower[, xi] <- agg$lower + fit$larger[, xi] <- agg$larger + fit$abs[, xi] <- agg$abs } } - if (is.null(comparison)) { - fit$coefficients <- fit$base$coefficients - fit$t <- fit$base$t - if (!is.null(fit$base$r.squared)) { - fit$r.squared <- fit$base$r.squared - fit$adj.r.squared <- fit$base$adj.r.squared - } - if (!is.null(fit$base$random.intercepts)) - fit$random.intercepts <- fit$base$random.intercepts - if (!is.null(fit$base$theta)) - fit$theta <- fit$base$theta - if (!is.null(fit$base$zi_coefficients)) - fit$zi_coefficients <- fit$base$zi_coefficients - if (!less_mem) fit$simple_fit <- fit$base$base_model - } else { - if (!less_mem) { - fit$simple_fits <- lapply(fit$base, `[[`, "base_model") - } + fit$coefficients <- fit$base$coefficients + fit$t <- fit$base$t + if (!is.null(fit$base$r.squared)) { + fit$r.squared <- fit$base$r.squared + fit$adj.r.squared <- fit$base$adj.r.squared } - - if (family == "binomial" && is.null(comparison)) { + if (!is.null(fit$base$random.intercepts)) + fit$random.intercepts <- fit$base$random.intercepts + if (!is.null(fit$base$theta)) + fit$theta <- fit$base$theta + if (!is.null(fit$base$zi_coefficients)) + fit$zi_coefficients <- fit$base$zi_coefficients + if (!less_mem) fit$simple_fit <- fit$base$base_model + + if (family == "binomial") { fit$confusion_matrix <- probabilistic_confusion_matrix( actual = pred[[dep]], predicted_prob = stats::fitted(fit$base$base_model), @@ -289,12 +211,10 @@ QAPglm <- function(formula, fit$times <- times fit$groups <- unique(unlist(groups)) fit$robust_se <- use_robust_errors - fit$comp <- comparison - fit$reference <- reference fit$pred <- pred fit$dep <- dep - if (family == "gaussian" && is.null(comparison)) { + if (family == "gaussian") { class(fit) <- "QAPRegression" } else { class(fit) <- "QAPGLM" @@ -320,9 +240,7 @@ QAPglmPermEst <- function(i, has_random., main_vars., data_vars., - parsed., - comp., - reference.) { + parsed.) { dep <- parsed.$dependent large <- is.list(matlist.[[dep]]) @@ -373,57 +291,23 @@ QAPglmPermEst <- function(i, xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL - if (is.null(comp.)) { - # A fit inside the permutation loop runs `times` times, so a fitter's - # convergence warning would print once per draw and drown the console. - # The count of draws that failed outright is reported by - # `aggregate_perm_results()`, which is the number the user needs. - perm_fit <- tryCatch( - suppressWarnings(fit_qap_model(mod = mod., - pred = pred, - family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random., - reference = reference.)), - error = function(e) NULL - ) - if (is.null(perm_fit)) return(NULL) - - return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit., xi = xi_arg)) - } - - xresL <- vector("list", length(comp.)) - names(xresL) <- names(comp.) - - for (k in seq_along(comp.)) { - predK <- pred[pred[[dep]] %in% comp.[[k]], ] - predK[[dep]] <- ifelse(predK[[dep]] == comp.[[k]][1], 0, 1) - - # A fit inside the permutation loop runs `times` times, so a fitter's - # convergence warning would print once per draw and drown the console. - # The count of draws that failed outright is reported by - # `aggregate_perm_results()`, which is the number the user needs. - perm_fit <- tryCatch( - suppressWarnings(fit_qap_model(mod = mod., - pred = predK, - family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random., - reference = reference.)), - error = function(e) NULL - ) - if (is.null(perm_fit)) return(NULL) - - xresL[[k]] <- compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit.[[k]], xi = xi_arg) - } - - return(xresL) + # A fit inside the permutation loop runs `times` times, so a fitter's + # convergence warning would print once per draw and drown the console. + # The count of draws that failed outright is reported by + # `aggregate_perm_results()`, which is the number the user needs. + perm_fit <- tryCatch( + suppressWarnings(fit_qap_model(mod = mod., + pred = pred, + family = family., + use_fixest = use_fixest., + fixest_se_cluster = fixest_se_cluster., + use_robust_errors = use_robust_errors., + main_vars = main_vars., + has_random = has_random.)), + error = function(e) NULL + ) + if (is.null(perm_fit)) return(NULL) + + return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, + fit., xi = xi_arg)) } diff --git a/R/qap_utils.R b/R/qap_utils.R index c2ee58c..2cff0b2 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -323,27 +323,11 @@ fit_qap_model <- function(...) { fixest_se_cluster = NULL, use_robust_errors = FALSE, main_vars = NULL, - has_random = FALSE, - reference = NULL) { + has_random = FALSE) { fit <- list() dep_var <- all.vars(mod)[1] nx <- length(main_vars) - if (family == "multinom") { - pred[[dep_var]] <- as.factor(pred[[dep_var]]) - if (!is.null(reference)) { - pred[[dep_var]] <- stats::relevel(pred[[dep_var]], ref = reference) - } - thisRequires("nnet", "for multinomial models") - base_model <- nnet::multinom(mod, data = pred, trace = FALSE) - fit$coefficients <- stats::coefficients(base_model) - fit$t <- stats::coefficients(base_model) / - summary(base_model)$standard.errors - fit$base_model <- base_model - return(fit) - } - - if (family == "zip") { if (has_random) { # The mixed variant needs {glmmTMB}, which is on feature/glmmtmb-mixed. diff --git a/man/regression.Rd b/man/regression.Rd index 4387e9f..4bb606b 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -46,14 +46,13 @@ reduces to \code{"outcome"} and says so. \item \code{strategy}: future plan, e.g. \code{"sequential"} (default), \code{"multisession"}. \item \code{family}: \code{"auto"} (default; gaussian for weighted networks, binomial for binary), \code{"gaussian"}, \code{"binomial"}, \code{"poisson"}, \code{"negbin"}, -\code{"zip"}, or \code{"multinom"}. +or \code{"zip"}. \item \code{directed}: logical, whether a tie from i to j differs from one from j to i. Read from \code{.data} unless given, and reported when read. \item \code{diag}: logical, include loops (default auto-detected). \item \code{seed}, \code{groups}, \code{ncores}: passed through to the engine. \item \code{use_robust_errors}: HC3 standard errors. \item \code{fixest_se_cluster}: cluster variable for fixest. -\item \code{reference}, \code{comparison}: multinomial / pairwise-comparison options. \item \code{random_intercept_nets} / \verb{_sender} / \verb{_receiver}: lme4-style REs. \item \code{less_mem}: drop the baseline model object from the return. }} @@ -84,8 +83,8 @@ The formula front-end is the familiar object. Internally the response and predictors are packed into matrices and handed to a QAP engine (ported from MrQAP) that supports: \itemize{ -\item gaussian, binomial, poisson, negbin, zero-inflated Poisson, and -multinomial families; +\item gaussian, binomial, poisson, negative binomial, and zero-inflated +Poisson families; \item two permutation schemes: \code{"predictor"} (Dekker's double semi-partialling) and \code{"outcome"}; \item random intercepts (lme4) and fixed effects (fixest); From 96fed9e893037c08999235ba9d27a547ffd43ebe Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:34:16 +0200 Subject: [PATCH 06/12] Branched off the fixest fixed effects, so a bar means one thing Moved to feature/fixest-fixed-effects. Reinstate with `git revert` of this commit. A bar in the formula meant two things. `getRHSNames()` read `y ~ a + b | c` as an lme4 random slope and rewrote it as `(b | c)`, while `parse_qap_formula()` read a bar without parentheses as a fixest fixed effect. Because the front end always adds the parentheses, fixed effects through the formula were unreachable, and {fixest} was only ever entered through `fixest_se_cluster`. With {fixest} on a branch a bar means an {lme4} random-effect term, and nothing else. `parse_qap_formula()` drops from three branches to one. Removes the `fixest_se_cluster` control and the fixest branch in fit_qap_model(). Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 1 - R/model_regression.R | 5 +- R/qap_css.R | 20 +---- R/qap_engine.R | 20 +---- R/qap_utils.R | 122 ++++++--------------------- man/regression.Rd | 3 +- tests/testthat/test-qap_estimators.R | 54 ------------ 7 files changed, 31 insertions(+), 194 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index af7480d..0c93126 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -39,7 +39,6 @@ Imports: reformulas Suggests: lme4, - fixest, MASS, pscl, testthat (>= 3.0.0) diff --git a/R/model_regression.R b/R/model_regression.R index f88bca9..32f7042 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -14,7 +14,7 @@ #' Poisson families; #' - two permutation schemes: `"predictor"` (Dekker's double semi-partialling) #' and `"outcome"`; -#' - random intercepts (lme4) and fixed effects (fixest); +#' - random intercepts (lme4); #' - robust (HC3) standard errors; #' - lists of networks, in which graphs that are missing any predictor are #' dropped with a warning and the remaining networks are pooled. @@ -54,7 +54,6 @@ #' - `diag`: logical, include loops (default auto-detected). #' - `seed`, `groups`, `ncores`: passed through to the engine. #' - `use_robust_errors`: HC3 standard errors. -#' - `fixest_se_cluster`: cluster variable for fixest. #' - `random_intercept_nets` / `_sender` / `_receiver`: lme4-style REs. #' - `less_mem`: drop the baseline model object from the return. #' @return An object of class `net_regression` inheriting from either @@ -142,7 +141,6 @@ net_regression <- function(formula, groups = ctrl$groups, strategy = ctrl$strategy, ncores = ctrl$ncores, - fixest_se_cluster = ctrl$fixest_se_cluster, random_intercept_nets = ctrl$random_intercept_nets, random_intercept_sender = ctrl$random_intercept_sender, random_intercept_receiver = ctrl$random_intercept_receiver, @@ -220,7 +218,6 @@ net_regression <- function(formula, groups = NULL, ncores = NULL, use_robust_errors = FALSE, - fixest_se_cluster = NULL, random_intercept_nets = FALSE, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, diff --git a/R/qap_css.R b/R/qap_css.R index a51c85d..fbbcb2d 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -93,8 +93,6 @@ QAPcssPermEst <- function(i, groups., fit., family., - use_fixest., - fixest_se_cluster., use_robust_errors., has_random., main_vars., @@ -177,8 +175,6 @@ QAPcssPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random.)), @@ -232,7 +228,6 @@ QAPcss <- function(formula, ncores = NULL, family = "gaussian", groups = NULL, - fixest_se_cluster = NULL, use_robust_errors = FALSE, random_intercept_nets = FALSE, random_intercept_sender = FALSE, @@ -242,7 +237,7 @@ QAPcss <- function(formula, if (!is.null(seed)) set.seed(seed) - parsed <- parse_qap_formula(formula, fixest_se_cluster) + parsed <- parse_qap_formula(formula) dep <- parsed$dependent main <- parsed$main data_vars <- intersect(parsed$all_data_vars, names(matlist)) @@ -273,13 +268,6 @@ QAPcss <- function(formula, rir = rir, rip = rip) mod_str <- paste(deparse(mod, width.cutoff = 500), collapse = " ") has_random <- grepl("\\(", mod_str) || parsed$has_random - use_fixest <- parsed$use_fixest - if (has_random && use_fixest) { - manynet::snet_warn( - c("Cannot combine {.pkg fixest} fixed effects with {.pkg lme4} random effects.", - i = "Using the random effects only.")) - use_fixest <- FALSE - } mod <- stats::as.formula(mod_str) if ((permute == "predictor") && (nx == 1)) permute <- "outcome" @@ -342,8 +330,6 @@ QAPcss <- function(formula, fit$base <- fit_qap_model(mod = mod, pred = pred, family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, main_vars = main, has_random = has_random) @@ -365,8 +351,6 @@ QAPcss <- function(formula, groups. = groups, fit. = fit$base, family. = family, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, has_random. = has_random, main_vars. = main, @@ -413,8 +397,6 @@ QAPcss <- function(formula, groups. = groups, fit. = fit$base, family. = family, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, has_random. = has_random, main_vars. = main, diff --git a/R/qap_engine.R b/R/qap_engine.R index 223e13a..a98fb61 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -18,7 +18,6 @@ QAPglm <- function(formula, groups = NULL, strategy = "sequential", ncores = NULL, - fixest_se_cluster = NULL, random_intercept_nets = FALSE, random_intercept_sender = FALSE, random_intercept_receiver = FALSE, @@ -27,7 +26,7 @@ QAPglm <- function(formula, if (!is.null(seed)) set.seed(seed) - parsed <- parse_qap_formula(formula, fixest_se_cluster) + parsed <- parse_qap_formula(formula) dep <- parsed$dependent main <- parsed$main data_vars <- intersect(parsed$all_data_vars, names(matlist)) @@ -43,13 +42,6 @@ QAPglm <- function(formula, mod <- build_internal_formula(formula, rin = rin, ris = ris, rir = rir) mod_str <- paste(deparse(mod, width.cutoff = 500), collapse = " ") has_random <- grepl("\\(", mod_str) || parsed$has_random - use_fixest <- parsed$use_fixest - if (has_random && use_fixest) { - manynet::snet_warn( - c("Cannot combine {.pkg fixest} fixed effects with {.pkg lme4} random effects.", - i = "Using the random effects only.")) - use_fixest <- FALSE - } mod <- stats::as.formula(mod_str) @@ -92,8 +84,6 @@ QAPglm <- function(formula, fit$base <- fit_qap_model(mod = mod, pred = pred, family = family, - use_fixest = use_fixest, - fixest_se_cluster = fixest_se_cluster, use_robust_errors = use_robust_errors, main_vars = main, has_random = has_random) @@ -127,8 +117,6 @@ QAPglm <- function(formula, groups. = groups, fit. = fit$base, family. = family, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, has_random. = has_random, main_vars. = main, @@ -166,8 +154,6 @@ QAPglm <- function(formula, groups. = groups, fit. = fit$base, family. = family, - use_fixest. = use_fixest, - fixest_se_cluster. = fixest_se_cluster, use_robust_errors. = use_robust_errors, has_random. = has_random, main_vars. = main, @@ -234,8 +220,6 @@ QAPglmPermEst <- function(i, groups., fit., family., - use_fixest., - fixest_se_cluster., use_robust_errors., has_random., main_vars., @@ -299,8 +283,6 @@ QAPglmPermEst <- function(i, suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., - use_fixest = use_fixest., - fixest_se_cluster = fixest_se_cluster., use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random.)), diff --git a/R/qap_utils.R b/R/qap_utils.R index 2cff0b2..b28aba7 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -8,46 +8,25 @@ #' Parse a QAP formula into its components #' @keywords internal #' @noRd -parse_qap_formula <- function(formula, fixest_se_cluster = NULL) { +parse_qap_formula <- function(formula) { dependent <- all.vars(formula)[1] + # A bar in the formula means one thing: an {lme4} random-effect term. The + # front end always writes it inside parentheses, and {fixest} -- which read a + # bare bar as a fixed effect -- is on feature/fixest-fixed-effects. formula_str <- paste(deparse(formula, width.cutoff = 500), collapse = " ") - has_pipe <- grepl("\\|", formula_str) - has_paren <- grepl("\\(", formula_str) - - if (has_pipe && has_paren) { - main <- all.vars(reformulas::nobars(formula))[-1] - fixed_effects <- NULL - use_fixest <- FALSE - has_random <- TRUE - all_data_vars <- main - } else if (has_pipe && !has_paren) { - main <- all.vars(formula[[3]][[2]]) - fixed_effects <- all.vars(formula[[3]][[3]]) - has_random <- FALSE - use_fixest <- TRUE - all_data_vars <- c(main, fixed_effects) - } else { - main <- all.vars(formula[-1]) - fixed_effects <- NULL - has_random <- FALSE - use_fixest <- !is.null(fixest_se_cluster) - all_data_vars <- main - } + has_random <- grepl("\\|", formula_str) - if (!is.null(fixest_se_cluster)) { - use_fixest <- TRUE - if (!(fixest_se_cluster %in% all_data_vars)) { - all_data_vars <- c(all_data_vars, fixest_se_cluster) - } + main <- if (has_random) { + all.vars(reformulas::nobars(formula))[-1] + } else { + all.vars(formula[-1]) } list(dependent = dependent, main = main, - fixed_effects = fixed_effects, has_random = has_random, - use_fixest = use_fixest, - all_data_vars = all_data_vars) + all_data_vars = main) } @@ -319,8 +298,6 @@ fit_qap_model <- function(...) { #' @keywords internal #' @noRd .fit_qap_model <- function(mod, pred, family, - use_fixest = FALSE, - fixest_se_cluster = NULL, use_robust_errors = FALSE, main_vars = NULL, has_random = FALSE) { @@ -352,70 +329,25 @@ fit_qap_model <- function(...) { } if (!has_random) { - if (use_fixest) { - thisRequires("fixest", "for fixed effects and clustered standard errors") - fe_family <- if (family == "negbin") "negbin" else family - base_model <- fixest::feglm(mod, data = pred, - family = fe_family, - cluster = fixest_se_cluster) - # {fixest} reports an intercept where no fixed effect is absorbed, and - # none where one is. Add the placeholder only in the second case; - # otherwise the coefficient vector carries two intercepts. - fe_coefs <- base_model$coefficients - fit$coefficients <- if ("(Intercept)" %in% names(fe_coefs)) { - fe_coefs - } else { - c("(Intercept)" = NA, fe_coefs) - } - resid <- stats::residuals(base_model) - - # `HC3()` and `vcov()` both return one standard error per estimated - # coefficient, so the placeholder is needed only where the intercept was - # absorbed and `fit$coefficients` carries an NA for it. - absorbed <- !("(Intercept)" %in% names(fe_coefs)) - if (use_robust_errors) { - xv <- as.matrix(pred[, main_vars, drop = FALSE]) - hc <- HC3(xv, resid) - fit$t <- if (absorbed) { - fit$coefficients / c(NA, hc[-1]) - } else { - fit$coefficients / hc - } - } else { - fe_se <- sqrt(diag(stats::vcov(base_model))) - fe_t <- fe_coefs / fe_se - fit$t <- if (absorbed) c("(Intercept)" = NA, fe_t) else fe_t - } - names(fit$t) <- names(fit$coefficients) - - if (family == "gaussian") { - r2s <- tryCatch(fixest::r2(base_model), error = function(e) NULL) - if (!is.null(r2s)) { - fit$r.squared <- r2s[["r2"]] - fit$adj.r.squared <- r2s[["ar2"]] - } - } + if (family == "gaussian") { + base_model <- stats::lm(mod, data = pred) + fit$r.squared <- summary(base_model)$r.squared + fit$adj.r.squared <- summary(base_model)$adj.r.squared + } else if (family == "negbin") { + thisRequires("MASS", "for negative binomial models") + base_model <- MASS::glm.nb(mod, data = pred) + fit$theta <- base_model$theta } else { - if (family == "gaussian") { - base_model <- stats::lm(mod, data = pred) - fit$r.squared <- summary(base_model)$r.squared - fit$adj.r.squared <- summary(base_model)$adj.r.squared - } else if (family == "negbin") { - thisRequires("MASS", "for negative binomial models") - base_model <- MASS::glm.nb(mod, data = pred) - fit$theta <- base_model$theta - } else { - base_model <- stats::glm(mod, data = pred, family = family) - } - fit$coefficients <- base_model$coefficients - resid <- stats::residuals(base_model) + base_model <- stats::glm(mod, data = pred, family = family) + } + fit$coefficients <- base_model$coefficients + resid <- stats::residuals(base_model) - if (use_robust_errors) { - xv <- as.matrix(pred[, main_vars, drop = FALSE]) - fit$t <- fit$coefficients / HC3(xv, resid) - } else { - fit$t <- summary(base_model)$coefficients[, 3] - } + if (use_robust_errors) { + xv <- as.matrix(pred[, main_vars, drop = FALSE]) + fit$t <- fit$coefficients / HC3(xv, resid) + } else { + fit$t <- summary(base_model)$coefficients[, 3] } } else { if (family == "gaussian") { diff --git a/man/regression.Rd b/man/regression.Rd index 4bb606b..71aa9d1 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -52,7 +52,6 @@ to i. Read from \code{.data} unless given, and reported when read. \item \code{diag}: logical, include loops (default auto-detected). \item \code{seed}, \code{groups}, \code{ncores}: passed through to the engine. \item \code{use_robust_errors}: HC3 standard errors. -\item \code{fixest_se_cluster}: cluster variable for fixest. \item \code{random_intercept_nets} / \verb{_sender} / \verb{_receiver}: lme4-style REs. \item \code{less_mem}: drop the baseline model object from the return. }} @@ -87,7 +86,7 @@ and handed to a QAP engine (ported from MrQAP) that supports: Poisson families; \item two permutation schemes: \code{"predictor"} (Dekker's double semi-partialling) and \code{"outcome"}; -\item random intercepts (lme4) and fixed effects (fixest); +\item random intercepts (lme4); \item robust (HC3) standard errors; \item lists of networks, in which graphs that are missing any predictor are dropped with a warning and the remaining networks are pooled. diff --git a/tests/testthat/test-qap_estimators.R b/tests/testthat/test-qap_estimators.R index b1a01b5..a647881 100644 --- a/tests/testthat/test-qap_estimators.R +++ b/tests/testthat/test-qap_estimators.R @@ -156,57 +156,3 @@ test_that("binomial and poisson random intercepts run", { expect_qap_shape(pois, COEFS3) }) - -# ---- fixed effects and clustered errors ------------------------------------ - -test_that("fixest reports one intercept, not two", { - skip_if_not_installed("fixest") - # `feglm()` reports an intercept where no fixed effect is absorbed. The engine - # used to prepend a placeholder regardless, giving two. - fit <- net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1, fixest_se_cluster = "sv")) - expect_qap_shape(fit, COEFS3) - expect_equal(sum(names(fit$coefficients) == "(Intercept)"), 1L) - expect_false(anyNA(fit$coefficients)) - expect_equal(length(fit$t), length(fit$coefficients)) -}) - -test_that("fixest coefficients match a direct feglm() fit", { - skip_if_not_installed("fixest") - g <- qap_net_gaussian() - ref <- qap_reference_data(FORM, g) - fe <- fixest::feglm(ref$formula, data = ref$pred, - family = "gaussian", cluster = "sv") - - fit <- net_regression(FORM, g, times = 10, - control = list(seed = 1, fixest_se_cluster = "sv")) - expect_equal(unname(fit$coefficients), unname(fe$coefficients)) -}) - -test_that("fixed effects and random effects together fall back to random", { - skip_if_not_installed("fixest") - skip_if_not_installed("lme4") - both <- suppressMessages(suppressWarnings( - net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1, fixest_se_cluster = "sv", - random_intercept_sender = TRUE)))) - random_only <- suppressMessages(suppressWarnings( - net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1, random_intercept_sender = TRUE)))) - expect_qap_shape(both, COEFS3) - # The fixed effects are dropped, so the fit is the random-effects one. - expect_equal(both$coefficients, random_only$coefficients) - expect_named(both$random.intercepts, "sv") -}) - -test_that("combining fixed and random effects warns", { - skip_if_not_installed("fixest") - skip_if_not_installed("lme4") - expect_snet_warning( - suppressMessages( - net_regression(FORM, qap_net_gaussian(), times = 10, - control = list(seed = 1, fixest_se_cluster = "sv", - random_intercept_sender = TRUE))), - "random effects") -}) - From f06fdb3c6a8d24cda5c27f5340156e4d61cb08ee Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:34:58 +0200 Subject: [PATCH 07/12] Recorded the branched extensions in NEWS Co-Authored-By: Claude Opus 5 --- NEWS.md | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/NEWS.md b/NEWS.md index b914f08..e98f996 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,11 @@ ## Package +- Branched off five model extensions, to settle the architecture first + - Each is on its own `feature/*` branch, and each strip is one commit that + `git revert` reinstates + - `Suggests` falls from eight modelling packages to three + - See the Github issues for the order they come back in - Renamed the engine's vocabulary to the front end's, so one word means one thing on both sides of the seam - `reps` is now `times`, everywhere including on the returned fit @@ -20,6 +25,22 @@ ## Regression +- Removed the `torch` GPU path (`feature/torch-gpu`) + - Gaussian only, duplicated for CSS, no test, and no hosted runner has a + CUDA device; `{torch}` in Suggests broke the CI build +- Removed the `gmm` estimator and the `estimator` control (`feature/gmm-estimator`) + - It warned that the coefficient covariance matrix was singular on every + family, on well-conditioned data +- Removed the mixed negbin and mixed zip paths (`feature/glmmtmb-mixed`) + - `{glmmTMB}` carries 62 recursive dependencies and must match `{TMB}` + - The standard `negbin` and `zip` paths are unaffected +- Removed `family = "multinom"` and the `comparison`/`reference` controls + (`feature/multinomial-comparison`) + - Unreachable from the front end, and its pairwise branch forked both + engines at 21 points +- Removed the `fixest_se_cluster` control (`feature/fixest-fixed-effects`) + - A bar in the formula now means an `{lme4}` random-effect term, and + nothing else; `parse_qap_formula()` drops from three branches to one - Fixed `net_regression()` failing on a two-mode network with more columns than rows (closing #4) - The validity mask was built as rows-by-rows, so a wider predictor extended From b35813ce308cce3c7e131e825bd077cc8aba4394 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 15:39:13 +0200 Subject: [PATCH 08/12] Documented the parked extensions and their issues Co-Authored-By: Claude Opus 5 --- .github/CONTRIBUTING.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index 23320c1..b11a16a 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -132,6 +132,25 @@ Make it clear when you are referring to functions by adding backticks and parent e.g. `a_function()`, and arguments by adding an equals sign, e.g. `argument=`. Argument values or variables can be in double quotation marks, e.g. "value". +## Parked extensions + +Five model extensions sit on `feature/*` branches while the architecture +settles, each tracked by a Github issue and each reinstated by reverting one +commit on `develop`: + +| Branch | Removes | Issue | +|---|---|---| +| `feature/multinomial-comparison` | `family = "multinom"`, and the `comparison`/`reference` controls | [#7](https://github.com/stocnet/infernet/issues/7) | +| `feature/fixest-fixed-effects` | the `fixest_se_cluster` control and the `{fixest}` branch | [#8](https://github.com/stocnet/infernet/issues/8) | +| `feature/glmmtmb-mixed` | mixed negbin and mixed zip | [#9](https://github.com/stocnet/infernet/issues/9) | +| `feature/gmm-estimator` | the `estimator` control and `R/qap_gmm.R` | [#10](https://github.com/stocnet/infernet/issues/10) | +| `feature/torch-gpu` | `R/qap_gpu.R` and the `use_gpu` control | [#11](https://github.com/stocnet/infernet/issues/11) | + +Do not reinstate one by reverting onto `develop` without reading its issue: +several need rewriting against the merged engine rather than reverting onto it. +Do not add a new model family that needs a new `Suggests` package until the +two engines are one, for the same reason these left. + ## Package architecture ### Project overview From 6d82785e13c767bc691425475b40e0c0a6feebed Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sat, 5 Sep 2026 16:10:59 +0200 Subject: [PATCH 09/12] Merged the two engines into one, with a shape strategy `QAPglm()` and `QAPcss()` were two 200-line functions that were 55% the same code, so every fix had to be made twice, and one of them was made in only one place. They are now `QAPengine()`, which fits both a dyadic network and a cognitive social structure. What the two shapes do differently is four functions in R/qap_shapes.R: how to vectorise a network into rows, how to draw a permutation, how to put a residualised predictor back, and which random intercepts exist. A fifth field says how many permutations to redraw before giving up: one for a dyadic network, 10,000 for a sparse CSS array. A random-intercept slot a shape does not list now aborts by name, so a perceiver intercept on a dyadic network says so rather than building a formula that will not parse. The engine files fall from 791 lines to 552. Also registers `print.QAPCSS()`, which was never an S3 method, and removes a dangling GMM line left in it by an earlier strip. Co-Authored-By: Claude Opus 5 --- .github/CONTRIBUTING.md | 31 ++- NAMESPACE | 1 + NEWS.md | 13 ++ R/model_regression.R | 5 +- R/qap_css.R | 348 +-------------------------- R/qap_engine.R | 349 ++++++++++++++-------------- R/qap_shapes.R | 106 +++++++++ tests/testthat/test-qap_shape_css.R | 89 +++++++ 8 files changed, 416 insertions(+), 526 deletions(-) create mode 100644 R/qap_shapes.R create mode 100644 tests/testthat/test-qap_shape_css.R diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index b11a16a..028e4aa 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -225,9 +225,10 @@ the internal engine ported from `MrQAP` is named `qap_*.R`. |---|---| | `model_tests.R` | the test family: `test_random()` (CUG), `test_configuration()`, `test_permutation()` (QAP), and `print.network_test()` | | `model_regression.R` | `net_regression()`, the formula front end (`convertToMatrixList()`, `getRHSNames()`, `specificationAdvice()`), and the `print.*` methods for its results | -| `qap_engine.R` | `QAPglm()` and `QAPglmPermEst()` — the matrix-level engine that performs the baseline fit and the permutation inference | +| `qap_engine.R` | `QAPengine()` and `QAPPermEst()` — the one matrix-level engine, for both a dyadic network and a cognitive social structure | +| `qap_shapes.R` | the four things the two shapes do differently, and nothing else | | `qap_utils.R` | formula parsing, input validation, `future` plumbing, matrix permutation (`RMPerm()`), the model-fitting dispatcher `fit_qap_model()`, and the permutation aggregators | -| `qap_css.R` | `QAPcss()` and `QAPcssPermEst()` — the parallel engine for cognitive social structures | +| `qap_css.R` | what a CSS needs that a dyadic network does not: a vectoriser for a three-dimensional array, and a print method | | `qap_gmm.R` | GMM moment conditions and residual functions for the `estimator = "gmm"` path | | `qap_gpu.R` | the optional `{torch}` batch OLS path, `gpu_batch_ols()` | | `qap_confusion.R` | probabilistic confusion matrices for binary outcomes | @@ -271,7 +272,31 @@ falls back to `"outcome"` and says so. Permuted coefficients and test statistics are then compared against the baseline by `compare_perm_to_baseline()` and reduced to `lower`/`larger`/`abs` p-value matrices by `aggregate_perm_results()`. -`QAPcss()` mirrors this same permute-refit-aggregate architecture for CSS data. +### One engine, two shapes + +`QAPengine()` fits a dyadic network and a cognitive social structure through the +same skeleton. They differ in four places and nowhere else, and those four live +in a *shape* returned by `.qap_shape()` +([R/qap_shapes.R](../R/qap_shapes.R)): + +| Field | Dyadic | Cognitive | +|---|---|---| +| `vectorise()` | `make_qap_data()`, one row per dyad | `make_css_data()`, one row per dyad per perceiver | +| `permute()` | `RMPerm()` | `RMPerm(CSS = TRUE)` | +| `unresidualise()` | `residuals_to_matrix()` | `residuals_to_array()` | +| `rand_slots` | sender, receiver, network | and perceiver | + +A fifth field, `max_trials`, says how many permutations to redraw before giving +up: one for a dyadic network, since a degenerate draw is simply dropped and +counted, and 10,000 for a CSS, whose sparse arrays often permute into an +outcome with a single value. + +Add a shape rather than a second engine. A random-intercept slot a shape does +not list cannot be requested, so a perceiver intercept on a dyadic network +aborts by name rather than producing a formula that will not parse. + +Before this merge the two were `QAPglm()` and `QAPcss()`, 55% the same code, and +every fix had to be made twice. One of them was made in only one place. The formula front end accepts these terms, and a new one should be added to `getRHSNames()` and `convertToMatrixList()` together: diff --git a/NAMESPACE b/NAMESPACE index b0cf7af..a3a85ec 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(print,QAPCSS) S3method(print,net_regression) S3method(print,network_test) export(net_regression) diff --git a/NEWS.md b/NEWS.md index e98f996..ecb6a71 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,15 @@ ## Package +- Merged the two engines into one, `QAPengine()` + - `QAPglm()` and `QAPcss()` were 55% the same code, so every fix had to be + made twice; one of them was made in only one place + - What the two shapes do differently is now four functions in + `R/qap_shapes.R`: how to vectorise, how to permute, how to put a + residualised predictor back, and which random intercepts exist + - A random intercept a shape does not have now aborts by name, so a + perceiver intercept on a dyadic network says so + - The engine files fall from 791 lines to 552, with no duplication left - Branched off five model extensions, to settle the architecture first - Each is on its own `feature/*` branch, and each strip is one commit that `git revert` reinstates @@ -59,6 +68,10 @@ ## Tests +- Added `test-qap_shape_css.R`, which fits a cognitive social structure through + the merged engine + - `net_regression()` has no CSS entry point yet, so the merge would + otherwise be untested on the shape it was merged for - Added a wide two-mode fixture and two regression tests for #4 - Added `test-qap_reporting.R`, which runs with `snet_verbosity = "verbose"` - Informational output is silent in every other test, so a message that diff --git a/R/model_regression.R b/R/model_regression.R index 32f7042..59d70da 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -129,9 +129,10 @@ net_regression <- function(formula, ctrl$diag <- isTRUE(manynet::is_complex(first_graph)) } - fit <- QAPglm( + fit <- QAPengine( formula = formula, - matlist = matlist, + matlist = matlist, + css = FALSE, family = ctrl$family, directed = ctrl$directed, diag = ctrl$diag, diff --git a/R/qap_css.R b/R/qap_css.R index fbbcb2d..b7164dd 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -1,8 +1,8 @@ -# CSS (Cognitive Social Structure) engine ------------------------------------ +# Cognitive social structures -------------------------------------------------- # -# 3D-array QAP engine. Not currently exposed through net_regression() but -# ported here so the machinery is available once manynet finalises its CSS -# representation. All functions internal. +# Internal. What a CSS needs that a dyadic network does not: a vectoriser for a +# sender-by-receiver-by-perceiver array, and a print method. Everything else it +# shares with the dyadic case, and lives in R/qap_engine.R. #' @keywords internal #' @noRd @@ -18,7 +18,6 @@ array_to_vector <- function(ar, directed., diag.) { return(v) } - #' @keywords internal #' @noRd make_css_data <- function(y, x, nets, diag, directed) { @@ -81,112 +80,6 @@ make_css_data <- function(y, x, nets, diag, directed) { return(list(pred = pred, valid = valid)) } - -#' @keywords internal -#' @noRd -QAPcssPermEst <- function(i, - matlist., - perm_var., - directed., - diag., - mod., - groups., - fit., - family., - use_robust_errors., - has_random., - main_vars., - data_vars., - parsed.) { - - dep <- parsed.$dependent - large <- is.list(matlist.[[dep]]) - - y_cat <- stats::na.omit(unique(as.vector(unlist(matlist.[[dep]])))) - - sufficient_data <- FALSE - trial <- 0 - max_trials <- 10000 - - while (!sufficient_data && trial < max_trials) { - trial <- trial + 1 - - d <- matlist. - if (is.null(perm_var.)) { - if (!large) { - d[[dep]] <- RMPerm(d[[dep]], groups., CSS = TRUE) - } else { - d[[dep]] <- lapply(d[[dep]], RMPerm, groups = groups., CSS = TRUE) - } - } else { - if (!large) { - d[[perm_var.]] <- RMPerm(d[[perm_var.]], groups., CSS = TRUE) - } else { - d[[perm_var.]] <- lapply(d[[perm_var.]], RMPerm, - groups = groups., CSS = TRUE) - } - } - - if (!large) { - x_list <- lapply(data_vars., function(v) d[[v]]) - names(x_list) <- data_vars. - pred <- make_css_data(y = d[[dep]], x = x_list, - nets = 1, - diag = diag., directed = directed.)$pred - } else { - pred_list <- vector("list", length(d[[dep]])) - for (gr in seq_along(d[[dep]])) { - xgr <- lapply(data_vars., function(v) d[[v]][[gr]]) - names(xgr) <- data_vars. - pred_list[[gr]] <- make_css_data(y = d[[dep]][[gr]], x = xgr, - nets = gr, - diag = diag., directed = directed.)$pred - } - pred <- do.call(rbind, pred_list) - } - - names(pred)[names(pred) == "yv"] <- dep - - y_ok <- length(stats::na.omit(unique(pred[[dep]]))) > 1 - - x_ok <- TRUE - num_preds <- pred[, data_vars.[data_vars. %in% names(pred)], drop = FALSE] - num_preds <- num_preds[, sapply(num_preds, is.numeric), drop = FALSE] - if (ncol(num_preds) > 0) { - x_ok <- all(sapply(num_preds, function(col) length(unique(col)) > 1)) - } - - sufficient_data <- y_ok && x_ok - } - - if (trial >= max_trials) { - manynet::snet_abort( - c("Cannot find a valid permutation after {max_trials} trials.", - i = "The network may be too sparse, or too many cells may be missing.")) - } - - xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL - - # A fit inside the permutation loop runs `times` times, so a fitter's - # convergence warning would print once per draw and drown the console. - # The count of draws that failed outright is reported by - # `aggregate_perm_results()`, which is the number the user needs. - perm_fit <- tryCatch( - suppressWarnings(fit_qap_model(mod = mod., - pred = pred, - family = family., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random.)), - error = function(e) NULL - ) - if (is.null(perm_fit)) return(NULL) - - return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit., xi = xi_arg)) -} - - # Coefficient-table helper used by print.QAPCSS #' @keywords internal #' @noRd @@ -214,235 +107,9 @@ glm_tab <- function(x) { cat("\n") } - -#' @keywords internal -#' @noRd -QAPcss <- function(formula, - matlist, - directed = TRUE, - diag = FALSE, - permute = "outcome", - times = 1000, - seed = NULL, - strategy = "sequential", - ncores = NULL, - family = "gaussian", - groups = NULL, - use_robust_errors = FALSE, - random_intercept_nets = FALSE, - random_intercept_sender = FALSE, - random_intercept_receiver = FALSE, - random_intercept_perceiver = FALSE, - less_mem = FALSE) { - - if (!is.null(seed)) set.seed(seed) - - parsed <- parse_qap_formula(formula) - dep <- parsed$dependent - main <- parsed$main - data_vars <- intersect(parsed$all_data_vars, names(matlist)) - nx <- length(main) - - validate_qap_input(matlist, parsed, css = TRUE) - large <- is.list(matlist[[dep]]) - - if (!large) { - y <- matlist[[dep]] - if (length(dim(y)) != 3) - manynet::snet_abort( - "The dependent variable {.val {dep}} must be a 3-dimensional array of sender, receiver, and perceiver.") - } else { - for (i in seq_along(matlist[[dep]])) { - if (length(dim(matlist[[dep]][[i]])) != 3) - manynet::snet_abort( - "Network {i} of the dependent variable {.val {dep}} must be a 3-dimensional array.") - } - } - - rin <- random_intercept_nets - rip <- random_intercept_perceiver - ris <- random_intercept_sender - rir <- random_intercept_receiver - - mod <- build_internal_formula(formula, rin = rin, ris = ris, - rir = rir, rip = rip) - mod_str <- paste(deparse(mod, width.cutoff = 500), collapse = " ") - has_random <- grepl("\\(", mod_str) || parsed$has_random - mod <- stats::as.formula(mod_str) - - if ((permute == "predictor") && (nx == 1)) permute <- "outcome" - if (!directed && (ris || rir)) { - manynet::snet_warn( - c("An undirected network has no senders or receivers.", - i = "Setting the sender and receiver random intercepts to {.val FALSE}.")) - ris <- rir <- FALSE - } - if (diag) - manynet::snet_warn( - "Results may not be valid where the diagonal is included.") - - rand_part <- "" - if (rin) rand_part <- paste(rand_part, "+ (1|nv)") - if (rip) rand_part <- paste(rand_part, "+ (1|pv)") - if (ris) rand_part <- paste(rand_part, "+ (1|sv)") - if (rir) rand_part <- paste(rand_part, "+ (1|rv)") - - if (!large) { - n <- dim(matlist[[dep]])[1] - if (!is.null(groups)) { - if (length(groups) != n) - manynet::snet_abort( - "{.arg groups} is of length {length(groups)}, but the network has {n} nodes.") - groups <- as.factor(groups) - } else { - groups <- as.factor(rep(1, n)) - } - } - - valid <- NULL; valid_list <- NULL - if (!large) { - x_list <- lapply(data_vars, function(v) matlist[[v]]) - names(x_list) <- data_vars - cssd <- make_css_data(y = matlist[[dep]], x = x_list, - nets = 1, - diag = diag, directed = directed) - pred <- cssd$pred - valid <- cssd$valid - } else { - pred_list <- vector("list", length(matlist[[dep]])) - valid_list <- vector("list", length(matlist[[dep]])) - for (gr in seq_along(matlist[[dep]])) { - xgr <- lapply(data_vars, function(v) matlist[[v]][[gr]]) - names(xgr) <- data_vars - cssd <- make_css_data(y = matlist[[dep]][[gr]], x = xgr, - nets = gr, - diag = diag, directed = directed) - pred_list[[gr]] <- cssd$pred - valid_list[[gr]] <- cssd$valid - } - pred <- do.call(rbind, pred_list) - } - - names(pred)[names(pred) == "yv"] <- dep - - fit <- list() - - fit$base <- fit_qap_model(mod = mod, - pred = pred, - family = family, - use_robust_errors = use_robust_errors, - main_vars = main, - has_random = has_random) - - old_plan <- setup_future_plan(strategy, ncores) - on.exit({ - future::plan(old_plan) - options(future.globals.maxSize = attr(old_plan, "old_maxSize")) - }, add = TRUE) - - if (permute == "outcome") { - res <- run_permutations( - times, QAPcssPermEst, - matlist. = matlist, - perm_var. = NULL, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = fit$base, - family. = family, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed - ) - - agg <- aggregate_perm_results(res, times) - fit$lower <- agg$lower - fit$larger <- agg$larger - fit$abs <- agg$abs - - } else if (permute == "predictor") { - n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) - fit$larger <- fit$abs <- fit$lower - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) - - for (xi in main) { - test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] - if (!is.numeric(test_val)) { - manynet::snet_warn( - c("Cannot residualise the non-numeric predictor {.val {xi}}.", - i = "Skipping double semi-partialling for this predictor.")) - next - } - - xR <- residualise_predictor(xi, pred, main, - has_random = has_random, - rand_formula = rand_part) - - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_array(xR, matlist[[xi]], valid, pred, - large, valid_list) - - res <- run_permutations( - times, QAPcssPermEst, - matlist. = matlist_resid, - perm_var. = xi, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = fit$base, - family. = family, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed - ) - - agg <- aggregate_perm_results(res, times) - fit$lower[, xi] <- agg$lower - fit$larger[, xi] <- agg$larger - fit$abs[, xi] <- agg$abs - } - } - - if (family == "binomial") { - fit$confusion_matrix <- probabilistic_confusion_matrix( - actual = pred[[dep]], - predicted_prob = stats::fitted(fit$base$base_model), - n_draws = 1000, seed = seed - ) - } - - fit$permute <- permute - fit$family <- family - fit$groups <- unique(unlist(groups)) - fit$diag <- diag - fit$directed <- directed - fit$times <- times - fit$random <- c(sender = ris, - receiver = rir, - perceiver = rip, - nets = rin) - fit$robust_se <- use_robust_errors - - if (!is.null(fit$base$theta)) - fit$theta <- fit$base$theta - if (!is.null(fit$base$zi_coefficients)) - fit$zi_coefficients <- fit$base$zi_coefficients - - class(fit) <- "QAPCSS" - return(fit) -} - - -#' @keywords internal -#' @noRd +# Registered so that a CSS fit prints as a model rather than as a list. The +# class is internal, so the method is registered without a help topic. +#' @exportS3Method print QAPCSS print.QAPCSS <- function(x, ...) { if (!any(x$random)) { @@ -451,7 +118,6 @@ print.QAPCSS <- function(x, ...) { cat("\nGeneralized Linear Mixed Network Model for CSS fit by REML\n\n") } - cat("Estimator: Generalized Method-of-Moments.\n") if (!is.null(x$theta)) cat("Negative binomial dispersion (theta):", format(round(x$theta, 4)), "\n") if (!is.null(x$zi_coefficients)) { diff --git a/R/qap_engine.R b/R/qap_engine.R index a98fb61..c853c1f 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -1,86 +1,104 @@ -# QAPglm engine -------------------------------------------------------------- +# QAP engine ------------------------------------------------------------------ # -# Internal. The matrix-level engine that drives net_regression(): performs -# the baseline fit and the QAP / QAP-DSP permutation inference on a pre-built -# list of matrices. Ported from MrQAP::QAPglm(). All parallelism uses the -# `future` framework via `run_permutations()`; no progressr integration. +# Internal. The matrix-level engine behind net_regression(): it performs the +# baseline fit and the permutation inference on a pre-built list of matrices, +# for both a dyadic network and a cognitive social structure. +# +# The two used to be `QAPglm()` and `QAPcss()`, two 200-line functions that were +# 55% the same code. Everything they shared is here; everything they did not is +# in R/qap_shapes.R. Ported from MrQAP. All parallelism uses the `future` +# framework via `run_permutations()`. #' @keywords internal #' @noRd -QAPglm <- function(formula, - matlist, - family = "gaussian", - directed = TRUE, - diag = FALSE, - permute = "predictor", - times = 1000, - seed = NULL, - groups = NULL, - strategy = "sequential", - ncores = NULL, - random_intercept_nets = FALSE, - random_intercept_sender = FALSE, - random_intercept_receiver = FALSE, - use_robust_errors = FALSE, - less_mem = FALSE) { +QAPengine <- function(formula, + matlist, + css = FALSE, + family = "gaussian", + directed = TRUE, + diag = FALSE, + permute = "predictor", + times = 1000, + seed = NULL, + groups = NULL, + strategy = "sequential", + ncores = NULL, + random_intercept_nets = FALSE, + random_intercept_sender = FALSE, + random_intercept_receiver = FALSE, + random_intercept_perceiver = FALSE, + use_robust_errors = FALSE, + less_mem = FALSE) { if (!is.null(seed)) set.seed(seed) + shape <- .qap_shape(css) parsed <- parse_qap_formula(formula) - dep <- parsed$dependent - main <- parsed$main - data_vars <- intersect(parsed$all_data_vars, names(matlist)) + dep <- parsed$dependent + main <- parsed$main + data_vars <- intersect(parsed$all_data_vars, names(matlist)) - validate_qap_input(matlist, parsed, css = FALSE) + validate_qap_input(matlist, parsed, css = css) large <- is.list(matlist[[dep]]) + # ---- random intercepts ---------------------------------------------------- + + requested <- c(nets = random_intercept_nets, + sender = random_intercept_sender, + receiver = random_intercept_receiver, + perceiver = random_intercept_perceiver) + # A slot a shape does not have cannot be asked for. A dyadic network has no + # perceiver, so a perceiver intercept there is a mistake worth naming. + unavailable <- names(requested)[requested & + !(names(requested) %in% names(shape$rand_slots))] + if (length(unavailable) > 0) { + label <- shape$label + manynet::snet_abort( + "A {label} has no {.val {unavailable}} random intercept{?s}.") + } + if (!directed && (requested[["sender"]] || requested[["receiver"]])) { + manynet::snet_warn( + c("An undirected network has no senders or receivers.", + i = "Setting the sender and receiver random intercepts to {.val FALSE}.")) + requested[c("sender", "receiver")] <- FALSE + } - rin <- random_intercept_nets - ris <- random_intercept_sender - rir <- random_intercept_receiver + slots <- shape$rand_slots + active <- slots[requested[names(slots)]] + # `paste0()` folds a zero-length argument into "", so an empty set of + # intercepts would otherwise produce the unparseable term " + (1|)". + rand_part <- if (length(active) == 0) { + "" + } else { + paste0(" + (1|", active, ")", collapse = "") + } - mod <- build_internal_formula(formula, rin = rin, ris = ris, rir = rir) - mod_str <- paste(deparse(mod, width.cutoff = 500), collapse = " ") - has_random <- grepl("\\(", mod_str) || parsed$has_random + mod <- stats::as.formula(paste( + paste(deparse(formula, width.cutoff = 500), collapse = " "), rand_part)) + has_random <- length(active) > 0 || parsed$has_random - mod <- stats::as.formula(mod_str) + if (diag && css) + manynet::snet_warn( + "Results may not be valid where the diagonal is included.") - if (!large) { - pred <- make_qap_data(y = matlist[[dep]], - x = matlist[data_vars], - g = groups, - diag = diag, - directed = directed, - net = 1, - perm = FALSE, - xi = NULL) - } else { - pred_list <- vector("list", length(matlist[[dep]])) - for (net in seq_along(matlist[[dep]])) { - x2 <- lapply(data_vars, function(v) matlist[[v]][[net]]) - names(x2) <- data_vars - g2 <- if (!is.null(groups)) groups[[net]] else NULL - pred_list[[net]] <- make_qap_data(y = matlist[[dep]][[net]], - x = x2, - g = g2, - diag = diag, - directed = directed, - net = net, - perm = FALSE, - xi = NULL) - } - pred <- do.call(rbind, pred_list) + # ---- vectorise ------------------------------------------------------------ + + if (!large && !is.null(groups)) { + n <- dim(matlist[[dep]])[1] + if (length(groups) != n) + manynet::snet_abort( + "{.arg groups} is of length {length(groups)}, but the network has {n} nodes.") + groups <- as.factor(groups) } + vec <- .vectorise_matlist(shape, matlist, dep, data_vars, groups, + diag, directed, large) + pred <- vec$pred names(pred)[names(pred) == "yv"] <- dep - fit <- list() - - rand_part <- "" - if (rin) rand_part <- paste(rand_part, "+ (1|nv)") - if (ris) rand_part <- paste(rand_part, "+ (1|sv)") - if (rir) rand_part <- paste(rand_part, "+ (1|rv)") + # ---- baseline ------------------------------------------------------------- + fit <- list() fit$base <- fit_qap_model(mod = mod, pred = pred, family = family, @@ -93,37 +111,29 @@ QAPglm <- function(formula, # outcome. Say so: the result would otherwise report a scheme nobody chose. if ((permute == "predictor") && (length(main) == 1)) { permute <- "outcome" - # `snet_info()` pastes its arguments, so give it separate strings rather - # than a named vector: a named vector loses its bullets and runs together. manynet::snet_info( "Permuting {.val outcome}, not {.val predictor}:", "with one predictor there is nothing to residualise it against.") } + # ---- permute -------------------------------------------------------------- + old_plan <- setup_future_plan(strategy, ncores) on.exit({ future::plan(old_plan) options(future.globals.maxSize = attr(old_plan, "old_maxSize")) }, add = TRUE) - if (permute == "outcome") { - res <- run_permutations( - times, QAPglmPermEst, - matlist. = matlist, - perm_var. = NULL, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = fit$base, - family. = family, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed - ) + perm_args <- list(shape. = shape, directed. = directed, diag. = diag, + mod. = mod, groups. = groups, fit. = fit$base, + family. = family, use_robust_errors. = use_robust_errors, + has_random. = has_random, main_vars. = main, + data_vars. = data_vars, parsed. = parsed) + if (permute == "outcome") { + res <- do.call(run_permutations, + c(list(times, QAPPermEst, matlist. = matlist, + perm_var. = NULL), perm_args)) agg <- aggregate_perm_results(res, times) fit$lower <- agg$lower fit$larger <- agg$larger @@ -137,30 +147,24 @@ QAPglm <- function(formula, fit$larger <- fit$abs <- fit$lower for (xi in main) { + test_val <- if (!large) matlist[[xi]] else matlist[[xi]][[1]] + if (!is.numeric(test_val)) { + manynet::snet_warn( + c("Cannot residualise the non-numeric predictor {.val {xi}}.", + i = "Skipping double semi-partialling for this predictor.")) + next + } + xR <- residualise_predictor(xi, pred, main, has_random = has_random, rand_formula = rand_part) - matlist_resid <- matlist - matlist_resid[[xi]] <- residuals_to_matrix(xR, matlist[[xi]], pred, large) - - res <- run_permutations( - times, QAPglmPermEst, - matlist. = matlist_resid, - perm_var. = xi, - directed. = directed, - diag. = diag, - mod. = mod, - groups. = groups, - fit. = fit$base, - family. = family, - use_robust_errors. = use_robust_errors, - has_random. = has_random, - main_vars. = main, - data_vars. = data_vars, - parsed. = parsed - ) + matlist_resid[[xi]] <- shape$unresidualise(xR, matlist[[xi]], pred, large, + vec$valid, vec$valid_list) + res <- do.call(run_permutations, + c(list(times, QAPPermEst, matlist. = matlist_resid, + perm_var. = xi), perm_args)) agg <- aggregate_perm_results(res, times) fit$lower[, xi] <- agg$lower fit$larger[, xi] <- agg$larger @@ -168,18 +172,16 @@ QAPglm <- function(formula, } } + # ---- assemble ------------------------------------------------------------- + + # Lift what a reader of the result needs out of the baseline fit, so that + # `fit$coefficients` works without reaching into `fit$base`. fit$coefficients <- fit$base$coefficients fit$t <- fit$base$t - if (!is.null(fit$base$r.squared)) { - fit$r.squared <- fit$base$r.squared - fit$adj.r.squared <- fit$base$adj.r.squared + for (el in c("r.squared", "adj.r.squared", "random.intercepts", + "theta", "zi_coefficients")) { + if (!is.null(fit$base[[el]])) fit[[el]] <- fit$base[[el]] } - if (!is.null(fit$base$random.intercepts)) - fit$random.intercepts <- fit$base$random.intercepts - if (!is.null(fit$base$theta)) - fit$theta <- fit$base$theta - if (!is.null(fit$base$zi_coefficients)) - fit$zi_coefficients <- fit$base$zi_coefficients if (!less_mem) fit$simple_fit <- fit$base$base_model if (family == "binomial") { @@ -194,102 +196,89 @@ QAPglm <- function(formula, fit$diag <- diag fit$family <- family fit$directed <- directed - fit$times <- times + fit$times <- times fit$groups <- unique(unlist(groups)) fit$robust_se <- use_robust_errors + fit$random <- requested[names(slots)] fit$pred <- pred fit$dep <- dep - if (family == "gaussian") { - class(fit) <- "QAPRegression" + + class(fit) <- if (css) { + "QAPCSS" + } else if (family == "gaussian") { + "QAPRegression" } else { - class(fit) <- "QAPGLM" + "QAPGLM" } - return(fit) + fit } +# One permutation: draw, vectorise, refit, and compare against the baseline. +# Returns NULL where the fit fails, which `aggregate_perm_results()` counts. #' @keywords internal #' @noRd -QAPglmPermEst <- function(i, - matlist., - perm_var., - directed., - diag., - mod., - groups., - fit., - family., - use_robust_errors., - has_random., - main_vars., - data_vars., - parsed.) { - - dep <- parsed.$dependent - large <- is.list(matlist.[[dep]]) - - d <- matlist. - if (is.null(perm_var.)) { - if (!large) { - d[[dep]] <- RMPerm(d[[dep]], groups.) +QAPPermEst <- function(i, + matlist., + perm_var., + shape., + directed., + diag., + mod., + groups., + fit., + family., + use_robust_errors., + has_random., + main_vars., + data_vars., + parsed.) { + + dep <- parsed.$dependent + large <- is.list(matlist.[[dep]]) + target <- if (is.null(perm_var.)) dep else perm_var. + + trial <- 0L + repeat { + trial <- trial + 1L + d <- matlist. + d[[target]] <- if (large) { + lapply(d[[target]], shape.$permute, groups = groups.) } else { - d[[dep]] <- lapply(d[[dep]], RMPerm, groups = groups.) + shape.$permute(d[[target]], groups.) } - } else { - if (!large) { - d[[perm_var.]] <- RMPerm(d[[perm_var.]], groups.) - } else { - d[[perm_var.]] <- lapply(d[[perm_var.]], RMPerm, groups = groups.) - } - } - if (!large) { - pred <- make_qap_data(y = d[[dep]], - x = d[data_vars.], - g = groups., - diag = diag., - directed = directed., - net = 1, - perm = FALSE, - xi = NULL) - } else { - pred_list <- vector("list", length(d[[dep]])) - for (net in seq_along(d[[dep]])) { - x2 <- lapply(data_vars., function(v) d[[v]][[net]]) - names(x2) <- data_vars. - g2 <- if (!is.null(groups.)) groups.[[net]] else NULL - pred_list[[net]] <- make_qap_data(y = d[[dep]][[net]], - x = x2, - g = g2, - diag = diag., - directed = directed., - net = net, - perm = FALSE, - xi = NULL) + vec <- .vectorise_matlist(shape., d, dep, data_vars., groups., + diag., directed., large) + pred <- vec$pred + names(pred)[names(pred) == "yv"] <- dep + + if (.sufficient_data(pred, dep, data_vars.)) break + if (trial >= shape.$max_trials) { + # A shape that takes the first draw lets the fit fail and be counted. + if (shape.$max_trials == 1L) break + manynet::snet_abort( + c("Cannot find a valid permutation after {trial} trials.", + i = "The network may be too sparse, or too many cells may be missing.")) } - pred <- do.call(rbind, pred_list) } - names(pred)[names(pred) == "yv"] <- dep - - xi_arg <- if (!is.null(perm_var.)) perm_var. else NULL - # A fit inside the permutation loop runs `times` times, so a fitter's # convergence warning would print once per draw and drown the console. # The count of draws that failed outright is reported by # `aggregate_perm_results()`, which is the number the user needs. perm_fit <- tryCatch( suppressWarnings(fit_qap_model(mod = mod., - pred = pred, - family = family., - use_robust_errors = use_robust_errors., - main_vars = main_vars., - has_random = has_random.)), + pred = pred, + family = family., + use_robust_errors = use_robust_errors., + main_vars = main_vars., + has_random = has_random.)), error = function(e) NULL ) if (is.null(perm_fit)) return(NULL) - return(compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, - fit., xi = xi_arg)) + compare_perm_to_baseline(perm_fit$coefficients, perm_fit$t, + fit., xi = perm_var.) } diff --git a/R/qap_shapes.R b/R/qap_shapes.R new file mode 100644 index 0000000..2305591 --- /dev/null +++ b/R/qap_shapes.R @@ -0,0 +1,106 @@ +# Data shapes ----------------------------------------------------------------- +# +# One engine fits both a dyadic network and a cognitive social structure. The +# two differ in four places and nowhere else: how a network becomes one row per +# observation, how a permutation is drawn, how a residualised predictor is put +# back, and which random intercepts exist. A shape holds those four, so that +# `QAPengine()` can be written once. +# +# Before this, the two were separate 200-line functions that were 55% the same +# code, and every fix had to be made twice. + +#' @keywords internal +#' @noRd +.qap_shape <- function(css = FALSE) { + if (css) .qap_shape_cognitive() else .qap_shape_dyadic() +} + +#' @keywords internal +#' @noRd +.qap_shape_dyadic <- function() { + list( + css = FALSE, + label = "network", + ndim = 2L, + # A sender and a receiver intercept, and one per network in a pooled list. + rand_slots = c(nets = "nv", sender = "sv", receiver = "rv"), + # The dyadic path takes the first draw. A degenerate permutation is dropped + # by the loop's own `tryCatch()`, and `aggregate_perm_results()` counts it. + max_trials = 1L, + vectorise = function(y, x, groups, diag, directed, net) { + pred <- make_qap_data(y = y, x = x, g = groups, diag = diag, + directed = directed, net = net) + list(pred = pred, valid = NULL) + }, + permute = function(m, groups) RMPerm(m, groups), + unresidualise = function(xR, original, pred, large, valid, valid_list) { + residuals_to_matrix(xR, original, pred, large) + } + ) +} + +#' @keywords internal +#' @noRd +.qap_shape_cognitive <- function() { + list( + css = TRUE, + label = "cognitive social structure", + ndim = 3L, + # A perceiver intercept as well, since each perceiver reports the whole + # network and their reports are not independent of each other. + rand_slots = c(nets = "nv", perceiver = "pv", sender = "sv", receiver = "rv"), + # A CSS array is sparse, so a permutation can leave an outcome with one + # value and nothing to fit. Redraw rather than discard: discarding biases + # the null towards the draws that happened to be dense. + max_trials = 10000L, + vectorise = function(y, x, groups, diag, directed, net) { + make_css_data(y = y, x = x, nets = net, diag = diag, directed = directed) + }, + permute = function(m, groups) RMPerm(m, groups, CSS = TRUE), + unresidualise = function(xR, original, pred, large, valid, valid_list) { + residuals_to_array(xR, original, valid, pred, large, valid_list) + } + ) +} + +# Applies a shape's vectoriser to one network or to a pooled list of them. +# Both engines wrote this loop out twice, once for each shape. +#' @keywords internal +#' @noRd +.vectorise_matlist <- function(shape, matlist, dep, data_vars, groups, + diag, directed, large) { + if (!large) { + x <- lapply(data_vars, function(v) matlist[[v]]) + names(x) <- data_vars + out <- shape$vectorise(y = matlist[[dep]], x = x, groups = groups, + diag = diag, directed = directed, net = 1) + return(list(pred = out$pred, valid = out$valid, valid_list = NULL)) + } + + n_nets <- length(matlist[[dep]]) + preds <- vector("list", n_nets) + valids <- vector("list", n_nets) + for (net in seq_len(n_nets)) { + x <- lapply(data_vars, function(v) matlist[[v]][[net]]) + names(x) <- data_vars + g <- if (!is.null(groups) && is.list(groups)) groups[[net]] else groups + out <- shape$vectorise(y = matlist[[dep]][[net]], x = x, groups = g, + diag = diag, directed = directed, net = net) + preds[[net]] <- out$pred + valids[[net]] <- out$valid + } + list(pred = do.call(rbind, preds), valid = NULL, valid_list = valids) +} + +# Whether a permuted sample can be fitted at all: the outcome must vary, and so +# must every numeric predictor. A constant column has no slope to estimate. +#' @keywords internal +#' @noRd +.sufficient_data <- function(pred, dep, data_vars) { + if (length(stats::na.omit(unique(pred[[dep]]))) <= 1) return(FALSE) + numeric_preds <- pred[, intersect(data_vars, names(pred)), drop = FALSE] + numeric_preds <- numeric_preds[, vapply(numeric_preds, is.numeric, logical(1)), + drop = FALSE] + if (ncol(numeric_preds) == 0) return(TRUE) + all(vapply(numeric_preds, function(col) length(unique(col)) > 1, logical(1))) +} diff --git a/tests/testthat/test-qap_shape_css.R b/tests/testthat/test-qap_shape_css.R new file mode 100644 index 0000000..3728f0c --- /dev/null +++ b/tests/testthat/test-qap_shape_css.R @@ -0,0 +1,89 @@ +# The engine fits two shapes, a dyadic network and a cognitive social structure, +# through one skeleton. These assert that the CSS half of that skeleton works, +# since the CSS entry point in `net_regression()` does not exist yet and the +# merge would otherwise be unverified on the shape it was merged for. + +css_fixture <- function(n = 8, seed = 7) { + set.seed(seed) + # One true network, which every perceiver sees with some noise. + truth <- matrix(stats::rbinom(n^2, 1, 0.35), n, n) + diag(truth) <- 0 + y <- array(0L, dim = c(n, n, n)) + for (p in seq_len(n)) { + y[, , p] <- pmin(truth + matrix(stats::rbinom(n^2, 1, 0.12), n, n), 1) + diag(y[, , p]) <- 0 + } + # A dyadic covariate, the same for every perceiver. + dyadic <- outer(stats::runif(n), stats::runif(n), "+") + x <- array(rep(as.vector(dyadic), n), dim = c(n, n, n)) + list(ties = y, cov = x) +} + +test_that("the engine fits a CSS and counts its observations", { + n <- 8 + fit <- QAPengine(ties ~ cov, css_fixture(n), css = TRUE, + family = "binomial", times = 20, seed = 1) + expect_s3_class(fit, "QAPCSS") + expect_named(fit$coefficients, c("(Intercept)", "cov")) + # Every cell of every perceiver's report, less the diagonals. + expect_equal(nrow(fit$pred), n * n * n - n * n) + expect_equal(dim(fit$lower), c(2L, 2L)) + expect_equal(rownames(fit$lower), c("perm_coefs", "perm_t")) +}) + +test_that("a CSS baseline matches glm() on the same observations", { + ml <- css_fixture() + fit <- QAPengine(ties ~ cov, ml, css = TRUE, family = "binomial", + times = 10, seed = 1) + ref <- stats::glm(ties ~ cov, data = fit$pred, family = stats::binomial()) + expect_equal(unname(fit$coefficients), unname(stats::coef(ref))) +}) + +test_that("a CSS carries sender, receiver, perceiver and network indices", { + fit <- QAPengine(ties ~ cov, css_fixture(), css = TRUE, + family = "binomial", times = 10, seed = 1) + expect_true(all(c("sv", "rv", "pv", "nv") %in% names(fit$pred))) + expect_s3_class(fit$pred$pv, "factor") +}) + +test_that("a perceiver random intercept is available to a CSS only", { + skip_if_not_installed("lme4") + fit <- suppressMessages(suppressWarnings( + QAPengine(ties ~ cov, css_fixture(), css = TRUE, family = "binomial", + times = 10, seed = 1, random_intercept_perceiver = TRUE))) + expect_named(fit$random.intercepts, "pv") + expect_true(fit$random[["perceiver"]]) +}) + +test_that("a dyadic network rejects a perceiver random intercept", { + m <- matrix(stats::rnorm(64), 8, 8) + diag(m) <- 0 + expect_error( + QAPengine(ties ~ cov, list(ties = m, cov = matrix(stats::rnorm(64), 8, 8)), + times = 5, random_intercept_perceiver = TRUE), + "perceiver") +}) + +test_that("a CSS rejects an outcome that is not three-dimensional", { + m <- matrix(stats::rnorm(64), 8, 8) + expect_error( + QAPengine(ties ~ cov, list(ties = m, cov = m), css = TRUE, times = 5), + "3-dimensional") +}) + +test_that("a CSS fit prints as a model", { + fit <- QAPengine(ties ~ cov, css_fixture(), css = TRUE, + family = "binomial", times = 10, seed = 1) + expect_output(print(fit), "CSS") + expect_output(print(fit), "Coefficients") + expect_invisible(print(fit)) +}) + +test_that("both shapes reproduce under the same seed", { + ml <- css_fixture() + a <- QAPengine(ties ~ cov, ml, css = TRUE, family = "binomial", + times = 20, seed = 42) + b <- QAPengine(ties ~ cov, ml, css = TRUE, family = "binomial", + times = 20, seed = 42) + expect_equal(a$lower, b$lower) +}) From 5824f77472baaae0df4e1ab7abd7176f68548c55 Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sun, 6 Sep 2026 08:21:57 +0200 Subject: [PATCH 10/12] Simplified news --- NEWS.md | 75 ++++++++++++++++++++------------------------------------- 1 file changed, 26 insertions(+), 49 deletions(-) diff --git a/NEWS.md b/NEWS.md index ecb6a71..c557051 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,38 +2,35 @@ ## Package -- Merged the two engines into one, `QAPengine()` - - `QAPglm()` and `QAPcss()` were 55% the same code, so every fix had to be - made twice; one of them was made in only one place - - What the two shapes do differently is now four functions in - `R/qap_shapes.R`: how to vectorise, how to permute, how to put a - residualised predictor back, and which random intercepts exist - - A random intercept a shape does not have now aborts by name, so a - perceiver intercept on a dyadic network says so - - The engine files fall from 791 lines to 552, with no duplication left -- Branched off five model extensions, to settle the architecture first - - Each is on its own `feature/*` branch, and each strip is one commit that - `git revert` reinstates - - `Suggests` falls from eight modelling packages to three - - See the Github issues for the order they come back in -- Renamed the engine's vocabulary to the front end's, so one word means one - thing on both sides of the seam - - `reps` is now `times`, everywhere including on the returned fit - - `mode` is now `directed`, a logical, and `"digraph"`/`"graph"` are gone; - `mode` is reserved for a nodeset, as in one-mode and two-mode - - `nullhyp` is now `permute`, and its values name what is shuffled: - `"predictor"` for Dekker's double semi-partialling, `"outcome"` for - permuting the dependent matrix alone - - `data` is retired as an identifier: it named the network in one half of - `R/model_regression.R` and the matrix list in the other, one letter away - from `.data` - - `matlist` is the named list of matrices the engine fits - - `net` is one coerced network, inside the formula front end - - `.data` remains the network the user passes in +- Branching off five model extensions reduces `Suggests` packages from eight to three - Updated CONTRIBUTING with the vocabulary table and the reporting rule ## Regression +- Fixed `net_regression()` failing on a two-mode network with more columns than + rows (closing #4) + - Validity was built rows-by-rows, so wider predictor extended it with `NA` + and the dyad count came back as `NA` + - Reported 448 by 12489 network now fits on all 5,595,072 dyads +- Standardised vocabulary in the engine to match the front end: + - Renamed `reps=` to `times=` including on the returned fit + - Renamed `method=`/`nullhyp=` to `permute=`, as method can be ambiguous + - `method = "qap"`/`nullhyp = "qapspp"` is now `permute = "predictor"` + - `method = "qapy"` is now `permute = "outcome"` + - Renamed `mode=` to `directed=`, reserving mode for one-mode and two-mode networks + - `mode = "undirected"` is now `directed = FALSE` + - `data` is retired as potentially confusing: + - `.data` remains the network the user passes in + - `matlist` is the named list of matrices the engine fits + - `net` is one coerced network, inside the formula front end +- Added `snet_info()` reporting of every default the model resolves for itself + - Family chosen from the outcome's values + - Directedness from the network + - `permute = "predictor"` falling back to `"outcome"` with one predictor +- Merged `QAPglm()` and `QAPcss()` engines into one, `QAPengine()` + - 55% the same code, reduces code from 791 lines to 552 + - Differences in treatment are now four functions: vectorisation, permutation, + returning residuals, and identifying random intercepts - Removed the `torch` GPU path (`feature/torch-gpu`) - Gaussian only, duplicated for CSS, no test, and no hosted runner has a CUDA device; `{torch}` in Suggests broke the CI build @@ -50,32 +47,12 @@ - Removed the `fixest_se_cluster` control (`feature/fixest-fixed-effects`) - A bar in the formula now means an `{lme4}` random-effect term, and nothing else; `parse_qap_formula()` drops from three branches to one -- Fixed `net_regression()` failing on a two-mode network with more columns than - rows (closing #4) - - The validity mask was built as rows-by-rows, so a wider predictor extended - it with `NA` and the dyad count came back as `NA` - - The reported 448 by 12489 network now fits, on all 5,595,072 dyads -- Renamed the `method` control to `permute` - - `method = "qap"` is now `permute = "predictor"`, and `method = "qapy"` is - now `permute = "outcome"` -- Renamed the `mode` control to `directed` - - `mode = "undirected"` is now `directed = FALSE` -- Added reporting of every default the model resolves for itself - - The family chosen from the outcome's values - - The directedness read from the network - - `permute = "predictor"` falling back to `"outcome"` with one predictor - - These use `snet_info()`, so `options(snet_verbosity = "verbose")` shows them ## Tests -- Added `test-qap_shape_css.R`, which fits a cognitive social structure through - the merged engine - - `net_regression()` has no CSS entry point yet, so the merge would - otherwise be untested on the shape it was merged for - Added a wide two-mode fixture and two regression tests for #4 - Added `test-qap_reporting.R`, which runs with `snet_verbosity = "verbose"` - - Informational output is silent in every other test, so a message that - `{cli}` cannot parse was invisible until it aborted; two shipped that way +- Added `test-qap_shape_css.R`, which fits a cognitive social structure # infernet 0.1.1 From 13b1a5f82e2f96fb5f4ad8675817a4721c68723d Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sun, 6 Sep 2026 08:23:28 +0200 Subject: [PATCH 11/12] Fixed closes syntax --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index c557051..a3d2412 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ ## Regression - Fixed `net_regression()` failing on a two-mode network with more columns than - rows (closing #4) + rows (closes #4) - Validity was built rows-by-rows, so wider predictor extended it with `NA` and the dyad count came back as `NA` - Reported 448 by 12489 network now fits on all 5,595,072 dyads From d655b35bca5440d945b0527c2fe38bb2c3cf8f0a Mon Sep 17 00:00:00 2001 From: James Hollway Date: Sun, 6 Sep 2026 08:49:33 +0200 Subject: [PATCH 12/12] Fixed `groups=` being refused on a two-mode network --- NEWS.md | 2 ++ R/model_regression.R | 15 --------------- R/qap_engine.R | 12 +++++++++--- tests/testthat/test-qap_shapes.R | 20 ++++++++++++++++++++ 4 files changed, 31 insertions(+), 18 deletions(-) diff --git a/NEWS.md b/NEWS.md index a3d2412..010385c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,8 @@ - Validity was built rows-by-rows, so wider predictor extended it with `NA` and the dyad count came back as `NA` - Reported 448 by 12489 network now fits on all 5,595,072 dyads +- Fixed `groups=` being refused on a two-mode network unless it matched the row + mode, though either mode may be the one that is blocked - Standardised vocabulary in the engine to match the front end: - Renamed `reps=` to `times=` including on the returned fit - Renamed `method=`/`nullhyp=` to `permute=`, as method can be ambiguous diff --git a/R/model_regression.R b/R/model_regression.R index 59d70da..72840b7 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -525,21 +525,6 @@ print.net_regression <- function(x, ..., } -.print_glm_table <- function(base, lower, larger, abs_mat, permute, print_b) { - cat("\n\nCoefficients:\n") - nc <- length(base$coefficients) - cmat <- matrix(NA, nrow = nc, ncol = 4) - cmat[, 1] <- format(round(as.numeric(base$coefficients), 4)) - cmat[, 2] <- format(lower[2, ]) - cmat[, 3] <- format(larger[2, ]) - cmat[, 4] <- format(abs_mat[2, ]) - if (permute == "predictor") cmat[1, 2:4] <- "*" - colnames(cmat) <- c("Estimate", "Pr(<=t)", "Pr(>=t)", "Pr(>=|t|)") - rownames(cmat) <- names(base$coefficients) - print.table(cmat) -} - - # ============================================================================ # Formula -> matrix-list front end # ============================================================================ diff --git a/R/qap_engine.R b/R/qap_engine.R index c853c1f..6b488e1 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -84,10 +84,16 @@ QAPengine <- function(formula, # ---- vectorise ------------------------------------------------------------ if (!large && !is.null(groups)) { - n <- dim(matlist[[dep]])[1] - if (length(groups) != n) + # A blocking factor names the nodes of one mode. A one-mode network has a + # single mode, but a two-mode network has two of different sizes, and either + # may be the one that is blocked. `.perm_order()` already permutes a mode + # freely when the factor cannot describe it, so accept a length that matches + # any side of the outcome and abort only when it matches none. + ns <- unique(dim(matlist[[dep]])) + if (!length(groups) %in% ns) manynet::snet_abort( - "{.arg groups} is of length {length(groups)}, but the network has {n} nodes.") + c("{.arg groups} is of length {length(groups)}.", + i = "It must match one of the outcome's node sets: {ns}.")) groups <- as.factor(groups) } diff --git a/tests/testthat/test-qap_shapes.R b/tests/testthat/test-qap_shapes.R index f208156..ea7ebbb 100644 --- a/tests/testthat/test-qap_shapes.R +++ b/tests/testthat/test-qap_shapes.R @@ -116,3 +116,23 @@ test_that("a missing dyad is dropped from the model", { control = list(seed = 1)) expect_equal(nrow(fit$pred), 20 * 19 - 1) }) + +# A blocking factor names the nodes of one mode. A two-mode network has two +# modes of different sizes, so a factor of the row length cannot also be of the +# column length. Requiring the row length refused a factor that blocks the +# columns, which `.perm_order()` handles. See the review of stocnet/infernet#13. +test_that("a two-mode network accepts a grouping factor for either mode", { + g <- qap_net_twomode_wide(nr = 12, nc = 40) + rows <- rep(c("a", "b"), length.out = 12) + cols <- rep(c("a", "b", "c", "d"), length.out = 40) + for (grp in list(rows, cols)) { + fit <- net_regression(. ~ ego(Att), g, times = 5, + control = list(seed = 1, groups = grp)) + expect_s3_class(fit, "net_regression") + expect_equal(nrow(fit$pred), 12 * 40) + } + expect_error( + net_regression(. ~ ego(Att), g, times = 5, + control = list(seed = 1, groups = rep("a", 7))), + "length 7") +})