diff --git a/.Rbuildignore b/.Rbuildignore index fb10167..c8bddc1 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -18,3 +18,8 @@ vignettes/precompile\.R ^CRAN-SUBMISSION$ ^README\.Rmd$ ^.mailmap$ +^\.positai$ +^\.claude$ +^CLAUDE\.md$ +^\.Rhistory$ +^\.RData$ diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index 5e6cf14..859486f 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -4,7 +4,11 @@ Contributions to `infernet`, whether in the form of issue identification, bug fixes, new code or documentation are encouraged and welcome. -## Git and Bitbucket +Please note that the `infernet` project is released with a +[Contributor Code of Conduct](CODE_OF_CONDUCT.md). +By contributing to this project, you agree to abide by its terms. + +## Git `stocnet` projects are maintained using the git version control system. A plain-English introduction to git can be found [here](https://blog.red-badger.com/2016/11/29/gitgithub-in-plain-english). @@ -18,32 +22,38 @@ but I recommend [Fork](https://git-fork.com) software for Mac and Windows. This allows mostly visual management of commits, diffs, branches, etc. There are various other git software packages available, but this one is fairly fully featured. -The Github page allows to access the issues assigned to you and check the commits. +The GitHub page allows to access the issues assigned to you and check the commits. You can also access the documents in the repository, although this won't be necessary after you have cloned it on your computer via Fork. -## Style +### Identifying issues -In terms of style, we are aiming for pleasant predictability in terms of user experience. -To that end, we have a regular syntax that users can rely on producing expected effects. +Please use the issues tracker on GitHub to identify any function-related issues. +You can use these issues to track progress on the issue and +to comment or continue a conversation on that issue. +The most useful issues are ones that precisely identify an error, +or propose a test that should pass but instead fails. +Examples for documentation are also most welcome. -## Fork +Issues that belong to another package in the family should be transferred there +rather than fixed here: `gh issue transfer stocnet/`. +See the division of labour below. ### Cloning + Once you have downloaded Fork, the first thing you have to do is to clone the remote repository on your computer. Before cloning, you will be able to choose on which `branch` you want to work: develop or main. -### Pull -This command allows you to `pull` changes from the remote repository to your local repository on Sourcetree. +### Pull + +This command allows you to `pull` changes from the remote repository to your local repository. Make sure you do that before starting working on your files so you have the newest versions. -When pulling, make sure you choose master or develop, +When pulling, make sure you choose main or develop, depending on the branch you decided to work with. Once you pulled, you have now all the new commits and files and you can start working on your assigned tasks. -Note that you can access and open the files either from the Finder or from Fork. -Some documents might be stored using Large File Storage (LFS) to save space on the repository. ### Commit and Push @@ -51,55 +61,529 @@ Once you have made modifications on a file and saved them, it will appear in you Here you can control one last time your file, write the commit message with the issue reference (see below) and commit. Once your commit is ready, you can `push` them to the origin/main repository. -Note that you can click the "push immediately" box in the commit window -if you don't want to do it in two steps. If you are working on a separate branch, it is important to select this branch when pushing to origin/main. -## Issues and tests +Commits may reference an existing GitHub issue number. +Where the issue number is preceded by `resolve`/`resolves`/`resolved`, +`close`/`closes`/`closed`, or `fix`/`fixes`/`fixed` (capitalised or not), +GitHub updates the status of the issue automatically. -Please use the issues tracker on Github to identify any function-related issues. -You can use these issues to track progress on the issue and -to comment or continue a conversation on that issue. -Currently issue tracking is only open to those involved in the project. +### Branching and CI + +- `main` is the release branch; `develop` is the working branch (clone/work on `develop`). +- PRs into `main` trigger [prchecks.yml](workflows/prchecks.yml): R CMD check + (macOS/Windows/Linux), binary build, codecov, lintr, spell check, + and PR metadata checks (DESCRIPTION version bump, PR title/description conventions). +- Merges/pushes to `main` trigger [pushrelease.yml](workflows/pushrelease.yml): + check, auto-bump version tag, GitHub release with binaries and release notes + taken from `NEWS.md`, then pkgdown site deploy. +- The PR metadata job requires that each PR into `main` bumps the `Version:` field in + `DESCRIPTION` by the appropriate increment, names that new version in the PR title, + and itemizes its changes in the PR description under `##` subsection titles matching + the `NEWS.md` conventions below. +- Development dependencies are declared in `DESCRIPTION` under `Config/Needs/build`, + `Config/Needs/check`, and `Config/Needs/website` rather than `Suggests` — + the workflows install them via `needs:` in `setup-r-dependencies`. +- A merge that touches [R/model_regression.R](../R/model_regression.R) or the + `R/qap_*.R` engine files deserves a parse check before it is pushed + (`Rscript -e 'devtools::load_all()'`). + The merge that created the current engine silently dropped three function + headers and left an orphan function body, so the package did not parse at all. + +## Style + +In terms of style, we are aiming for "pleasant predictability" in terms of user experience. +To that end, we have a regular syntax that users can rely on producing expected effects. +Functions in the same family (`test_*()`, etc.) should share +argument order and naming, so that behaviour is guessable across the family. + +We are also aiming for "declarative simplicity", +where functions carry as few arguments as possible to reduce the documentation burden, +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). + +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, +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". + +## Package architecture + +### Project overview + +`infernet` is an R package (part of the [stocnet](https://github.com/stocnet) ecosystem) +providing the *inferential layer* for network analysis: +conditional uniform graph (CUG) and quadratic assignment procedure (QAP) tests of network +statistics, and multiple regression QAP (MRQAP) for network data. +Because it builds on `{manynet}`, every function accepts matrices, edgelists, +`{igraph}`, `{network}`, `{tidygraph}` or `stocnet` objects, +and one-mode or two-mode networks alike. + +`infernet` combines two lines of work: + +- the formula front end, multimodal treatment, + and ease of use of `{migraph}`'s `net_regression()`, and +- the estimator range, missing-data handling, and cognitive social structure (CSS) + treatment of Robert Krause's `MrQAP`, which is where the `R/qap_*.R` engine + files were ported from. + +Division of labour to keep in mind when adding functions: + +- `{manynet}`: network classes/coercion (`as_*()`), making and manipulating networks, + and network-level logical tests (e.g. `is_directed()`, `is_twomode()`). +- `{netrics}`: everything analytic — marks, measures, memberships, motifs — + at the node, tie, and network level. +- `{autograph}`: drawing graphs and plotting analytic, modelling, or diagnostic results, + along with deep (often institutional) theming. *All* plot methods should live there. +- `{infernet}` (this package): testing and modelling, e.g. CUG/QAP/MRQAP. +- `{migraph}`: the software companion to *Multimodal Political Networks*, holding the + `mpn_*` datasets and the diffusion models. + `{migraph}` still carries older copies of `net_regression()` and the `test_*()` family; + those are what `{infernet}` supersedes, so a fix made there needs porting here, + and a fix made here does not need porting back. +- `{goldfish}`: another stocnet package for estimating network event models + such as relational event models and dynamic network actor models. + Note that these two packages, `{infernet}` and `{goldfish}`, + should reuse similar vocabulary and syntax where possible + to improve the "pleasant predictability" of each package. + +### Common commands + +This is a standard R package developed with `devtools`/`roxygen2`. +Run these from an R console with the working directory set to the package root +(or via `Rscript -e`). + +- Load package for interactive development: `devtools::load_all()` +- Regenerate docs & NAMESPACE after editing roxygen comments: `devtools::document()` +- Run full test suite: `devtools::test()` +- Run a single test file: `devtools::test(filter = "net_regression")` + (matches `test-net_regression.R`), or `testthat::test_file("tests/testthat/test-net_regression.R")` +- Full package check (mirrors CI): `devtools::check()` or `rcmdcheck::rcmdcheck()` +- Lint: `lintr::lint_package()` +- Spell check: `spelling::spell_check_package()` +- Code coverage: `covr::package_coverage()` +- Rebuild `README.md` from `README.Rmd`: `devtools::build_readme()` +- Check every topic is in the pkgdown index: `pkgdown::check_pkgdown()` +- Build pkgdown site locally: `pkgdown::build_site()` + +There is no non-R build system — no package.json/Makefile. +Roxygen is configured with `markdown = TRUE`; +`NAMESPACE` and all `man/*.Rd` files are generated — never hand-edit them. +Likewise `README.md` is generated from `README.Rmd` — edit the `.Rmd` and re-knit. + +### File organization + +`R/` files are grouped by theme rather than one file per function. +The user-facing files are named `model_*.R`; +the internal engine ported from `MrQAP` is named `qap_*.R`. + +| File | Contains | +|---|---| +| `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_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_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 | +| `qap_misc.R` | small combining and reshaping helpers | +| `infernet-package.R` | package-level doc, global variables, and the shim silencing R CMD check's unused-import note | +| `zzz.R` | the attach-time greeting and the cached stocnet version check | + +Only `model_*.R` holds exported functions. +Everything in `qap_*.R` is internal, and carries `@keywords internal` and `@noRd`. +Keep it that way: the engine's argument names are the `MrQAP` ones, and +exporting them would freeze an interface we still intend to tidy. + +### The `net_regression()` pipeline + +`net_regression()` ([R/model_regression.R](../R/model_regression.R)) is the single +regression entry point. Its control flow is: + +1. Merge the user's `control` list over `.default_control()`, + so an unnamed entry falls back to the default rather than to `NULL`. +2. Dispatch on the input shape. + A single network goes to `convertToMatrixList()`; + a list of networks goes to `.prepare_list_of_graphs()`, which converts each + 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 + `diag` from the network with `manynet::is_directed()` and `manynet::is_complex()`. +4. Call `QAPglm()`, which parses the formula, fits the baseline model once + via `fit_qap_model()`, then runs `reps` 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. +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. + +The formula front end accepts these terms, and a new one should be added +to `getRHSNames()` and `convertToMatrixList()` together: + +- `ego(attr)` — the sender's value of a nodal attribute, +- `alter(attr)` — the receiver's value, +- `same(attr)` — 1 where sender and receiver share an attribute value, +- `dist(attr)` — the absolute difference in a numeric attribute, +- `sim(attr)` — the proportional similarity in a numeric attribute, +- `tertius(attr, fn)` — an aggregate of an attribute over a node's other ties, +- a plain name — another network, used as a dyadic covariate. + +A formula is meant to be reusable across models, +so a term must mean the same thing whichever family or engine consumes it. +Where a term cannot apply, say so through `snet_abort()` rather than +silently dropping it. + +### Function body conventions + +Test functions consistently: + +1. Compute the observed statistic by applying the user-supplied `FUN` to `.data`. +2. Generate `times` random, configuration-preserving, or permuted networks + via the corresponding `{manynet}` generator (`generate_random()`, + `generate_configuration()`, `to_permuted()`), rebinding node attributes with + `manynet::bind_node_attributes()` where the statistic needs them. +3. Recompute `FUN` over each simulated network. +4. Return a `network_test` object recording the test type, observed value, + simulated distribution, one- and two-tailed p-values, + and the network's properties (`is_directed()`, `is_complex()`). + +Properties of the dependent network — modes, directedness, loops — must always be +respected in permutations and analysis; +one-mode and two-mode cases are branched on explicitly rather than projected away. +All `manynet`, `netrics`, `furrr`, and `future` calls use explicit `::` namespacing +(with per-file `@importFrom` roxygen tags for NAMESPACE generation). +Plot methods belong in `{autograph}`, not here. + +Every `print.*` method returns `invisible(x)`. +A print method that returns the result of its last `cat()` prints `NULL` when its +value is used, which is what `print.network_test()` did before this was fixed. + +### Parallelism + +Every simulation-heavy function takes: + +- `times` — the number of simulations (default `1000`; 1,000–10,000 for publication). +- `strategy` — a `{future}` plan name (default `"sequential"`; + `"multisession"`/`"multicore"` for multiple cores), + set with `future::plan(strategy)` and restored via `on.exit()`. + In `net_regression()` this is passed through the `control` list. + +The `test_*()` family maps simulations with `furrr::future_map*()` using +`furrr::furrr_options(seed = TRUE)` for reproducible parallel RNG. +The regression engine wraps the same machinery in `setup_future_plan()` and +`run_permutations()` ([R/qap_utils.R](../R/qap_utils.R)); use those rather than +calling `future::plan()` from a new engine function. + +Progress reporting is not a separate argument. +It is read from `options(snet_verbosity)`, so that one option governs how +talkative every stocnet package is. + +### Input shapes + +Most defects reported against this package are not wrong arithmetic. +They are a network shape the function did not expect. +Before finishing a function, run it on a signed, a weighted, a directed, a two-mode, +a multiplex, a multilevel and a longitudinal network, and decide each case deliberately. +Document each decision in the roxygen block with an `@section` named for the shape, +e.g. `@section Two-mode networks:`. + +Missing data deserves particular care here, since handling it well is a reason +this package exists. +An unobserved dyad is `NA`, never a zero tie and never a source's numeric code. +`make_qap_data()` drops NA and diagonal cells before fitting, so a predictor +that introduces `NA` silently shrinks the sample: +report what was dropped rather than letting the count change without comment. + +### Console messaging + +All user-facing messages go through the `snet_*()` wrappers exported by `{manynet}`, +rather than base `message()`/`stop()`/`warning()` or `{cli}` calls directly: + +| Wrapper | Use for | +|---|---| +| `snet_abort()` | errors: the function cannot proceed | +| `snet_warn()` | the function proceeds, but the user should know something | +| `snet_info()` | notable information about what was done, e.g. a defaulted argument or the method dispatched to | +| `snet_minor_info()` | incidental detail | +| `snet_success()` | confirmation that a requested operation completed | +| `snet_unavailable()` | not-yet-implemented features | + +Every wrapper except `snet_abort()` is silenced by +`options(snet_verbosity = "quiet")`, which is the *default* — +so informational output must never be load-bearing, +and errors must carry everything the user needs to act. +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.")`. +- 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. +- Use `{cli}`'s pluralisation rather than hand-written branches: + `snet_warn("Dropped {length(dropped)} network{?s}.")`. + +Prefer "`{.arg times}` must be a positive whole number" over "invalid input". +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\")}."))`. + +The model specification advice printed by `specificationAdvice()` is +informational, not load-bearing, so it belongs at `snet_info()`. + +### Dependencies + +`infernet` `Depends` on `{manynet}` (network classes, coercion and logical tests) +and `{netrics}` (measures), and `Imports` the `{future}`/`{furrr}`/`{purrr}` +stack plus `{reformulas}` for formula surgery. + +Everything that only one estimator needs is in `Suggests`: +`lme4` and `glmmTMB` (random effects), `fixest` (fixed effects), +`gmm` (GMM estimation), `MASS` (negative binomial), `pscl` (zero-inflated Poisson), +`nnet` (multinomial), and `torch` (the GPU path). +A code path that depends on a suggested package must guard with +`requireNamespace()` and abort with an actionable message, +and its tests must `skip_if_not_installed()`. +Keeping these optional is deliberate: the common case — a gaussian or binomial +MRQAP — must install and run with no compiler and no heavy dependency tree. + +The declared minimum of each `stocnet` dependency is the version on CRAN, +so that CI can install it. +Where `infernet` needs something that only a newer, unreleased `{manynet}` has, +reach it through a shim rather than by raising the minimum, +and resolve the name at call time from the namespace. +Test for the function rather than for the version string, +because a pre-release development build can carry the version +without yet exporting the function. + +### Tests -The most useful issues are ones that precisely identify an error, -or propose a test that should pass but instead fails. This package uses the `testthat` package for testing functions. Please see the [testthat website](https://testthat.r-lib.org) for more details. +`testthat` edition 3 with parallel execution is configured in `DESCRIPTION` +(`Config/testthat/parallel: true`). +`Config/testthat/start-first` should prioritise the test files that take longest to run. -## Bug fixing or adding new code +Tests in `tests/testthat/` mirror the `R/` files for the exported functions +(`test-net_regression.R`, `test-model_tests.R`), +and are grouped by contract for the engine: -Independent or assigned code contributions are most welcome. -When writing new code, please follow -[standard R guidelines](https://www.r-bloggers.com/🖊-r-coding-style-guide/). -It can help to use packages such as `lintr`, `goodpractice` and `formatR` -to ensure these are followed. +| File | Asserts | +|---|---| +| `test-qap_estimators.R` | each `family` and `estimator` combination, against the equivalent standard fit | +| `test-qap_shapes.R` | the dyads that reach the model, for each shape of network | +| `test-qap_reproducibility.R` | that a seed reproduces a run, sequentially and in parallel | +| `test-qap_control.R` | the `control` list and the choice of null hypothesis | -Currently, commits can only be pushed to Bitbucket where they reference an existing issue. -If no issue exists for the code you have developed, please add an issue first before pushing. -Once the issue exists, you will need to mention the issue number (preceded by a hash symbol: #) -in the commit description: +[tests/testthat/helper-infernet.R](../tests/testthat/helper-infernet.R) +holds the shared fixtures and expectations: -``` Resolved #31 by adding a new function that does things, also updated documentation ``` +- `qap_net_gaussian()`, `qap_net_binary()`, `qap_net_count()`, `qap_net_zip()`, + `qap_net_undirected()`, `qap_net_twomode()` — each seeded, and each carrying + real signal, so that every family converges and the comparison is not testing + noise against noise. +- `qap_reference_data()` — rebuilds the dyad-level data frame the engine fits, + so a baseline coefficient can be compared against `lm()`, `glm()`, + `MASS::glm.nb()`, `pscl::zeroinfl()`, `lme4::lmer()` or `fixest::feglm()` + on identical data. +- `expect_qap_shape()` — the shape contract every estimator meets, whatever it + fits underneath: named coefficients, and `lower`/`larger`/`abs` as + two-row matrices of proportions with matching dimnames. -Where the issue hash (i.e. #31) is preceded by -`resolve`, `resolves`, `resolved`, `close`, `closes`, `closed`, `fix`, `fixes`, or `fixed` -(capitalised or not), -Github will automatically updated the status of the issue(s) mentioned. +Four things are worth asserting for every estimator that is added: -Our current syntactical standard is to mention the issue first and then -provide a short description of what the committed changes do -in relation to that issue. -Any ancillary changes can be mentioned after a comma. +1. That the baseline coefficients match those of the equivalent standard fit on + the same dyad-level data, using `qap_reference_data()`. + The permutation inference is what is novel; the point estimates are not, + and they should agree. +2. That the result meets `expect_qap_shape()`. + Most engine defects found so far surfaced as a name or a dimension, not as a + wrong number: a backticked coefficient name broke double semi-partialling, + and a stray placeholder intercept broke the `{fixest}` path. +3. That a given `seed` reproduces the same p-values. + Permutation results are only comparable across runs if the RNG is, + and `furrr_options(seed = TRUE)` and `future.seed = TRUE` are what make that + true in parallel. +4. That the estimator is reached at all. + Several paths in `fit_qap_model()` are selected by a combination of `family`, + `estimator` and the random/fixed effects flags, + so a test that does not name that combination does not cover it. -## Documentation +An estimator that needs a package from `Suggests` takes +`skip_if_not_installed()`, so the suite still passes where that package is +absent. Do not let a path go untested because the package is missing locally: +install it, and check that the test runs before you rely on the skip. -A final way of contributing to the package is in developing the -vignettes/articles that illustrate the value added in the package. -Please contact me with any proposals here. +Note that `skip_if_not_installed()` is weaker than it looks. CI installs every +`Suggests`, so the skip does not fire there, and an installed package is not +always a working one: `{torch}` installs as an R package before its Lantern +backend is downloaded, and `torch::cuda_is_available()` then throws rather than +returning `FALSE`. Guard on the capability, not on the package. -Please note that the `infernet` project is released with a -[Contributor Code of Conduct](CODE_OF_CONDUCT.md). -By contributing to this project, you agree to abide by its terms. +A test must not depend on a `{manynet}` feature newer than the CRAN version, or +it passes here and fails on CI. + +Note that `options(snet_verbosity)` is unset under `R CMD check`, +because manynet's `.onAttach` only sets it in an interactive session. +Never write a test that depends on `snet_info()` output. + +Count the dyads rather than checking that a call returns. +A directed network of *n* nodes contributes *n*(*n*-1) dyads, an undirected one +*n*(*n*-1)/2, and a two-mode one every cell of its incidence matrix. +Each of those was wrong at some point, and each looked like a working model. + +A fitter's warning raised inside the permutation loop is held back, since it +would print once per draw; `aggregate_perm_results()` reports the number of +draws that failed outright. A test therefore should not expect a convergence +warning from a permutation, only from the baseline. + +The aim is to work towards comprehensive coverage, +so each change should be fully covered by tests. +However, we also need to keep an eye on the clock: +CRAN complains if tests take too long, +so use small fixtures, low `times`, and `skip_on_cran()` for taxing tests. +`# nocov start` and `# nocov end` can be used to exclude lines or functions +that are too difficult to cover. + +### Documentation + +Roxygen is configured with `markdown = TRUE`; +`NAMESPACE` and all `man/*.Rd` files are generated — never hand-edit them. +Run `devtools::document()` after changing any roxygen comment. + +- Related functions share one roxygen block via `@name`/`@rdname`, + matching the file organisation above. +- Document each argument once. + `net_regression()` takes its options through a `control` list, so document + each entry of that list as a bullet under `@param control` rather than as its + own `@param`. A stray `@param` for something that is no longer a formal + argument is an R CMD check warning, and the duplicated blocks that the + engine merge introduced were exactly that. +- Every exported function needs a runnable `@examples` block: + examples are run by R CMD check, and they are also the fastest documentation for users. + Prefer the bundled `ison_*`/`fict_*` networks from `{manynet}` over ad hoc + constructions, and keep `times` small so the example is fast. +- Use the native pipe `|>` in examples, never `%>%`. + `{migraph}` shipped examples that called `%>%` after the re-export was + removed, and every one of them failed R CMD check. +- Cite the source of a method with `@references` in the ecosystem's format + (authors, year, title, journal, and `\doi{}` where available), + so that users can trace an implementation back to its definition. +- Documented behaviour and implemented behaviour must agree. + When you change a default, search the roxygen for it too. + +### README and website + +The README offers a landing page for new users, both on the GitHub repository +as well as on the website. +As such, it should make a compelling case for the value added of the package, +and not drift out of date. +Note that `README.md` is generated from `README.Rmd` — edit `README.Rmd` and re-knit +(`devtools::build_readme()`), never edit `README.md` directly. + +The website is created by pkgdown from [pkgdown/_pkgdown.yml](../pkgdown/_pkgdown.yml), +and is deployed automatically when changes reach `main`. +Please make sure that the pkgdown website will build correctly before opening a PR: + +```r +pkgdown::check_pkgdown() # every topic is in the index +pkgdown::build_site(preview = FALSE) # everything else +``` + +The most common failure is a new exported function that is not picked up under the +function overview (the `reference:` section of `_pkgdown.yml`) — +pkgdown requires *every* exported topic to appear there exactly once, or it will not build. +A helper that users are not meant to call takes `@keywords internal` instead. +These `reference:` titles are also the headings used in `NEWS.md` (see below), +so keep the two in step. + +### `NEWS.md` conventions + +`NEWS.md` groups each version's changes under `##` headings that mirror the website +function overview (`pkgdown/_pkgdown.yml` `reference:` titles). +Lead with `## Package` (package-wide/website/infrastructure changes), +then `## Tests` (the `test_*()` family) and `## Regression` (`net_regression()` +and the engine behind it). +Each heading appears at most once per version. + +Start each bullet with a verb matching the change type: + +- `Added ...` — new functionality +- `Fixed ...` — bug fixes; if it relates to a GitHub issue, suffix with `(closing #123)` +- `Renamed ... to ...` — function or data name migrations +- `Improved ...` — functional updates to existing behaviour +- `Updated ...` — documentation changes + +If a cited GitHub issue was **not** authored by @jhollway, thank the author with an +`@`-tag in the bullet. + +#### Grouping + +Group first, and only then write the bullets. +The more entries a version holds, the more this matters. + +- Cluster related changes as indented sub-bullets under a lead bullet. +- Where several changes concern one function, lead with an `Improved ...` bullet naming + the function, and put the individual `Fixed ...`/`Added ...` points beneath it, + so the cluster groups by function rather than by change type. +- Under such a lead bullet, do not name the function again in the sub-bullets, + since the lead bullet already carries it. +- Where one decision runs across many functions, lead with the decision rather than + with each function. +- Sub-bullets indent by two spaces, and nest at most one level further (four spaces). + +#### Writing the bullets + +`NEWS.md` is read by users scanning for what changed, not by reviewers reading prose, +so each bullet is a headline rather than a sentence, +so avoid over-punctuation or over-explanation. +Details can be added to the function documentation, if necessary. + +- No full stop at the end of a bullet +- Keep every bullet to one line of fewer than 81 characters ideally + (a few more or less is fine) + - If a bullet wraps, it holds too much: shorten it, + or split it into a lead bullet and sub-bullets +- One clause where possible, and at most one comma + - Use a semicolon for a short second clause, e.g. "old spelling still works but warns" + - Use a sub-bullet where the second clause needs more room than that +- Name the function or object in backticks and say what changed to it, + dropping scaffolding like "This change ...", "In order to ...", or "as part of an effort to" +- Keep the *what*, and add the *why* only where the behaviour would otherwise look arbitrary +- No trailing rationale, no restating the same change twice in different words, + and no marketing adjectives such as "comprehensive" or "robust" +- A sub-bullet does not need a verb: it can state the consequence, + the previous behaviour, or an example call +- Cut a sub-bullet that only restates what the lead bullet already implies +- Where several bullets describe parallel changes, reuse the sentence structure, + so that a reader sees the parallelism at a glance +- Use one word for one thing throughout a version's entries, + rather than varying the wording for effect + +For example, instead of: + +> Fixed a bug where, in some cases, `print.network_test()` was not returning its +> input invisibly, which meant that `x <- print(test)` assigned NULL. + +write: + +> Fixed `print.network_test()` to return its input invisibly + +and instead of: + +> Added a new control option, `use_gpu`, which is a useful option that allows the +> permutations to be run in batch on the GPU using torch. + +write: +> Added `use_gpu` control for batch OLS permutations on the GPU diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index b3a6137..6cf530b 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -2,18 +2,7 @@ # Checklist: -- PR form - - [ ] Description above itemizes changes under subtitles, e.g. "## Data"" - - [ ] Any closed, fixed, or related issues are referenced and explained in the description above, e.g. "Fixed #0 by adding A" - - [ ] Package builds on my OS without issues -- PR checks all pass for latest commit - - [ ] CodeFactor check: Package improves or maintains good style - - [ ] Package builds on Mac - - [ ] Package builds on Windows - - [ ] Package builds on Linux - - [ ] CodeCov check: Package improves or maintains good test coverage - Documentation - - [ ] Any new or modified functions or data have roxygen style documentation in their .R scripts - - [ ] Longer functions are commented inline or broken down into helper functions so that it is easier to debug in the future + - [ ] Longer functions are commented inline or broken down into helper functions to help debugging +- PR form - [ ] PR description above and the NEWS.md file are aligned - - [ ] DESCRIPTION file version is bumped by the appropriate increment (major, minor, patch) diff --git a/.github/workflows/prchecks.yml b/.github/workflows/prchecks.yml index 9e541a0..05e3120 100644 --- a/.github/workflows/prchecks.yml +++ b/.github/workflows/prchecks.yml @@ -24,11 +24,13 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v7 - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + use-public-rspm: true - uses: r-lib/actions/setup-pandoc@v2 @@ -40,6 +42,7 @@ jobs: any::rcmdcheck - uses: r-lib/actions/check-r-package@v2 + env: _R_CHECK_FORCE_SUGGESTS_: false with: upload-snapshots: true @@ -53,13 +56,16 @@ jobs: shell: Rscript {0} - name: Save binary artifact - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: ${{ matrix.config.asset_name }} path: build/ - name: Calculate code coverage - run: Rscript -e "covr::codecov()" + if: runner.os == 'macOS' + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + run: Rscript -e 'covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))' - name: Lint run: lintr::lint_package() @@ -68,3 +74,49 @@ jobs: - name: Spell check run: spelling::spell_check_package() shell: Rscript {0} + + pr-metadata: + name: PR metadata checks + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v7 + with: + fetch-depth: 0 + + - name: Check DESCRIPTION version is bumped + env: + BASE_SHA: ${{ github.event.pull_request.base.sha }} + run: | + BASE_VERSION=$(git show "$BASE_SHA:DESCRIPTION" | grep -m1 '^Version:' | sed 's/Version: *//') + HEAD_VERSION=$(grep -m1 '^Version:' DESCRIPTION | sed 's/Version: *//') + echo "Base version: $BASE_VERSION" + echo "Head version: $HEAD_VERSION" + if [ "$BASE_VERSION" = "$HEAD_VERSION" ]; then + echo "::error::DESCRIPTION Version has not been bumped from $BASE_VERSION" + exit 1 + fi + HIGHEST=$(printf '%s\n%s\n' "$BASE_VERSION" "$HEAD_VERSION" | sort -V | tail -n1) + if [ "$HIGHEST" != "$HEAD_VERSION" ]; then + echo "::error::DESCRIPTION Version $HEAD_VERSION is not greater than base version $BASE_VERSION" + exit 1 + fi + + - name: Check PR title mentions the new version number + env: + PR_TITLE: ${{ github.event.pull_request.title }} + run: | + HEAD_VERSION=$(grep -m1 '^Version:' DESCRIPTION | sed 's/Version: *//') + if [[ "$PR_TITLE" != *"$HEAD_VERSION"* ]]; then + echo "::error::PR title does not mention the new version number ($HEAD_VERSION): \"$PR_TITLE\"" + exit 1 + fi + + - name: Check PR description itemizes changes under subsection titles + env: + PR_BODY: ${{ github.event.pull_request.body }} + run: | + if ! grep -qE '^##[[:space:]]' <<< "$PR_BODY"; then + echo "::error::PR description does not itemize changes under subsection titles (e.g. \"## Data\")" + exit 1 + fi diff --git a/.github/workflows/pushrelease.yml b/.github/workflows/pushrelease.yml index 9961552..da6afa5 100644 --- a/.github/workflows/pushrelease.yml +++ b/.github/workflows/pushrelease.yml @@ -24,7 +24,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v7 - uses: r-lib/actions/setup-r@v2 with: @@ -55,14 +55,16 @@ jobs: shell: Rscript {0} - name: Save binary artifact - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: ${{ matrix.config.asset_name }} path: build/ - name: Calculate code coverage - if: runner.os == 'macOS-latest' - run: Rscript -e "covr::codecov()" + if: runner.os == 'macOS' + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + run: Rscript -e 'covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))' release: name: Bump version and release @@ -73,7 +75,7 @@ jobs: contents: write steps: - name: Checkout one - uses: actions/checkout@v4 + uses: actions/checkout@v7 with: fetch-depth: '0' - name: Bump version and push tag @@ -85,7 +87,7 @@ jobs: DEFAULT_BUMP: patch RELEASE_BRANCHES: main - name: Checkout two - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Extract version run: | @@ -93,7 +95,7 @@ jobs: echo "PACKAGE_NAME=$(grep '^Package' DESCRIPTION | sed 's/.*: *//')" >> $GITHUB_ENV - name: Download binaries - uses: actions/download-artifact@v4 + uses: actions/download-artifact@v8 - name: Rename binaries release shell: bash @@ -104,13 +106,27 @@ jobs: cp ./winOS/${{ env.PACKAGE_NAME }}_${{ env.PACKAGE_VERSION }}*.zip . echo "Renamed files" ls infernet_* - + + - name: Extract release notes from NEWS.md + shell: bash + run: | + # Take the NEWS.md section for this version, else the topmost section. + awk -v ver="${{ env.PACKAGE_NAME }} ${{ env.PACKAGE_VERSION }}" ' + /^# / { if (found) exit; if (substr($0, 3) == ver) { found = 1; next } } + found { print } + ' NEWS.md > RELEASE_NOTES.md + if [ ! -s RELEASE_NOTES.md ]; then + awk 'NR > 1 && /^# / { exit } NR > 1 { print }' NEWS.md > RELEASE_NOTES.md + fi + cat RELEASE_NOTES.md + - name: Create Release and Upload Assets id: create_release uses: softprops/action-gh-release@v2 with: tag_name: ${{ steps.newtag.outputs.tag }} name: Release ${{ steps.newtag.outputs.tag }} + body_path: RELEASE_NOTES.md draft: false prerelease: false fail_on_unmatched_files: true @@ -120,7 +136,7 @@ jobs: infernet_*.tar.gz infernet_*.zip env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} pkgdown: name: Build and deploy website @@ -130,11 +146,11 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v7 - uses: r-lib/actions/setup-r@v2 - - uses: r-lib/actions/setup-pandoc@v1 + - uses: r-lib/actions/setup-pandoc@v2 - uses: r-lib/actions/setup-r-dependencies@v2 with: @@ -143,7 +159,7 @@ jobs: any::rcmdcheck any::pkgdown any::rsconnect - needs: check + needs: website - name: Install package run: R CMD INSTALL . diff --git a/.gitignore b/.gitignore index 02df9f1..580f8ab 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,4 @@ CRAN-SUBMISSION toadd/* .data.csv cache/* +.positai diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..d4bde28 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,29 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## Project overview + +`infernet` is an R package (part of the [stocnet](https://github.com/stocnet) ecosystem) providing the *inferential layer* for network analysis: +conditional uniform graph (CUG) and quadratic assignment procedure (QAP) tests of network statistics, and multiple regression QAP (MRQAP). +It combines a formula front end with a range of estimator range, missing-data handling, and treatments for multimodal, multilevel, and cognitive social structure networks. +It builds on `{manynet}` (network classes and coercion) and `{netrics}` (measures), +so its functions accept matrices, edgelists, `{igraph}`, `{network}`, `{tidygraph}`, or `stocnet` objects, +and one-mode or two-mode networks alike. +Plotting results belongs in `{autograph}`. + +Full package documentation — common dev commands, file organization, the `net_regression()` pipeline, +function body and parallelism conventions, console messaging, dependency practice, test conventions, `NEWS.md` conventions, and branching/CI — lives in [.github/CONTRIBUTING.md](.github/CONTRIBUTING.md). +Read it before adding or restructuring functions, +or when you need the exact `devtools`/testing commands for this repo. + +Note that `{migraph}` still carries older copies of `net_regression()` and the `test_*()` family, +which `{infernet}` supersedes. +A fix made in `{migraph}` needs porting here; a fix made here does not need porting back. + +## Where to make changes + +- Make all changes on the `develop` branch. `develop` is the working branch; `main` is the release branch. +- Do not commit to `main` directly. Do not open a new feature branch for a fix unless the user asks. +- Keep changes reviewable. Leave them as an uncommitted working-tree diff on `develop`, or as local commits on `develop`. Commit only when the user asks. +- Do not push to `origin/develop` without explicit approval. `develop` is shared. diff --git a/DESCRIPTION b/DESCRIPTION index 5ecb80e..d9a01f2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,25 @@ Package: infernet -Title: Basic inferential models for many different types of networks -Version: 0.1.0 -Date: 2025-06-23 -Description: TBD. +Title: Inferential Models for Many Different Types of Networks +Version: 0.1.1 +Description: A set of tools for testing networks. + It includes functions for univariate and multivariate + conditional uniform graph and quadratic assignment procedure testing, + and multiple regression quadratic assignment procedure (MRQAP) models + for a range of families, with random and fixed effects, + missing data, and cognitive social structures. + Built on the 'manynet' package, all functions operate with matrices, + edge lists, and 'igraph', 'network', and 'tidygraph' objects, + and on one-mode and two-mode (bipartite) networks. URL: https://stocnet.github.io/infernet/ BugReports: https://github.com/stocnet/infernet/issues License: MIT + file LICENSE Language: en-GB Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.3 Depends: - R (>= 3.6.0) + R (>= 4.1.0), + manynet (>= 2.3.1), + netrics (>= 1.0.1) Authors@R: c(person(given = "James", family = "Hollway", @@ -24,16 +32,11 @@ Authors@R: comment = c(ORCID = "0000-0003-4288-4732")) ) Imports: - cli, - dplyr, furrr, future, future.apply, - manynet, purrr, - reformulas, - stats, - utils + reformulas Suggests: lme4, nnet, @@ -56,3 +59,5 @@ Config/Needs/website: pkgdown Config/testthat/parallel: true Config/testthat/edition: 3 +Config/testthat/start-first: qap_estimators, model_tests +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index 2711a77..b0cf7af 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,11 +6,12 @@ export(net_regression) export(test_configuration) export(test_permutation) export(test_random) -importFrom(cli,cli_div) -importFrom(cli,cli_end) -importFrom(cli,cli_inform) -importFrom(manynet,bind_node_attributes) -importFrom(manynet,generate_configuration) -importFrom(manynet,generate_random) -importFrom(manynet,is_complex) -importFrom(manynet,is_directed) +importFrom(manynet, + bind_node_attributes, + generate_configuration, + generate_random, + is_complex, + is_directed, + snet_info +) +importFrom(netrics,net_by_heterophily) diff --git a/NEWS.md b/NEWS.md index 6d45571..61cc9c9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,74 @@ +# infernet 0.1.1 + +## Package + +- Updated `DESCRIPTION` + - Raised the R minimum to 4.1.0, since the examples use the native pipe + - Pinned the minimum versions of `{manynet}` and `{netrics}` +- Updated CONTRIBUTING to document the architecture and the house conventions +- Added `README.Rmd` and `pkgdown/_pkgdown.yml` for this package +- Improved the loading messages +- Updated the Github Actions workflows + - Added the PR metadata checks for the version bump and the PR title and body + - Release notes are now taken from the matching `NEWS.md` section + - Updated the action versions in `prchecks` and `pushrelease` +- Improved console messaging to use the `snet_*()` wrappers from `{manynet}` + - Informational output is now silent by default, and follows `snet_verbosity` + - Errors name what is missing and what is available + - Added `thisRequires()`, which names the install command for a suggested package + +## Tests + +- Improved `test_permutation()` by dropping two unused computations +- Updated the `tests` documentation to describe `test_configuration()` +- Updated the examples to use the native pipe `|>` + +## Regression + +- Added a test suite for the regression engine, in four files + - `test-qap_estimators.R` compares each estimator's baseline against the + equivalent `lm()`, `glm()`, `MASS`, `pscl`, `lme4`, or `fixest` fit + - `test-qap_shapes.R` counts the dyads that reach the model for a directed, + an undirected, a two-mode, a pooled, and a partly missing network + - `test-qap_reproducibility.R` fixes the seed contract, sequential and parallel + - `test-qap_control.R` covers the control list and the null-hypothesis choice + - `helper-infernet.R` holds the seeded fixtures and `expect_qap_shape()` +- Fixed a two-mode network being read as a square one + - An 18x14 incidence matrix produced 306 dyads rather than 252 + - `RMPerm()` now permutes the rows and the columns of a rectangular matrix + independently, rather than erroring on the shorter side + - `dist()` and `sim()` now read each mode separately, as `ego()` and + `alter()` already did +- Fixed an undirected network contributing each dyad twice, shrinking the standard error +- Fixed `family = "zip"` failing during permutation + - Several estimators returned backticked coefficient names, which double + semi-partialling could not look up + - Coefficient names are now cleaned on the one path every estimator takes +- Fixed the `{fixest}` path reporting two intercepts where none was absorbed +- Fixed crossed sender and receiver intercepts aborting the run + - Residualising falls back to no random intercepts, with a warning, where the + mixed fit is singular +- Fixed `use_gpu = TRUE` aborting where `{torch}` or CUDA is unavailable +- Improved `control` to reject a name it does not take, and offer the nearest +- Improved the permutation loop to hold back a fitter's convergence warnings + - These printed once per draw; the count of failed draws is still reported +- Improved `lower`, `larger`, and `abs` to carry the same row names under both + null hypotheses +- Improved `HC3()` by dropping a `gc()` call that ran once per permutation +- Improved the missing-predictor error to name the sender, receiver, and + network indices that a formula can also use +- Fixed `tertius()` rejecting a quoted summary function + - `tertius(x, "mean")`, the documented spelling, now works alongside + `tertius(x, mean)` +- Added tests for the `tertius()` spellings and for the missing-attribute error +- Fixed `R/model_regression.R`, which a bad merge left unable to parse + - Restored `.default_control()`, `.is_list_of_graphs()`, and the head of + `.prepare_list_of_graphs()` + - Removed `vectorise_list()`, and the copies of `logit_moments()` and + `logit_resid()` that duplicate `R/qap_gmm.R` + # infernet 0.1.0 ## Package - Initialised package - diff --git a/R/infernet-package.R b/R/infernet-package.R new file mode 100644 index 0000000..4376a4c --- /dev/null +++ b/R/infernet-package.R @@ -0,0 +1,34 @@ +# nocov start + +#' @keywords internal +"_PACKAGE" + +## usethis namespace: start +## usethis namespace: end +NULL + +# Checks for a package in Suggests, and aborts with the install command where it +# is missing. Deliberately not a prompt: `utils::askYesNo()` reads from stdin, +# and a permutation run started from a script would stall on it. +#' @keywords internal +#' @noRd +thisRequires <- function(pkgname, why) { + if (!requireNamespace(pkgname, quietly = TRUE)) { + manynet::snet_abort( + c(paste0("The {.pkg ", pkgname, "} package is required ", why, "."), + i = paste0("Install it with {.run install.packages(\"", pkgname, "\")}."))) + } + invisible(TRUE) +} + +# Suppress R CMD check note +# Namespace in Imports field not imported from: PKG +# All declared Imports should be used. +#' @importFrom netrics net_by_heterophily +ignore_unused_imports <- function() { + # This function exists only to reference functions and suppress R CMD check notes about unused imports. + netrics::net_by_heterophily + NULL +} + +# nocov end diff --git a/R/model_regression.R b/R/model_regression.R index 6703225..5b639d9 100644 --- a/R/model_regression.R +++ b/R/model_regression.R @@ -80,7 +80,7 @@ #' \doi{10.1007/s11336-007-9016-1}. #' @examples #' \dontrun{ -#' networkers <- manynet::ison_networkers %>% +#' networkers <- manynet::ison_networkers |> #' manynet::to_subgraph(Discipline == "Sociology") #' model1 <- net_regression( #' weight ~ ego(Citations) + alter(Citations) + sim(Citations), @@ -93,11 +93,7 @@ net_regression <- function(formula, times = 1000, control = list()) { - ctrl <- .default_control() - if (length(control) > 0) { - ctrl[names(control)] <- control - } - ctrl$method <- match.arg(ctrl$method, choices = c("qap", "qapy")) + ctrl <- .resolve_control(control) if (.is_list_of_graphs(.data)) { prepared <- .prepare_list_of_graphs(formula, .data) @@ -163,6 +159,38 @@ net_regression <- function(formula, # ---- default control ------------------------------------------------------- +# Merges the user's list over the defaults. A name that is not a control is +# rejected rather than added silently: a misspelt name would otherwise leave the +# option it was meant to set at its default, with nothing to say so. +#' @keywords internal +#' @noRd +.resolve_control <- function(control = list()) { + ctrl <- .default_control() + if (length(control) == 0) { + ctrl$method <- match.arg(ctrl$method, choices = c("qap", "qapy")) + return(ctrl) + } + if (is.null(names(control)) || any(!nzchar(names(control)))) { + manynet::snet_abort("Every entry of {.arg control} must be named.") + } + unknown <- setdiff(names(control), names(ctrl)) + if (length(unknown) > 0) { + near <- vapply(unknown, function(u) { + d <- utils::adist(u, names(ctrl), ignore.case = TRUE)[1, ] + if (min(d) <= max(2, nchar(u) %/% 3)) names(ctrl)[which.min(d)] else NA_character_ + }, character(1)) + msg <- c("{.arg control} does not take {.val {unknown}}.") + if (any(!is.na(near))) { + msg <- c(msg, i = "Did you mean {.val {unname(near[!is.na(near)])}}?") + } + msg <- c(msg, i = "Available controls: {.val {names(ctrl)}}.") + manynet::snet_abort(msg) + } + ctrl[names(control)] <- control + ctrl$method <- match.arg(ctrl$method, choices = c("qap", "qapy")) + ctrl +} + .default_control <- function() { list( method = c("qap", "qapy"), @@ -220,8 +248,9 @@ net_regression <- function(formula, keep <- setdiff(seq_along(glist), fail_idx) if (length(keep) == 0) { - stop("None of the supplied networks could be converted for the given ", - "formula. Reasons:\n ", paste(fail_reason, collapse = "\n ")) + manynet::snet_abort( + c("None of the supplied networks could be converted for this formula.", + stats::setNames(unique(fail_reason), rep("x", length(unique(fail_reason)))))) } ref_names <- names(ml_list[[keep[1]]]$mydata) @@ -241,13 +270,11 @@ net_regression <- function(formula, } else { paste0("[", dropped, "]") } - warning("Dropping ", length(dropped), - " network(s) missing one or more predictors: ", - paste(dropped_labels, collapse = ", "), - call. = FALSE) + manynet::snet_warn( + "Dropping {length(dropped)} network{?s} missing one or more predictors: {.val {dropped_labels}}.") } if (length(keep) == 0) { - stop("All supplied networks were dropped due to missing predictors.") + manynet::snet_abort("All the supplied networks were dropped for missing predictors.") } kept_mls <- ml_list[keep] @@ -520,6 +547,22 @@ print.net_regression <- function(x, ..., convertToMatrixList <- function(formula, .data, advise = TRUE) { data <- manynet::as_tidygraph(.data) DV <- manynet::as_matrix(data) + # 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. + side_matrices <- function(attrib, DV, twomode, type) { + if (twomode) { + rows <- matrix(attrib[!type], nrow(DV), ncol(DV)) + cols <- matrix(attrib[type], nrow(DV), ncol(DV), byrow = TRUE) + } else { + rows <- matrix(attrib, nrow(DV), ncol(DV)) + cols <- matrix(attrib, nrow(DV), ncol(DV), byrow = TRUE) + } + list(rows = rows, cols = cols) + } + twomode <- manynet::is_twomode(data) + node_type <- if (twomode) manynet::node_attribute(data, "type") else NULL + names_form <- getRHSNames(formula) .check_formula_vars(names_form$IVnames, data) if (advise) specificationAdvice(names_form$IVnames, data) @@ -576,24 +619,28 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { out } else if (IV[[elem]][1] == "dist") { if (is.character(manynet::node_attribute(data, IV[[elem]][2]))) { - stop("Distance undefined for factors.") + manynet::snet_abort( + c("{.fn dist} is undefined for a categorical attribute.", + i = "Try {.fn same} instead.")) } - rows <- matrix(manynet::node_attribute(data, IV[[elem]][2]), - nrow(DV), ncol(DV)) - cols <- matrix(manynet::node_attribute(data, IV[[elem]][2]), - nrow(DV), ncol(DV), byrow = TRUE) + sides <- side_matrices(manynet::node_attribute(data, IV[[elem]][2]), + DV, twomode, node_type) + rows <- sides$rows + cols <- sides$cols out <- abs(rows - cols) out <- list(out) names(out) <- paste(IV[[elem]], collapse = " ") out } else if (IV[[elem]][1] == "sim") { if (is.character(manynet::node_attribute(data, IV[[elem]][2]))) { - stop("Similarity undefined for factors. Try `same()` instead.") + manynet::snet_abort( + c("{.fn sim} is undefined for a categorical attribute.", + i = "Try {.fn same} instead.")) } - rows <- matrix(manynet::node_attribute(data, IV[[elem]][2]), - nrow(DV), ncol(DV)) - cols <- matrix(manynet::node_attribute(data, IV[[elem]][2]), - nrow(DV), ncol(DV), byrow = TRUE) + sides <- side_matrices(manynet::node_attribute(data, IV[[elem]][2]), + DV, twomode, node_type) + rows <- sides$rows + cols <- sides$cols denom <- max(abs(rows - cols), na.rm = TRUE) if (!is.finite(denom) || denom == 0) denom <- 1 out <- abs(1 - abs(rows - cols) / denom) @@ -609,6 +656,9 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { if (is.na(IV[[elem]][3])) { IV[[elem]][3] <- "mean" } + # The deparsed term keeps the quotation marks of `tertius(x, "mean")`, + # so strip them; otherwise only the unquoted spelling is recognised. + IV[[elem]][3] <- gsub('^"|"$', "", IV[[elem]][3]) out <- t(vapply(seq_len(nrow(DV)), function(x) { if (IV[[elem]][3] == "mean") { @@ -616,7 +666,8 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { } else if (IV[[elem]][3] == "sum") { colSums(val[-x, ], na.rm = TRUE) } else { - stop("tertius summary function not recognised") + manynet::snet_abort( + "{.fn tertius} takes {.val mean} or {.val sum}, not {.val {IV[[elem]][3]}}.") } }, FUN.VALUE = numeric(ncol(DV)))) @@ -631,7 +682,7 @@ convertToMatrixList <- function(formula, .data, advise = TRUE) { names(out) <- IV[[elem]][1] out } else { - stop("Predictor '", IV[[elem]][1], "' not found in the network.") + manynet::snet_abort("Predictor {.val {IV[[elem]][1]}} not found in the network.") } } }) @@ -743,7 +794,11 @@ getRHSNames <- function(formula) { .check_formula_vars <- function(IVnames, data) { node_fns <- c("ego", "alter", "same", "dist", "sim", "tertius") node_attrs <- manynet::net_node_attributes(data) - tie_attrs <- manynet::net_tie_attributes(data) + # 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()) missing_node <- character(0) missing_tie <- character(0) @@ -763,22 +818,26 @@ getRHSNames <- function(formula) { } if (length(missing_node) > 0) { - stop("Node attribute(s) not found: ", - paste(shQuote(unique(missing_node)), collapse = ", "), - ".\n Available: ", - paste(shQuote(node_attrs), collapse = ", "), - call. = FALSE) + manynet::snet_abort( + c("Node attribute{?s} {.val {unique(missing_node)}} not found.", + i = "Available node attributes: {.val {node_attrs}}.")) } if (length(missing_tie) > 0) { - stop("Tie attribute / predictor(s) not found: ", - paste(shQuote(unique(missing_tie)), collapse = ", "), - ".\n Available tie attributes: ", - paste(shQuote(tie_attrs), collapse = ", "), - call. = FALSE) + available <- manynet::net_tie_attributes(data) + structurals <- .structural_vars() + manynet::snet_abort( + c("Tie attribute or predictor{?s} {.val {unique(missing_tie)}} not found.", + i = "Available tie attributes: {.val {available}}.", + i = "Sender, receiver, and network indices are also available as {.val {structurals}}.")) } invisible(TRUE) } +# Columns the engine builds for every dyad, rather than reads off the network. +#' @keywords internal +#' @noRd +.structural_vars <- function() c("sv", "rv", "nv", "pv") + #' @keywords internal #' @noRd @@ -813,9 +872,9 @@ specificationAdvice <- function(formula, data) { if (length(suggests) > 1) { suggests <- paste0(suggests, collapse = ", ") } - cat(paste("When testing for homophily,", - "it is recommended to include all more fundamental effects.\n", - "Try adding", suggests, "to the model specification.\n\n")) + manynet::snet_info( + "When testing for homophily, include all the more fundamental effects.", + "Try adding {suggests} to the model specification.") } } } diff --git a/R/model_tests.R b/R/model_tests.R index 9139743..908b110 100644 --- a/R/model_tests.R +++ b/R/model_tests.R @@ -1,19 +1,28 @@ # Tests of network measures #### #' Tests of network measures -#' +#' @name tests #' @description #' These functions conduct tests of any network-level statistic: #' #' - `test_random()` performs a conditional uniform graph (CUG) test #' of a measure against a distribution of measures on random networks #' of the same dimensions. +#' - `test_configuration()` performs a CUG test against a distribution of +#' measures on random networks that preserve the degree sequence +#' of the original network. #' - `test_permutation()` performs a quadratic assignment procedure (QAP) test #' of a measure against a distribution of measures on permutations #' of the original network. #' -#' @name tests #' @inheritParams regression +#' @param strategy If `{furrr}` is installed, +#' then multiple cores can be used to accelerate the function. +#' By default `"sequential"`, +#' but if multiple cores available, +#' then `"multisession"` or `"multicore"` may be useful. +#' Generally this is useful only when `times` > 1000. +#' See [`{furrr}`](https://furrr.futureverse.org) for more. #' @family models #' @param FUN A graph-level statistic function to test. #' @param ... Additional arguments to be passed on to FUN, @@ -23,17 +32,17 @@ NULL #' @rdname tests #' @importFrom manynet generate_random bind_node_attributes is_directed is_complex #' @examples -#' marvel_friends <- to_unsigned(ison_marvel_relationships) -#' marvel_friends <- to_giant(marvel_friends) %>% +#' marvel_friends <- fict_marvel |> to_uniplex("relationship") |> +#' to_unsigned() |> to_giant() |> #' to_subgraph(PowerOrigin == "Human") -#' (cugtest <- test_random(marvel_friends, manynet::net_heterophily, attribute = "Attractive", +#' (cugtest <- test_random(marvel_friends, net_by_heterophily, attribute = "Attractive", #' times = 200)) #' # plot(cugtest) #' @export test_random <- function(.data, FUN, ..., times = 1000, - strategy = "sequential", - verbose = FALSE){ + strategy = "sequential"){ + verbose <- ifelse(is.null(getOption("snet_verbosity")), FALSE, getOption("snet_verbosity") == "verbose") args <- unlist(list(...)) if (!is.null(args)) { obsd <- FUN(.data, args) @@ -44,12 +53,12 @@ test_random <- function(.data, FUN, ..., on.exit(future::plan(oplan), add = TRUE) rands <- furrr::future_map(1:times, manynet::generate_random, n = .data, .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) if (length(args) > 0) { rands <- furrr::future_map(rands, manynet::bind_node_attributes, object2 = .data, .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) } if (!is.null(args)) { simd <- furrr::future_map_dbl(rands, @@ -76,8 +85,8 @@ test_random <- function(.data, FUN, ..., #' @export test_configuration <- function(.data, FUN, ..., times = 1000, - strategy = "sequential", - verbose = FALSE){ + strategy = "sequential"){ + verbose <- ifelse(is.null(getOption("snet_verbosity")), FALSE, getOption("snet_verbosity") == "verbose") args <- unlist(list(...)) if (!is.null(args)) { obsd <- FUN(.data, args) @@ -86,14 +95,15 @@ test_configuration <- function(.data, FUN, ..., } oplan <- future::plan(strategy) on.exit(future::plan(oplan), add = TRUE) - rands <- furrr::future_map(1:times, manynet::generate_configuration, n = .data, + rands <- furrr::future_map(1:times, + ~ manynet::generate_configuration(.data), .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) if (length(args) > 0) { rands <- furrr::future_map(rands, manynet::bind_node_attributes, object2 = .data, .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) } if (!is.null(args)) { simd <- furrr::future_map_dbl(rands, @@ -118,38 +128,36 @@ test_configuration <- function(.data, FUN, ..., #' @rdname tests #' @examples #' # (qaptest <- test_permutation(marvel_friends, -#' # manynet::net_heterophily, attribute = "Attractive", +#' # net_by_heterophily, attribute = "Attractive", #' # times = 200)) #' # plot(qaptest) #' @export test_permutation <- function(.data, FUN, ..., times = 1000, - strategy = "sequential", - verbose = FALSE){ + strategy = "sequential"){ + verbose <- ifelse(is.null(getOption("snet_verbosity")), FALSE, getOption("snet_verbosity") == "verbose") args <- unlist(list(...)) if (!is.null(args)) { obsd <- FUN(.data, args) } else { obsd <- FUN(.data) } - n <- manynet::net_dims(.data) - d <- manynet::net_density(.data) oplan <- future::plan(strategy) on.exit(future::plan(oplan), add = TRUE) rands <- furrr::future_map(1:times, function(x) manynet::to_permuted(.data), .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) if (!is.null(args)) { simd <- furrr::future_map_dbl(rands, FUN, args, .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) } else { simd <- furrr::future_map_dbl(rands, FUN, .progress = verbose, - .options = furrr::furrr_options(seed = T)) + .options = furrr::furrr_options(seed = TRUE)) } out <- list(test = "QAP", testval = obsd, @@ -171,5 +179,6 @@ print.network_test <- function(x, ..., cat("Observed Value:", x$testval, "\n") cat("Pr(X>=Obs):", x$pgteobs, "\n") cat("Pr(X<=Obs):", x$plteobs, "\n\n") + invisible(x) } diff --git a/R/qap_css.R b/R/qap_css.R index ef4e19c..7b14fe4 100644 --- a/R/qap_css.R +++ b/R/qap_css.R @@ -192,14 +192,20 @@ QAPcssPermEst <- function(i, } if (trial >= max_trials) { - stop("Cannot find valid permutation after ", max_trials, " 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 if (is.null(comp.)) { + # A fit inside the permutation loop runs `reps` 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( - fit_qap_model(mod = mod., + suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., estimator = estimator., @@ -208,7 +214,7 @@ QAPcssPermEst <- function(i, use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random., - reference = reference.), + reference = reference.)), error = function(e) NULL ) if (is.null(perm_fit)) return(NULL) @@ -224,8 +230,12 @@ 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 + # 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( - fit_qap_model(mod = mod., + suppressWarnings(fit_qap_model(mod = mod., pred = predK, family = family., estimator = estimator., @@ -234,7 +244,7 @@ QAPcssPermEst <- function(i, use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random., - reference = reference.), + reference = reference.)), error = function(e) NULL ) if (is.null(perm_fit)) return(NULL) @@ -339,12 +349,13 @@ QAPcss <- function(formula, if (!large) { y <- data[[dep]] if (length(dim(y)) != 3) - stop("data[['", dep, "']] must be a 3-dimensional array ", - "[sender, receiver, perceiver].") + 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) - stop("data[['", dep, "']][[", i, "]] must be a 3D array.") + manynet::snet_abort( + "Network {i} of the dependent variable {.val {dep}} must be a 3-dimensional array.") } } @@ -359,29 +370,36 @@ QAPcss <- function(formula, has_random <- grepl("\\(", mod_str) || parsed$has_random use_fixest <- parsed$use_fixest if (has_random && use_fixest) { - warning("Cannot combine fixest FE and lme4 random effects. ", - "Using lme4 only.") + 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 (has_random && family == "multinom") { - warning("Random intercepts not implemented for multinomial. ", - "Using standard nnet::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") { - warning("Robust SEs not implemented for multinomial.") + manynet::snet_warn( + "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)) { - warning("Undirected mode: sender/receiver random intercepts set to FALSE.") + 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) warning("Results may not be valid when diagonal is used.") + 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)") @@ -393,7 +411,8 @@ QAPcss <- function(formula, n <- dim(data[[dep]])[1] if (!is.null(groups)) { if (length(groups) != n) - stop("groups length (", length(groups), ") != N (", 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)) @@ -484,8 +503,9 @@ QAPcss <- function(formula, for (xi in main) { test_val <- data[[xi]] if (!is.numeric(test_val)) { - warning("Cannot residualise non-numeric predictor '", xi, - "'. Skipping qapspp for this variable.") + 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, @@ -597,8 +617,9 @@ QAPcss <- function(formula, for (xi in main) { test_val <- if (!large) data[[xi]] else data[[xi]][[1]] if (!is.numeric(test_val)) { - warning("Cannot residualise non-numeric predictor '", xi, - "'. Skipping qapspp for this variable.") + manynet::snet_warn( + c("Cannot residualise the non-numeric predictor {.val {xi}}.", + i = "Skipping double semi-partialling for this predictor.")) next } diff --git a/R/qap_engine.R b/R/qap_engine.R index f582ed9..8290bf9 100644 --- a/R/qap_engine.R +++ b/R/qap_engine.R @@ -50,8 +50,9 @@ QAPglm <- function(formula, has_random <- grepl("\\(", mod_str) || parsed$has_random use_fixest <- parsed$use_fixest if (has_random && use_fixest) { - warning("Cannot combine fixest FE and lme4 random effects. ", - "Using lme4 random effects only.") + manynet::snet_warn( + c("Cannot combine {.pkg fixest} fixed effects with {.pkg lme4} random effects.", + i = "Using the random effects only.")) use_fixest <- FALSE } @@ -129,8 +130,19 @@ QAPglm <- function(formula, if ((nullhyp == "qapspp") && (length(main) == 1)) nullhyp <- "qapy" - if (use_gpu && family == "gaussian" && !has_random && !use_fixest && - is.null(comparison) && !large) { + # 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 + } + + if (use_gpu) { if (nullhyp == "qapy") { gpu_res <- gpu_batch_ols(data = data, @@ -147,10 +159,10 @@ QAPglm <- function(formula, } else if (nullhyp == "qapspp") { n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) + 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 - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) for (xi in main) { xR <- residualise_predictor(xi, pred, main, @@ -228,10 +240,10 @@ QAPglm <- function(formula, } else if (nullhyp == "qapspp") { if (is.null(comparison)) { n_coefs <- length(fit$base$coefficients) - fit$lower <- matrix(NA, nrow = 2, ncol = n_coefs) + 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 - colnames(fit$lower) <- colnames(fit$larger) <- - colnames(fit$abs) <- names(fit$base$coefficients) } else { fit$lower <- fit$larger <- fit$abs <- vector("list", length(comparison)) @@ -239,10 +251,10 @@ QAPglm <- function(formula, 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$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]] - colnames(fit$lower[[k]]) <- colnames(fit$larger[[k]]) <- - colnames(fit$abs[[k]]) <- names(fit$base[[k]]$coefficients) } } @@ -422,8 +434,12 @@ 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 + # 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( - fit_qap_model(mod = mod., + suppressWarnings(fit_qap_model(mod = mod., pred = pred, family = family., estimator = estimator., @@ -432,7 +448,7 @@ QAPglmPermEst <- function(i, use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random., - reference = reference.), + reference = reference.)), error = function(e) NULL ) if (is.null(perm_fit)) return(NULL) @@ -448,8 +464,12 @@ 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 + # 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( - fit_qap_model(mod = mod., + suppressWarnings(fit_qap_model(mod = mod., pred = predK, family = family., estimator = estimator., @@ -458,7 +478,7 @@ QAPglmPermEst <- function(i, use_robust_errors = use_robust_errors., main_vars = main_vars., has_random = has_random., - reference = reference.), + reference = reference.)), error = function(e) NULL ) if (is.null(perm_fit)) return(NULL) diff --git a/R/qap_gpu.R b/R/qap_gpu.R index 8ad1ad5..92cdf97 100644 --- a/R/qap_gpu.R +++ b/R/qap_gpu.R @@ -9,13 +9,10 @@ gpu_batch_ols <- function(data, parsed, mode, diag, groups, reps, baseline_fit, perm_var = NULL, batch_size = 500, device = "cuda") { - if (!requireNamespace("torch", quietly = TRUE)) { - stop("The 'torch' package is required for GPU acceleration. ", - "Install it with: install.packages('torch')") - } + thisRequires("torch", "for GPU acceleration") if (device == "cuda" && !torch::cuda_is_available()) { - message("CUDA not available. Falling back to CPU torch.") + manynet::snet_info("CUDA is not available, so falling back to CPU {.pkg torch}.") device <- "cpu" } @@ -167,12 +164,10 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, baseline_fit, perm_var = NULL, batch_size = 500, device = "cuda") { - if (!requireNamespace("torch", quietly = TRUE)) { - stop("The 'torch' package is required for GPU acceleration.") - } + thisRequires("torch", "for GPU acceleration") if (device == "cuda" && !torch::cuda_is_available()) { - message("CUDA not available. Falling back to CPU torch.") + manynet::snet_info("CUDA is not available, so falling back to CPU {.pkg torch}.") device <- "cpu" } @@ -314,5 +309,10 @@ gpu_batch_ols_css <- function(data, parsed, mode, diag, groups, reps, #' @noRd gpu_available <- function() { if (!requireNamespace("torch", quietly = TRUE)) return(FALSE) - torch::cuda_is_available() + # {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/R/qap_misc.R b/R/qap_misc.R index 9cb9823..d362a9b 100644 --- a/R/qap_misc.R +++ b/R/qap_misc.R @@ -103,8 +103,9 @@ df_to_mat <- function(df, if (loops) n_s * n_r else n_s * n_r - min(n_s, n_r) } if (anyNA(df[var_names]) || nrow(df) != expected) { - warning("Incomplete dyadic data: some cells will be NA.", - "\nCheck the data and consider coding matrices manually.") + manynet::snet_warn( + c("Incomplete dyadic data, so some cells will be {.val NA}.", + i = "Check the data, or code the matrices manually.")) } make_structure <- function(var) { diff --git a/R/qap_utils.R b/R/qap_utils.R index 7a02a39..2afd0c5 100644 --- a/R/qap_utils.R +++ b/R/qap_utils.R @@ -78,13 +78,13 @@ build_internal_formula <- function(formula, validate_qap_input <- function(data, parsed, css = FALSE) { dep <- parsed$dependent if (!(dep %in% names(data))) { - stop("Dependent variable '", dep, "' not found in data.") + 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))) { - stop("Predictor '", v, "' not found in data.") + manynet::snet_abort("Predictor {.val {v}} not found in the data.") } } @@ -93,21 +93,25 @@ validate_qap_input <- function(data, parsed, css = FALSE) { if (!css) { if (!large) { - if (!is.matrix(y)) stop("data[['", dep, "']] must be a matrix.") + if (!is.matrix(y)) + manynet::snet_abort("The dependent variable {.val {dep}} must be a matrix.") } else { for (i in seq_along(y)) { if (!is.matrix(y[[i]])) - stop("data[['", dep, "']][[", i, "]] must be a matrix.") + manynet::snet_abort( + "Network {i} of the dependent variable {.val {dep}} must be a matrix.") } } } else { if (!large) { if (length(dim(y)) != 3) - stop("data[['", dep, "']] must be a 3-dimensional array.") + manynet::snet_abort( + "The dependent variable {.val {dep}} must be a 3-dimensional array.") } else { for (i in seq_along(y)) { if (length(dim(y[[i]])) != 3) - stop("data[['", dep, "']][[", i, "]] must be a 3D array.") + manynet::snet_abort( + "Network {i} of the dependent variable {.val {dep}} must be a 3-dimensional array.") } } } @@ -147,6 +151,18 @@ run_permutations <- function(reps, FUN, ...) { # ---- matrix permutation ----------------------------------------------------- +# One permutation of the node order, respecting a blocking factor where its +# length matches the mode being permuted. A `groups` vector of the wrong length +# cannot describe this mode, so that mode permutes freely rather than silently +# recycling the factor, which is what `split()` did before. +#' @keywords internal +#' @noRd +.perm_order <- function(n, groups = NULL) { + if (is.null(groups) || length(groups) != n) return(sample(seq_len(n))) + groups <- as.character(groups) + unsplit(lapply(split(seq_len(n), groups), FUN = sample), groups) +} + #' @keywords internal #' @noRd RMPerm <- function(m, groups = NULL, CSS = FALSE) { @@ -162,8 +178,20 @@ RMPerm <- function(m, groups = NULL, CSS = FALSE) { } if (length(dim(m)) == 2) { - o <- unsplit(lapply(split(1:dim(m)[1], groups), FUN = sample), groups) - p <- matrix(data = m[o, o], nrow = dim(m)[1], ncol = dim(m)[2]) + nr <- dim(m)[1] + nc <- dim(m)[2] + if (nr == nc) { + o <- unsplit(lapply(split(1:nr, groups), FUN = sample), groups) + p <- matrix(data = m[o, o], nrow = nr, ncol = nc) + } else { + # A two-mode incidence matrix has two nodesets, so the rows and the + # columns permute independently. Permuting both by one order, as the + # square case does, indexes past the shorter side and errors. + or <- .perm_order(nr, groups) + oc <- .perm_order(nc, groups) + p <- matrix(data = m[or, oc], nrow = nr, ncol = nc) + } + dimnames(p) <- dimnames(m) } else if (CSS) { p <- array(dim = c(dim(m)[1], dim(m)[2], dim(m)[3])) o <- unsplit(lapply(split(1:dim(m)[2], groups), FUN = sample), groups) @@ -193,9 +221,21 @@ make_qap_data <- function(y, x, g = NULL, diag = FALSE, mode = "digraph", x[[xi]] <- RMPerm(x[[xi]], g) } - n <- dim(y)[1] - valid <- matrix(TRUE, n, n) - if (!diag) diag(valid) <- FALSE + # The dependent matrix is square for a one-mode network and rectangular for a + # two-mode one. Reading `nc` from the matrix rather than assuming `nr` is what + # keeps a two-mode network from being read as a square one, which used to + # wrap past the last column and invent dyads. + nr <- dim(y)[1] + nc <- dim(y)[2] + square <- identical(nr, nc) + + valid <- matrix(TRUE, nr, nc) + if (!diag && square) diag(valid) <- FALSE + # An undirected one-mode network holds each dyad twice, once on each side of + # 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 for (var in seq_len(nx)) { valid[is.na(x[[var]])] <- FALSE @@ -210,23 +250,22 @@ make_qap_data <- function(y, x, g = NULL, diag = FALSE, mode = "digraph", } if (sum(vv) == 0) { - stop("No valid dyads remain after removing NA and diagonal cells for ", - "network ", net, ". Check that your predictors and outcome have ", - "non-missing values for overlapping node pairs.", - call. = FALSE) + manynet::snet_abort( + c("No valid dyads remain in network {net} after dropping NA and diagonal cells.", + i = "Check that the predictors and the outcome are observed for the same node pairs.")) } pred <- data.frame( - location = as.vector(matrix(seq_len(n^2), n, n))[vv], + location = as.vector(matrix(seq_len(nr * nc), nr, nc))[vv], yv = as.vector(y)[vv] ) pred$nv <- as.factor(net) - sv <- matrix(seq_len(n), n, n) + sv <- matrix(seq_len(nr), nr, nc) sv[!valid] <- NA pred$sv <- as.vector(sv)[vv] - rv <- t(matrix(seq_len(n), n, n)) + rv <- matrix(seq_len(nc), nr, nc, byrow = TRUE) rv[!valid] <- NA pred$rv <- as.vector(rv)[vv] @@ -252,17 +291,34 @@ HC3 <- function(X, e) { } h <- apply(XO, 1, hf, XTXINV = XTXINV) om <- e^2 / (1 - h)^2 - x <- sqrt(diag(t(t(XTXINV %*% t(XO)) * om) %*% XO %*% XTXINV)) - gc() - return(x) + # No gc() here: this runs once per permutation, and forcing a collection + # thousands of times costs far more than the memory it returns. + sqrt(diag(t(t(XTXINV %*% t(XO)) * om) %*% XO %*% XTXINV)) } # ---- baseline + perm fitting ------------------------------------------------ +# 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, +# so the names are cleaned here, where every path passes through exactly once. +#' @keywords internal +#' @noRd +fit_qap_model <- function(...) { + fit <- .fit_qap_model(...) + for (el in c("coefficients", "t", "zi_coefficients")) { + if (!is.null(fit[[el]]) && !is.null(names(fit[[el]]))) { + names(fit[[el]]) <- gsub("`", "", names(fit[[el]]), fixed = TRUE) + } + } + fit +} + #' @keywords internal #' @noRd -fit_qap_model <- function(mod, pred, family, +.fit_qap_model <- function(mod, pred, family, estimator = "standard", use_fixest = FALSE, fixest_se_cluster = NULL, @@ -279,8 +335,7 @@ fit_qap_model <- function(mod, pred, family, if (!is.null(reference)) { pred[[dep_var]] <- stats::relevel(pred[[dep_var]], ref = reference) } - if (!requireNamespace("nnet", quietly = TRUE)) - stop("Package 'nnet' is required for multinomial models.") + 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) / @@ -290,8 +345,7 @@ fit_qap_model <- function(mod, pred, family, } if (estimator == "gmm") { - if (!requireNamespace("gmm", quietly = TRUE)) - stop("Package 'gmm' is required for GMM estimation.") + thisRequires("gmm", "for GMM estimation") y_vec <- pred[[dep_var]] x_mat <- cbind(1, as.matrix(pred[, main_vars, drop = FALSE])) @@ -326,8 +380,9 @@ fit_qap_model <- function(mod, pred, family, resid <- zip_resid(base_model) has_extra_param <- TRUE } else { - stop("GMM estimator is available for binomial, poisson, negbin, ", - "and zip families.") + 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 @@ -360,8 +415,7 @@ fit_qap_model <- function(mod, pred, family, if (family == "zip" && estimator == "standard") { if (has_random) { - if (!requireNamespace("glmmTMB", quietly = TRUE)) - stop("Package 'glmmTMB' is required for mixed ZIP models.") + thisRequires("glmmTMB", "for mixed zero-inflated Poisson models") base_model <- glmmTMB::glmmTMB(mod, data = pred, family = stats::poisson(), ziformula = ~1) @@ -376,8 +430,7 @@ fit_qap_model <- function(mod, pred, family, fit$random.intercepts[[rV]] <- re[[rV]][, 1] } } else { - if (!requireNamespace("pscl", quietly = TRUE)) - stop("Package 'pscl' is required for zero-inflated Poisson models.") + thisRequires("pscl", "for zero-inflated Poisson models") base_model <- pscl::zeroinfl(mod, data = pred, dist = "poisson") fit$coefficients <- base_model$coefficients$count resid <- stats::residuals(base_model, type = "response") @@ -395,23 +448,38 @@ fit_qap_model <- function(mod, pred, family, if (!has_random) { if (use_fixest) { - if (!requireNamespace("fixest", quietly = TRUE)) - stop("Package 'fixest' is required when fixest_se_cluster or fixed effects with | ") + 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) - fit$coefficients <- c("(Intercept)" = NA, base_model$coefficients) + # {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 <- fit$coefficients / c(NA, hc[-1]) + fit$t <- if (absorbed) { + fit$coefficients / c(NA, hc[-1]) + } else { + fit$coefficients / hc + } } else { fe_se <- sqrt(diag(stats::vcov(base_model))) - fit$t <- c("(Intercept)" = NA, - base_model$coefficients / fe_se) + fe_t <- fe_coefs / fe_se + fit$t <- if (absorbed) c("(Intercept)" = NA, fe_t) else fe_t } names(fit$t) <- names(fit$coefficients) @@ -428,8 +496,7 @@ fit_qap_model <- function(mod, pred, family, fit$r.squared <- summary(base_model)$r.squared fit$adj.r.squared <- summary(base_model)$adj.r.squared } else if (family == "negbin") { - if (!requireNamespace("MASS", quietly = TRUE)) - stop("Package 'MASS' is required for negative binomial models.") + thisRequires("MASS", "for negative binomial models") base_model <- MASS::glm.nb(mod, data = pred) fit$theta <- base_model$theta } else { @@ -447,12 +514,10 @@ fit_qap_model <- function(mod, pred, family, } } else { if (family == "gaussian") { - if (!requireNamespace("lme4", quietly = TRUE)) - stop("Package 'lme4' is required for random effects.") + thisRequires("lme4", "for random effects") base_model <- lme4::lmer(mod, data = pred) } else if (family == "negbin") { - if (!requireNamespace("glmmTMB", quietly = TRUE)) - stop("Package 'glmmTMB' is required for mixed negative binomial models.") + thisRequires("glmmTMB", "for mixed negative binomial models") base_model <- glmmTMB::glmmTMB(mod, data = pred, family = glmmTMB::nbinom2()) fit$coefficients <- glmmTMB::fixef(base_model)$cond @@ -473,8 +538,7 @@ fit_qap_model <- function(mod, pred, family, fit$base_model <- base_model return(fit) } else { - if (!requireNamespace("lme4", quietly = TRUE)) - stop("Package 'lme4' is required for random effects.") + thisRequires("lme4", "for random effects") base_model <- lme4::glmer(mod, data = pred, family = family, control = lme4::glmerControl( calc.derivs = FALSE, @@ -500,14 +564,6 @@ fit_qap_model <- function(mod, pred, family, } fit$base_model <- base_model - - if (!is.null(fit$coefficients) && !is.null(names(fit$coefficients))) { - names(fit$coefficients) <- gsub("`", "", names(fit$coefficients), fixed = TRUE) - } - if (!is.null(fit$t) && !is.null(names(fit$t))) { - names(fit$t) <- gsub("`", "", names(fit$t), fixed = TRUE) - } - return(fit) } @@ -540,10 +596,13 @@ compare_perm_to_baseline <- function(perm_coefs, perm_t, base_fit, aggregate_perm_results <- function(results, reps) { results <- Filter(Negate(is.null), results) n_valid <- length(results) - if (n_valid == 0) stop("All permutations failed to converge.") + if (n_valid == 0) + manynet::snet_abort( + c("All {reps} permutations failed to converge.", + i = "Try a simpler model, another {.arg family}, or fewer predictors.")) if (n_valid < reps) { - warning(reps - n_valid, " of ", reps, - " permutations failed and were excluded.") + manynet::snet_warn( + "{reps - n_valid} of {reps} permutation{?s} failed and {?was/were} excluded.") } resL <- unlist(results, recursive = FALSE) list( @@ -578,9 +637,24 @@ residualise_predictor <- function(xi, pred, main_vars, if (!has_random) { xm <- stats::lm(modx, data = pred) } else { - if (!requireNamespace("lme4", quietly = TRUE)) - stop("Package 'lme4' is required for random effects.") - xm <- lme4::lmer(modx, data = pred) + thisRequires("lme4", "for random effects") + # Residualising is a step towards the null distribution, not a result the + # user reads, so a degenerate mixed fit here must not abort the whole run. + # Crossed sender and receiver intercepts on a predictor are often singular, + # and `lmer()` then stops with "Downdated VtV is not positive definite". + xm <- tryCatch(suppressWarnings(lme4::lmer(modx, data = pred)), + error = function(e) NULL) + if (is.null(xm)) { + manynet::snet_warn( + c("Could not residualise {.val {xi}} with random intercepts.", + i = "Residualising it without them instead.")) + xm <- stats::lm(stats::as.formula( + paste(bq(xi), "~ 1", + if (length(others) > 0) + paste("+", paste(vapply(others, bq, character(1)), + collapse = " + ")) else "")), + data = pred) + } } stats::residuals(xm) } diff --git a/R/zzz.R b/R/zzz.R index d6b4dbb..944a999 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,59 +1,41 @@ -#' @importFrom cli cli_div cli_inform cli_end +# nocov start + +# The stocnet version check lives in {migraph}, which loads and checks the whole +# stack at once. This package is one of the packages it checks, so it does no +# version check of its own. +#' @importFrom manynet snet_info .onAttach <- function(...) { - # suppressMessages(suppressPackageStartupMessages(library("manynet", warn.conflicts = FALSE))) if (!interactive()) return() - - # options(manynet_verbosity = getOption("manynet_verbosity", "verbose")) + options(stocnet_theme = getOption("stocnet_theme", "default")) - # options(cli.theme = manynet_console_theme()) - # options(cli.progress_clear = TRUE) - - # pkgs <- as.data.frame(utils::available.packages(utils::contrib.url(getOption("repos")))) - # - # cran_version <- pkgs[pkgs$Package == "manynet","Version"] local_version <- utils::packageVersion("infernet") - cli::cli_inform("You are using {.pkg infernet} version {.version {local_version}}.", - class = "packageStartupMessage") - old.list <- as.data.frame(utils::old.packages()) - behind_cran <- "infernet" %in% old.list$Package - + manynet::snet_info("You are using {.infr infernet} version {.version {local_version}}.") + greet_startup_cli <- function() { tips <- c( - "i" = "There are lots of ways to contribute to {.pkg infernet} at {.url https://github.com/stocnet/infernet/}.", - "i" = "Please let us know any bugs, issues, or feature requests at {.url https://github.com/stocnet/infernet/issues}. It's really helpful!", - # "i" = "To eliminate package startup messages, use: `suppressPackageStartupMessages(library({.pkg autograph}))`.", - # "i" = "Changing the theme of all your graphs is straightforward with `set_manynet_theme()`", - # "i" = "If there are too many messages in the console, run `options(manynet_verbosity = 'quiet')`", - "i" = "Visit the website to learn more: {.url https://stocnet.github.io/infernet/}.", - "i" = "We recommend the 'Function Overview' page online to discover new analytic opportunities: {.url https://stocnet.github.io/infernet/reference/index.html}.", - # "i" = "Star me at {.url https://github.com/users/follow?target=jhollway}.", - # "i" = "You can list all the tutorials available in {.pkg manynet} using {.fn run_tute}, and run them too!", - "i" = "Discover all the {.emph stocnet} R packages at {.url https://github.com/stocnet/}." + "i" = "Share bugs, issues, or feature requests at {.url https://github.com/stocnet/infernet/issues}.", + "i" = "If too many messages appear in the console, run {.run base::options(snet_verbosity = 'quiet')}", + "i" = "Explore changes since the last version with {.run [news(package = 'infernet')](utils::news(package = 'infernet'))}.", + # "i" = "Test any network statistic against a null distribution with {.fn test_random}, {.fn test_configuration}, or {.fn test_permutation}.", + # "i" = "Regress a network on nodal and dyadic covariates with {.fn net_regression}.", + # "i" = "Write {.code ego()}, {.code alter()}, {.code same()}, {.code dist()}, {.code sim()}, or {.code tertius()} in a formula to build a predictor from a nodal attribute.", + "i" = "Speed up a long run with {.code control = list(strategy = 'multisession')}.", + # "i" = "Measures to test are in {.tric netrics}; plots of results are in {.auto autograph}.", + "i" = "Visit {.url https://stocnet.github.io/infernet/} to learn more.", + "i" = "Discover new functions at {.url https://stocnet.github.io/infernet/reference/index.html}.", + "i" = "Discover {.emph stocnet} R packages at {.url https://github.com/stocnet/}." ) - cli::cli_inform(sample(tips, 1), class = "packageStartupMessage") + manynet::snet_info(sample(tips, 1)) } - if (interactive()) { - if (behind_cran) { - msg <- "A new version of infernet is available with bug fixes and new features." - packageStartupMessage(msg, "\nWould you like to install it?") - if (utils::menu(c("Yes", "No")) == 1) { - utils::update.packages("infernet") - } - } else { - greet_startup_cli() - # packageStartupMessage(paste(strwrap(tip), collapse = "\n")) - } - } + greet_startup_cli() } +# nocov end + # Global variables #### # defining global variables more centrally utils::globalVariables(c(".data")) - - - - diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..83fb201 --- /dev/null +++ b/README.Rmd @@ -0,0 +1,103 @@ +--- +output: github_document +--- + + + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "man/figures/README-", + out.width = "100%" +) +``` + +# infernet +infernet logo + + +[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) +![CRAN/METACRAN](https://img.shields.io/cran/v/infernet) +![GitHub release (latest by date)](https://img.shields.io/github/v/release/stocnet/infernet) +![GitHub Release Date](https://img.shields.io/github/release-date/stocnet/infernet) +[![Codecov test coverage](https://codecov.io/gh/stocnet/infernet/branch/main/graph/badge.svg)](https://app.codecov.io/gh/stocnet/infernet?branch=main) + + +## About the package + +`{infernet}` is the inferential layer of the [stocnet](https://github.com/stocnet) +ecosystem. It offers two things: + +- **Tests of network statistics.** `test_random()` runs a conditional uniform +graph (CUG) test, `test_configuration()` conditions on the degree sequence, +and `test_permutation()` runs a quadratic assignment procedure (QAP) test. +Each takes any graph-level statistic and compares it against a simulated +null distribution. +- **Network regression.** `net_regression()` fits a multiple regression +quadratic assignment procedure (MRQAP) model, using either Dekker et al's +double semi-partialling or a permutation of the dependent network alone. + +The package combines the formula interface of `{migraph}` with the estimators, +missing-data handling, and cognitive social structure treatment of +Robert Krause's `MrQAP`. + +## A formula you can reuse + +Models are specified with a formula, so the same specification can be moved +between networks and between model types: + +```{r example, eval = FALSE} +library(infernet) + +networkers <- manynet::ison_networkers |> + manynet::to_subgraph(Discipline == "Sociology") + +net_regression(weight ~ ego(Citations) + alter(Citations) + sim(Citations), + networkers, times = 200) +``` + +Alongside plain references to other networks, which enter as dyadic +covariates, the right-hand side accepts: + +| Term | Constructs a matrix of | +|---|---| +| `ego(attr)` | the sender's value of a nodal attribute | +| `alter(attr)` | the receiver's value | +| `same(attr)` | 1 where sender and receiver share an attribute value | +| `dist(attr)` | the absolute difference in a numeric attribute | +| `sim(attr)` | the proportional similarity in a numeric attribute | +| `tertius(attr, fn)` | an aggregate of an attribute over a node's other ties | + +Further options are passed through a `control` list: the model family, +the null hypothesis, random or fixed effects, robust standard errors, +the parallel strategy, and an optional `{torch}` path for running +permutations on a GPU. + +## Installation + +### Development + +`{infernet}` is not yet on CRAN. +The latest binary releases for all major OSes -- Windows, Mac, and Linux -- +can be found [here](https://github.com/stocnet/infernet/releases/latest). +Download the appropriate binary for your operating system, +and install using an adapted version of the following commands: + +- For Windows: `install.packages("~/Downloads/infernet_winOS.zip", repos = NULL)` +- For Mac: `install.packages("~/Downloads/infernet_macOS.tgz", repos = NULL)` +- For Unix: `install.packages("~/Downloads/infernet_linuxOS.tar.gz", repos = NULL)` + +To install from source, +please install the `{remotes}` package from CRAN and then: + +- For latest stable version: +`remotes::install_github("stocnet/infernet")` +- For latest development version: +`remotes::install_github("stocnet/infernet@develop")` + +## Funding details + +Development on this package has been funded by the Swiss National Science Foundation (SNSF) +[Grant Number 188976](https://data.snf.ch/grants/grant/188976): +"Power and Networks and the Rate of Change in Institutional Complexes" (PANARCHIC). diff --git a/README.md b/README.md new file mode 100644 index 0000000..d002ab1 --- /dev/null +++ b/README.md @@ -0,0 +1,105 @@ + + + +# infernet + +infernet logo + + + +[![Lifecycle: +experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) +![CRAN/METACRAN](https://img.shields.io/cran/v/infernet) ![GitHub +release (latest by +date)](https://img.shields.io/github/v/release/stocnet/infernet) +![GitHub Release +Date](https://img.shields.io/github/release-date/stocnet/infernet) +[![Codecov test +coverage](https://codecov.io/gh/stocnet/infernet/branch/main/graph/badge.svg)](https://app.codecov.io/gh/stocnet/infernet?branch=main) + + +## About the package + +`{infernet}` is the inferential layer of the +[stocnet](https://github.com/stocnet) ecosystem. It offers two things: + +- **Tests of network statistics.** `test_random()` runs a conditional + uniform graph (CUG) test, `test_configuration()` conditions on the + degree sequence, and `test_permutation()` runs a quadratic assignment + procedure (QAP) test. Each takes any graph-level statistic and + compares it against a simulated null distribution. +- **Network regression.** `net_regression()` fits a multiple regression + quadratic assignment procedure (MRQAP) model, using either Dekker et + al’s double semi-partialling or a permutation of the dependent network + alone. + +It offers these capabilities for one-mode and two-mode networks, and for unimodal, +directed, undirected, weighted, or cognitive social structure networks alike. +It accepts matrices, edgelists, +`{igraph}`, `{network}`, `{tidygraph}`, or `stocnet` objects, +and handles missing-data gracefully. +It can run permutations in parallel, and can optionally use a GPU via `{torch}`. + +## A formula you can reuse + +Models are specified with a formula, so the same specification can be +moved between networks and between model types: + +``` r +library(infernet) + +networkers <- manynet::ison_networkers |> + manynet::to_subgraph(Discipline == "Sociology") + +net_regression(weight ~ ego(Citations) + alter(Citations) + sim(Citations), + networkers, times = 200) +``` + +Alongside plain references to other networks, which enter as dyadic +covariates, the right-hand side accepts: + +| Term | Constructs a matrix of | +|---------------------|-------------------------------------------------------| +| `ego(attr)` | the sender’s value of a nodal attribute | +| `alter(attr)` | the receiver’s value | +| `same(attr)` | 1 where sender and receiver share an attribute value | +| `dist(attr)` | the absolute difference in a numeric attribute | +| `sim(attr)` | the proportional similarity in a numeric attribute | +| `tertius(attr, fn)` | an aggregate of an attribute over a node’s other ties | + +Further options are passed through a `control` list: the model family, +the null hypothesis, random or fixed effects, robust standard errors, +the parallel strategy, and an optional `{torch}` path for running +permutations on a GPU. + +## Installation + +### Development + +`{infernet}` is not yet on CRAN. The latest binary releases for all +major OSes – Windows, Mac, and Linux – can be found +[here](https://github.com/stocnet/infernet/releases/latest). Download +the appropriate binary for your operating system, and install using an +adapted version of the following commands: + +- For Windows: + `install.packages("~/Downloads/infernet_winOS.zip", repos = NULL)` +- For Mac: + `install.packages("~/Downloads/infernet_macOS.tgz", repos = NULL)` +- For Unix: + `install.packages("~/Downloads/infernet_linuxOS.tar.gz", repos = NULL)` + +To install from source, please install the `{remotes}` package from CRAN +and then: + +- For latest stable version: + `remotes::install_github("stocnet/infernet")` +- For latest development version: + `remotes::install_github("stocnet/infernet@develop")` + +## Funding details + +Development on this package has been funded by the Swiss National +Science Foundation (SNSF) [Grant Number +188976](https://data.snf.ch/grants/grant/188976): “Power and Networks +and the Rate of Change in Institutional Complexes” (PANARCHIC). diff --git a/cran-comments.md b/cran-comments.md new file mode 100644 index 0000000..da2adfd --- /dev/null +++ b/cran-comments.md @@ -0,0 +1,22 @@ +## Submission + +This is a first submission of `infernet` to CRAN. + +## Test environments + +* local R installation, aarch64-apple-darwin20, R 4.6.1 +* macOS (on Github Actions), R release +* Microsoft Windows Server 2022 (on Github Actions), R release +* Ubuntu 24.04 (on Github Actions), R release + +## R CMD check results + +0 errors | 0 warnings | 0 notes + +## Notes for reviewers + +* The package Depends on `manynet` and `netrics`, both on CRAN. +* All estimator-specific packages (`lme4`, `glmmTMB`, `fixest`, `gmm`, `MASS`, + `pscl`, `nnet`, `torch`) are in Suggests, and every code path that needs one + guards with `requireNamespace()`. +* Examples keep the number of permutations low so that they run quickly. diff --git a/inst/infernet.png b/inst/infernet.png new file mode 100644 index 0000000..7997505 Binary files /dev/null and b/inst/infernet.png differ diff --git a/man/figures/logo.png b/man/figures/logo.png new file mode 100644 index 0000000..7997505 Binary files /dev/null and b/man/figures/logo.png differ diff --git a/man/infernet-package.Rd b/man/infernet-package.Rd new file mode 100644 index 0000000..08ec645 --- /dev/null +++ b/man/infernet-package.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/infernet-package.R +\docType{package} +\name{infernet-package} +\alias{infernet} +\alias{infernet-package} +\title{infernet: Inferential Models for Many Different Types of Networks} +\description{ +\if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}} + +A set of tools for testing networks. It includes functions for univariate and multivariate conditional uniform graph and quadratic assignment procedure testing, and multiple regression quadratic assignment procedure (MRQAP) models for a range of families, with random and fixed effects, missing data, and cognitive social structures. Built on the 'manynet' package, all functions operate with matrices, edge lists, and 'igraph', 'network', and 'tidygraph' objects, and on one-mode and two-mode (bipartite) networks. +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://stocnet.github.io/infernet/} + \item Report bugs at \url{https://github.com/stocnet/infernet/issues} +} + +} +\author{ +\strong{Maintainer}: James Hollway \email{james.hollway@graduateinstitute.ch} (\href{https://orcid.org/0000-0002-8361-9647}{ORCID}) (IHEID) [contributor] + +Authors: +\itemize{ + \item James Hollway \email{james.hollway@graduateinstitute.ch} (\href{https://orcid.org/0000-0002-8361-9647}{ORCID}) (IHEID) [contributor] +} + +Other contributors: +\itemize{ + \item Robert Krause (\href{https://orcid.org/0000-0003-4288-4732}{ORCID}) [contributor] +} + +} +\keyword{internal} diff --git a/man/regression.Rd b/man/regression.Rd index 2fe38b2..0b5c0bd 100644 --- a/man/regression.Rd +++ b/man/regression.Rd @@ -27,7 +27,7 @@ dyadic covariate. }} \item{.data}{A manynet-consistent network (see -\code{\link[manynet:coerce_graph]{manynet::as_tidygraph()}}), or a list of such networks. When a list is +\code{\link[manynet:as_tidygraph]{manynet::as_tidygraph()}}), or a list of such networks. When a list is supplied the model is fit jointly; graphs missing any predictor are dropped with a warning.} @@ -96,7 +96,7 @@ dropped with a warning and the remaining networks are pooled. } \examples{ \dontrun{ -networkers <- manynet::ison_networkers \%>\% +networkers <- manynet::ison_networkers |> manynet::to_subgraph(Discipline == "Sociology") model1 <- net_regression( weight ~ ego(Citations) + alter(Citations) + sim(Citations), @@ -118,7 +118,7 @@ conditions." \doi{10.1007/s11336-007-9016-1}. } \seealso{ -Other models: +Other models: \code{\link{tests}} } \concept{models} diff --git a/man/tests.Rd b/man/tests.Rd index 938a09e..57c3f7b 100644 --- a/man/tests.Rd +++ b/man/tests.Rd @@ -7,36 +7,15 @@ \alias{test_permutation} \title{Tests of network measures} \usage{ -test_random( - .data, - FUN, - ..., - times = 1000, - strategy = "sequential", - verbose = FALSE -) +test_random(.data, FUN, ..., times = 1000, strategy = "sequential") -test_configuration( - .data, - FUN, - ..., - times = 1000, - strategy = "sequential", - verbose = FALSE -) +test_configuration(.data, FUN, ..., times = 1000, strategy = "sequential") -test_permutation( - .data, - FUN, - ..., - times = 1000, - strategy = "sequential", - verbose = FALSE -) +test_permutation(.data, FUN, ..., times = 1000, strategy = "sequential") } \arguments{ \item{.data}{A manynet-consistent network (see -\code{\link[manynet:coerce_graph]{manynet::as_tidygraph()}}), or a list of such networks. When a list is +\code{\link[manynet:as_tidygraph]{manynet::as_tidygraph()}}), or a list of such networks. When a list is supplied the model is fit jointly; graphs missing any predictor are dropped with a warning.} @@ -47,6 +26,14 @@ e.g. the name of the attribute.} \item{times}{Integer. Number of permutations for the null distribution. 1000 is the default; publication-ready work usually needs 1000-10000.} + +\item{strategy}{If \code{{furrr}} is installed, +then multiple cores can be used to accelerate the function. +By default \code{"sequential"}, +but if multiple cores available, +then \code{"multisession"} or \code{"multicore"} may be useful. +Generally this is useful only when \code{times} > 1000. +See \href{https://furrr.futureverse.org}{\code{{furrr}}} for more.} } \description{ These functions conduct tests of any network-level statistic: @@ -54,25 +41,28 @@ These functions conduct tests of any network-level statistic: \item \code{test_random()} performs a conditional uniform graph (CUG) test of a measure against a distribution of measures on random networks of the same dimensions. +\item \code{test_configuration()} performs a CUG test against a distribution of +measures on random networks that preserve the degree sequence +of the original network. \item \code{test_permutation()} performs a quadratic assignment procedure (QAP) test of a measure against a distribution of measures on permutations of the original network. } } \examples{ -marvel_friends <- to_unsigned(ison_marvel_relationships) -marvel_friends <- to_giant(marvel_friends) \%>\% +marvel_friends <- fict_marvel |> to_uniplex("relationship") |> + to_unsigned() |> to_giant() |> to_subgraph(PowerOrigin == "Human") -(cugtest <- test_random(marvel_friends, manynet::net_heterophily, attribute = "Attractive", +(cugtest <- test_random(marvel_friends, net_by_heterophily, attribute = "Attractive", times = 200)) # plot(cugtest) # (qaptest <- test_permutation(marvel_friends, -# manynet::net_heterophily, attribute = "Attractive", +# net_by_heterophily, attribute = "Attractive", # times = 200)) # plot(qaptest) } \seealso{ -Other models: +Other models: \code{\link{regression}} } \concept{models} diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml new file mode 100644 index 0000000..152a50b --- /dev/null +++ b/pkgdown/_pkgdown.yml @@ -0,0 +1,48 @@ +url: https://stocnet.github.io/infernet/ +development: + mode: auto +template: + bootstrap: 5 + bootswatch: superhero +authors: + James Hollway: + href: https://jameshollway.com +navbar: + structure: + left: + - home + - intro + - reference + - news + right: + - search + - github + - cran + components: + home: + icon: fa-home fa-lg + href: index.html + aria-label: Go to home + reference: + text: Function Overview + href: reference/index.html + news: + text: News + href: news/index.html + github: + icon: "fab fa-github fa-lg" + href: https://github.com/stocnet/infernet + aria-label: View on Github + cran: + icon: "fab fa-r-project" + href: https://cloud.r-project.org/package=infernet + aria-label: View on CRAN +reference: + - title: "Tests" + desc: "Functions for testing network statistics against a null distribution:" + contents: + - starts_with("test") + - title: "Regression" + desc: "Functions for regressing networks on nodal and dyadic covariates:" + contents: + - regression diff --git a/tests/testthat.R b/tests/testthat.R index b1f9df7..f5a8c68 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,6 +1,4 @@ library(testthat) library(infernet) -stocnet_theme("default") test_check("infernet") -devtools::test_coverage(pkg = "infernet", type = "tests") diff --git a/tests/testthat/helper-infernet.R b/tests/testthat/helper-infernet.R new file mode 100644 index 0000000..54f16f1 --- /dev/null +++ b/tests/testthat/helper-infernet.R @@ -0,0 +1,133 @@ +# Shared fixtures and helpers. +# +# Every fixture is seeded and carries real signal, so that each family converges +# and the reference comparison is not testing noise against noise. +# `times` stays small everywhere: these tests check the estimator, not the +# precision of the null distribution. + +# ---- fixtures --------------------------------------------------------------- + +qap_net_gaussian <- function(n = 30, seed = 101) { + set.seed(seed) + age <- stats::runif(n, 20, 60) + m <- outer(age, age, function(a, b) 0.05 * a - 0.03 * b) + + matrix(stats::rnorm(n^2, sd = 0.5), n, n) + diag(m) <- 0 + manynet::mutate(manynet::as_tidygraph(m), + Age = age, + Cit = stats::rpois(n, 5), + Grp = rep(c("a", "b", "c"), length.out = n)) +} + +qap_net_binary <- function(n = 30, seed = 102) { + set.seed(seed) + age <- stats::runif(n, 20, 60) + p <- stats::plogis(outer(age, age, function(a, b) 0.06 * (a - 40) - 0.04 * (b - 40))) + m <- matrix(stats::rbinom(n^2, 1, p), n, n) + diag(m) <- 0 + manynet::mutate(manynet::as_tidygraph(m), + Age = age, + Grp = rep(c("a", "b", "c"), length.out = n)) +} + +qap_net_count <- function(n = 30, seed = 103) { + set.seed(seed) + age <- stats::runif(n, 20, 60) + lambda <- exp(outer(age, age, function(a, b) 0.02 * (a - 40) - 0.01 * (b - 40))) + m <- matrix(stats::rpois(n^2, lambda), n, n) + diag(m) <- 0 + manynet::mutate(manynet::as_tidygraph(m), Age = age) +} + +qap_net_zip <- function(n = 30, seed = 104) { + set.seed(seed) + age <- stats::runif(n, 20, 60) + lambda <- exp(outer(age, age, function(a, b) 0.02 * (a - 40) - 0.01 * (b - 40))) + m <- matrix(stats::rpois(n^2, lambda) * stats::rbinom(n^2, 1, 0.7), n, n) + diag(m) <- 0 + manynet::mutate(manynet::as_tidygraph(m), Age = age) +} + +qap_net_undirected <- function(n = 24, seed = 105) { + set.seed(seed) + age <- stats::runif(n, 20, 60) + m <- outer(age, age, function(a, b) 0.02 * (a + b)) + + matrix(stats::rnorm(n^2, sd = 0.5), n, n) + m <- (m + t(m)) / 2 + diag(m) <- 0 + manynet::mutate(manynet::as_tidygraph(m, twomode = FALSE), Age = age) +} + +qap_net_twomode <- function(seed = 106) { + set.seed(seed) + sw <- manynet::ison_southern_women + manynet::mutate(sw, Att = stats::runif(manynet::net_nodes(sw))) +} + +# ---- reference data --------------------------------------------------------- + +# Rebuilds the dyad-level data frame that the engine fits, so that a baseline +# 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) { + 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" + } + pred <- make_qap_data(y = ml$mydata[[parsed$dependent]], + x = ml$mydata[parsed$main], + diag = diag, mode = mode) + names(pred)[names(pred) == "yv"] <- parsed$dependent + list(pred = pred, formula = ml$formula, parsed = parsed) +} + +# ---- expectations ----------------------------------------------------------- + +# The p-value matrices are the product of the permutation loop, so their shape +# is the contract every estimator has to meet, whatever it fits underneath. +expect_qap_shape <- function(fit, coefs) { + testthat::expect_s3_class(fit, "net_regression") + testthat::expect_named(fit$coefficients, coefs) + for (el in c("lower", "larger", "abs")) { + m <- fit[[el]] + testthat::expect_equal(dim(m), c(2L, length(coefs)), + info = paste("dim of", el)) + testthat::expect_equal(rownames(m), c("perm_coefs", "perm_t"), + info = paste("rownames of", el)) + testthat::expect_equal(colnames(m), coefs, info = paste("colnames of", el)) + finite <- m[!is.na(m)] + testthat::expect_true(all(finite >= 0 & finite <= 1), + info = paste(el, "outside [0, 1]")) + } + invisible(fit) +} + +# ---- version-tolerant warnings ---------------------------------------------- + +# `manynet::snet_warn()` only began raising a catchable warning condition in +# manynet 2.3.2. Before that it printed a cli alert that the default +# `snet_verbosity = "quiet"` suppressed, and signalled nothing. CI installs the +# CRAN version, so a test that expects a warning has to know which it has. +# Probe the behaviour rather than the version string: a development build can +# carry the version without the behaviour. +snet_warn_signals <- function() { + isTRUE(tryCatch({ + manynet::snet_warn("probe") + FALSE + }, warning = function(w) TRUE)) +} + +# Asserts the warning where manynet raises one, and otherwise just evaluates the +# expression, so that whatever the test asserts about the returned value still +# runs. The behaviour a warning accompanies is always asserted separately. +expect_snet_warning <- function(object, regexp) { + if (snet_warn_signals()) { + testthat::expect_warning(object, regexp) + } else { + testthat::skip(paste("manynet", utils::packageVersion("manynet"), + "does not signal snet_warn() conditions")) + } +} diff --git a/tests/testthat/test-model_tests.R b/tests/testthat/test-model_tests.R new file mode 100644 index 0000000..360dc9f --- /dev/null +++ b/tests/testthat/test-model_tests.R @@ -0,0 +1,64 @@ +# Making sure the tests family of functions works as intended. +# Ported from {migraph}, which these functions supersede. + +marvel_friends <- manynet::to_uniplex(manynet::fict_marvel, "relationship") |> + manynet::to_giant() |> manynet::to_unsigned() |> + manynet::to_subgraph(PowerOrigin == "Human") + +cugtest <- test_random(marvel_friends, + netrics::net_by_heterophily, + attribute = "Attractive", + times = 200) +cugtest2 <- test_random(marvel_friends, + netrics::net_by_betweenness, + times = 200) + +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$diag) + expect_equal(cugtest$cmode, "edges") + expect_type(cugtest$plteobs, "double") + expect_type(cugtest$pgteobs, "double") + expect_equal(cugtest$reps, 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 + expect_equal(round(cugtest2$plteobs), 1) + expect_equal(round(cugtest2$pgteobs), 0) + expect_s3_class(cugtest2, "network_test") +}) + +qaptest <- test_permutation(marvel_friends, + netrics::net_by_heterophily, + attribute = "Attractive", + times = 200) + +test_that("test_permutation works", { + expect_equal(as.numeric(qaptest$testval), -0.85714, tolerance = 0.001) + 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_s3_class(qaptest, "network_test") +}) + +test_that("test_configuration works", { + testthat::skip_on_os("linux") + configtest <- test_configuration(marvel_friends, + netrics::net_by_heterophily, + attribute = "Attractive", + times = 200) + expect_s3_class(configtest, "network_test") + expect_equal(as.numeric(configtest$testval), -0.85714, tolerance = 0.001) + expect_type(configtest$plteobs, "double") # NB: Stochastic + expect_type(configtest$pgteobs, "double") # NB: Stochastic + expect_length(configtest$testdist, 200) # NB: Stochastic +}) + +test_that("print.network_test prints and returns its input invisibly", { + expect_output(print(cugtest), "CUG Test Results") + expect_invisible(print(cugtest)) + expect_identical(withVisible(print(cugtest))$value, cugtest) +}) diff --git a/tests/testthat/test-net_regression.R b/tests/testthat/test-net_regression.R index 8b1344e..bf61f99 100644 --- a/tests/testthat/test-net_regression.R +++ b/tests/testthat/test-net_regression.R @@ -104,15 +104,21 @@ test_that("net_regression fits on a list of graphs", { # ---- list-of-graphs: drop graphs missing a predictor, with warning --------- -test_that("graphs missing a predictor are dropped with a warning", { +test_that("graphs missing a predictor are dropped", { g1 <- make_weighted_net(n = 8, seed = 1) g2 <- manynet::as_tidygraph(matrix(stats::rnorm(8^2), 8, 8)) gs <- list(A = g1, B = g2) - expect_warning( - fit <- net_regression(weight ~ sim(Age), gs, times = 10), - regexp = "Dropping" - ) + fit <- suppressWarnings(net_regression(weight ~ sim(Age), gs, times = 10)) expect_s3_class(fit, "net_regression") + expect_length(unique(fit$pred$nv), 1L) +}) + +test_that("dropping a graph warns", { + g1 <- make_weighted_net(n = 8, seed = 1) + g2 <- manynet::as_tidygraph(matrix(stats::rnorm(8^2), 8, 8)) + expect_snet_warning( + net_regression(weight ~ sim(Age), list(A = g1, B = g2), times = 10), + "Dropping") }) @@ -140,3 +146,33 @@ test_that("method = 'qapy' runs and flags the nullhyp on the fit", { control = list(method = "qapy")) expect_equal(fit$nullhyp, "qapy") }) + + +# ---- tertius --------------------------------------------------------------- + +test_that("tertius() accepts a quoted and an unquoted summary function", { + g <- make_weighted_net() + quoted <- net_regression(weight ~ tertius(Age, "mean"), g, times = 10) + unquoted <- net_regression(weight ~ tertius(Age, mean), g, times = 10) + bare <- net_regression(weight ~ tertius(Age), g, times = 10) + expect_equal(quoted$coefficients, unquoted$coefficients) + expect_equal(quoted$coefficients, bare$coefficients) +}) + +test_that("tertius() sum differs from mean, and rejects anything else", { + g <- make_weighted_net() + mean_fit <- net_regression(weight ~ tertius(Age, "mean"), g, times = 10) + sum_fit <- net_regression(weight ~ tertius(Age, "sum"), g, times = 10) + expect_false(isTRUE(all.equal(mean_fit$coefficients, sum_fit$coefficients))) + expect_error(net_regression(weight ~ tertius(Age, "median"), g, times = 10), + "mean") +}) + + +# ---- messaging ------------------------------------------------------------- + +test_that("a missing attribute names what is available", { + g <- make_weighted_net() + expect_error(net_regression(weight ~ ego(Nope), g, times = 10), + "Age") +}) diff --git a/tests/testthat/test-qap_control.R b/tests/testthat/test-qap_control.R new file mode 100644 index 0000000..fcb98bd --- /dev/null +++ b/tests/testthat/test-qap_control.R @@ -0,0 +1,77 @@ +# `control` is merged over the defaults by name, so a name that is not a control +# would be added silently and the option it was meant to set would keep its +# default, with nothing to say so. + +FORM <- weight ~ ego(Age) + alter(Age) + +test_that("an unknown control name is rejected, with the nearest match", { + g <- qap_net_gaussian(n = 15) + expect_error( + net_regression(FORM, g, times = 5, control = list(strateggy = "sequential")), + "strategy") + expect_error( + net_regression(FORM, g, times = 5, control = list(nonsense = 1)), + "nonsense") +}) + +test_that("an unnamed control entry is rejected", { + g <- qap_net_gaussian(n = 15) + expect_error( + net_regression(FORM, g, times = 5, control = list("sequential")), + "named") +}) + +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()$strategy, "sequential") + expect_equal(.resolve_control()$family, "auto") +}) + +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") +}) + +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("method = 'qapy' 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") + # 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_true(all(is.na(spp$lower[, "(Intercept)"]))) +}) + +test_that("a single predictor falls back from qapspp to qapy", { + g <- qap_net_gaussian(n = 20) + fit <- net_regression(weight ~ ego(Age), g, times = 10, + control = list(seed = 1, method = "qap")) + # Double semi-partialling residualises a predictor against the others, and + # with one predictor there are none. + expect_equal(fit$nullhyp, "qapy") +}) + +test_that("mode 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, + control = list(seed = 1, + mode = "undirected"))$mode, + "undirected") + loops <- net_regression(FORM, g, times = 5, + control = list(seed = 1, diag = TRUE)) + expect_true(loops$diag) + expect_equal(nrow(loops$pred), 15 * 15) +}) diff --git a/tests/testthat/test-qap_estimators.R b/tests/testthat/test-qap_estimators.R new file mode 100644 index 0000000..fe9d8c4 --- /dev/null +++ b/tests/testthat/test-qap_estimators.R @@ -0,0 +1,254 @@ +# 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. +# +# 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 +# inference is the novel part, the point estimates are not. Second, the result +# meets the shape contract in `expect_qap_shape()`. + +COEFS3 <- c("(Intercept)", "ego Age", "alter Age", "sim Age") +FORM <- weight ~ ego(Age) + alter(Age) + sim(Age) +FORM_B <- . ~ ego(Age) + alter(Age) + sim(Age) + + +# ---- gaussian -------------------------------------------------------------- + +test_that("gaussian baseline matches lm() on the same dyads", { + g <- qap_net_gaussian() + ref <- qap_reference_data(FORM, g) + lm_fit <- stats::lm(ref$formula, data = ref$pred) + + fit <- net_regression(FORM, g, times = 10, control = list(seed = 1)) + expect_qap_shape(fit, COEFS3) + expect_equal(unname(fit$coefficients), unname(stats::coef(lm_fit))) + expect_equal(unname(fit$t), unname(summary(lm_fit)$coefficients[, 3])) + expect_equal(fit$r.squared, summary(lm_fit)$r.squared) + expect_equal(fit$adj.r.squared, summary(lm_fit)$adj.r.squared) + expect_equal(fit$family, "gaussian") +}) + +test_that("gaussian with HC3 keeps the coefficients and changes the t values", { + g <- qap_net_gaussian() + plain <- net_regression(FORM, g, times = 10, control = list(seed = 1)) + robust <- net_regression(FORM, g, times = 10, + control = list(seed = 1, use_robust_errors = TRUE)) + expect_qap_shape(robust, COEFS3) + expect_equal(robust$coefficients, plain$coefficients) + expect_false(isTRUE(all.equal(robust$t, plain$t))) + expect_true(robust$robust_se) +}) + + +# ---- binomial and poisson -------------------------------------------------- + +test_that("binomial baseline matches glm() on the same dyads", { + g <- qap_net_binary() + ref <- qap_reference_data(FORM_B, g) + glm_fit <- stats::glm(ref$formula, data = ref$pred, + family = stats::binomial()) + + fit <- net_regression(FORM_B, g, times = 10, + control = list(seed = 1, family = "binomial")) + expect_qap_shape(fit, COEFS3) + expect_equal(unname(fit$coefficients), unname(stats::coef(glm_fit))) + expect_equal(unname(fit$t), unname(summary(glm_fit)$coefficients[, 3])) +}) + +test_that("family = 'auto' picks binomial for a binary outcome", { + auto <- net_regression(FORM_B, qap_net_binary(), times = 10, + control = list(seed = 1)) + named <- net_regression(FORM_B, qap_net_binary(), times = 10, + control = list(seed = 1, family = "binomial")) + expect_equal(auto$family, "binomial") + expect_equal(auto$coefficients, named$coefficients) +}) + +test_that("poisson baseline matches glm() on the same dyads", { + g <- qap_net_count() + ref <- qap_reference_data(FORM, g) + glm_fit <- stats::glm(ref$formula, data = ref$pred, family = stats::poisson()) + + fit <- net_regression(FORM, g, times = 10, + control = list(seed = 1, family = "poisson")) + expect_qap_shape(fit, COEFS3) + expect_equal(unname(fit$coefficients), unname(stats::coef(glm_fit))) +}) + + +# ---- negative binomial and zero-inflated Poisson --------------------------- + +test_that("negbin baseline matches MASS::glm.nb() on the same dyads", { + skip_if_not_installed("MASS") + g <- qap_net_count() + ref <- qap_reference_data(FORM, g) + nb <- suppressWarnings(MASS::glm.nb(ref$formula, data = ref$pred)) + + fit <- suppressWarnings( + net_regression(FORM, g, times = 10, + control = list(seed = 1, family = "negbin"))) + expect_qap_shape(fit, COEFS3) + expect_equal(unname(fit$coefficients), unname(stats::coef(nb))) + expect_equal(fit$theta, nb$theta) +}) + +test_that("zip baseline matches pscl::zeroinfl() and names its coefficients", { + skip_if_not_installed("pscl") + g <- qap_net_zip() + ref <- qap_reference_data(FORM, g) + zi <- pscl::zeroinfl(ref$formula, data = ref$pred, dist = "poisson") + + fit <- net_regression(FORM, g, times = 10, + control = list(seed = 1, family = "zip")) + # A backticked name here used to break double semi-partialling, which looks a + # column up by the unquoted predictor name. + expect_qap_shape(fit, COEFS3) + expect_false(any(grepl("`", names(fit$coefficients), fixed = TRUE))) + expect_equal(unname(fit$coefficients), unname(zi$coefficients$count)) + expect_equal(unname(fit$zi_coefficients), unname(zi$coefficients$zero)) +}) + + +# ---- 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", { + skip_if_not_installed("lme4") + g <- qap_net_gaussian() + ref <- qap_reference_data(FORM, g) + mixed <- suppressMessages(suppressWarnings( + lme4::lmer(build_internal_formula(ref$formula, ris = TRUE), + data = ref$pred))) + + fit <- suppressMessages(suppressWarnings( + net_regression(FORM, g, times = 10, + control = list(seed = 1, random_intercept_sender = TRUE)))) + expect_qap_shape(fit, COEFS3) + expect_equal(unname(fit$coefficients), + unname(summary(mixed)$coefficients[, 1])) + expect_named(fit$random.intercepts, "sv") +}) + +test_that("crossed sender and receiver intercepts do not abort the run", { + skip_if_not_installed("lme4") + # Residualising a predictor against both intercepts is often singular, and + # `lmer()` then stops with "Downdated VtV is not positive definite". That is + # a step towards the null distribution, so it falls back rather than aborting. + fit <- suppressMessages(suppressWarnings( + net_regression(FORM, qap_net_gaussian(), times = 10, + control = list(seed = 1, + random_intercept_sender = TRUE, + random_intercept_receiver = TRUE)))) + expect_qap_shape(fit, COEFS3) +}) + +test_that("binomial and poisson random intercepts run", { + skip_if_not_installed("lme4") + bin <- suppressMessages(suppressWarnings( + net_regression(FORM_B, qap_net_binary(), times = 10, + control = list(seed = 1, family = "binomial", + random_intercept_sender = TRUE)))) + expect_qap_shape(bin, COEFS3) + pois <- suppressMessages(suppressWarnings( + net_regression(FORM, qap_net_count(), times = 10, + control = list(seed = 1, family = "poisson", + random_intercept_sender = TRUE)))) + 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") +}) + + +# ---- 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_reproducibility.R b/tests/testthat/test-qap_reproducibility.R new file mode 100644 index 0000000..afb1060 --- /dev/null +++ b/tests/testthat/test-qap_reproducibility.R @@ -0,0 +1,64 @@ +# A permutation result that cannot be reproduced cannot be published, so the +# seed is part of the interface rather than an implementation detail. + +FORM <- weight ~ ego(Age) + alter(Age) + sim(Age) + +test_that("the same seed gives the same p-values", { + g <- qap_net_gaussian() + a <- net_regression(FORM, g, times = 30, control = list(seed = 99)) + b <- net_regression(FORM, g, times = 30, control = list(seed = 99)) + expect_equal(a$lower, b$lower) + expect_equal(a$larger, b$larger) + expect_equal(a$abs, b$abs) + expect_equal(a$coefficients, b$coefficients) +}) + +test_that("a different seed gives a different null distribution", { + # Pure noise, so that the tallies land inside (0, 1) and can differ. A strong + # effect drives every p-value to 0 or 1 under any seed, which would make this + # assertion pass for the wrong reason. + set.seed(77) + n <- 25 + m <- matrix(stats::rnorm(n^2), n, n) + diag(m) <- 0 + g <- manynet::mutate(manynet::as_tidygraph(m), Age = stats::runif(n, 20, 60)) + a <- net_regression(FORM, g, times = 60, control = list(seed = 99)) + b <- net_regression(FORM, g, times = 60, control = list(seed = 7)) + # The observed model does not depend on the seed; only the null does. + expect_equal(a$coefficients, b$coefficients) + expect_false(isTRUE(all.equal(a$lower, b$lower))) +}) + +test_that("an outer set.seed() reproduces a run with no seed control", { + g <- qap_net_gaussian() + set.seed(5); a <- net_regression(FORM, g, times = 30) + set.seed(5); b <- net_regression(FORM, g, times = 30) + expect_equal(a$lower, b$lower) +}) + +test_that("a parallel plan gives the same answer as a sequential one", { + skip_on_cran() + g <- qap_net_gaussian() + seq <- net_regression(FORM, g, times = 30, control = list(seed = 99)) + par <- net_regression(FORM, g, times = 30, + control = list(seed = 99, strategy = "multisession")) + # `furrr_options(seed = TRUE)` and `future.seed = TRUE` are what make this + # true; without them the plan would change the result. + expect_equal(par$lower, seq$lower) + expect_equal(par$coefficients, seq$coefficients) +}) + +test_that("the future plan is restored after a run", { + before <- class(future::plan()) + net_regression(FORM, qap_net_gaussian(), times = 10, + control = list(seed = 1, strategy = "multisession")) + expect_equal(class(future::plan()), before) +}) + +test_that("the test family reproduces under an outer seed", { + g <- qap_net_gaussian(n = 18) + set.seed(3); a <- test_random(g, netrics::net_by_density, times = 30) + set.seed(3); b <- test_random(g, netrics::net_by_density, times = 30) + expect_equal(a$testdist, b$testdist) + expect_equal(a$pgteobs, b$pgteobs) +}) diff --git a/tests/testthat/test-qap_shapes.R b/tests/testthat/test-qap_shapes.R new file mode 100644 index 0000000..f821c93 --- /dev/null +++ b/tests/testthat/test-qap_shapes.R @@ -0,0 +1,96 @@ +# Properties of the dependent network -- modes, directedness, loops -- must be +# respected in the data the engine fits and in the permutations it draws. +# Each shape below was wrong at some point, so each is asserted by counting the +# dyads that reach the model rather than by checking that the call returns. + +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_equal(nrow(fit$pred), 20 * 19) +}) + +test_that("an undirected network contributes each dyad once", { + g <- qap_net_undirected(n = 24) + expect_false(manynet::is_directed(g)) + 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_equal(nrow(fit$pred), 24 * 23 / 2) +}) + +test_that("a two-mode network contributes every cell of the incidence matrix", { + sw <- qap_net_twomode() + dims <- manynet::net_dims(sw) + fit <- net_regression(. ~ ego(Att) + alter(Att), sw, times = 10, + control = list(seed = 1)) + # An incidence matrix is rectangular, and has no diagonal to drop. Reading it + # as square wrapped past the last column and invented dyads. + expect_equal(nrow(fit$pred), dims[1] * dims[2]) + expect_named(fit$coefficients, c("(Intercept)", "ego Att", "alter Att")) +}) + +test_that("RMPerm() permutes a rectangular matrix without erroring", { + m <- matrix(seq_len(18 * 14), 18, 14) + p <- RMPerm(m) + expect_equal(dim(p), c(18L, 14L)) + # A permutation relabels the nodes, so it moves cells but keeps the multiset. + expect_setequal(as.vector(p), as.vector(m)) +}) + +test_that("RMPerm() keeps the row and column order aligned for a square matrix", { + set.seed(4) + m <- matrix(seq_len(36), 6, 6) + p <- RMPerm(m) + expect_equal(dim(p), c(6L, 6L)) + expect_setequal(as.vector(p), as.vector(m)) + # One order for both margins, so the diagonal stays the diagonal. + expect_setequal(diag(p), diag(m)) +}) + +test_that("dist() and sim() read each mode of a two-mode network separately", { + sw <- qap_net_twomode() + ml <- convertToMatrixList(. ~ ego(Att) + alter(Att) + dist(Att) + sim(Att), + sw, advise = FALSE) + ego <- ml$mydata[["ego Att"]] + alt <- ml$mydata[["alter Att"]] + expect_equal(ml$mydata[["dist Att"]], abs(ego - alt)) + denom <- max(abs(ego - alt)) + expect_equal(ml$mydata[["sim Att"]], abs(1 - abs(ego - alt) / denom)) +}) + +test_that("a list of networks is pooled, and one missing a predictor is dropped", { + good <- list(qap_net_gaussian(n = 18, seed = 1), + qap_net_gaussian(n = 18, seed = 2)) + fit <- net_regression(weight ~ ego(Age), good, times = 10, + control = list(seed = 1)) + expect_equal(nrow(fit$pred), 2 * 18 * 17) + expect_length(unique(fit$pred$nv), 2L) + + bare <- manynet::as_tidygraph(matrix(stats::rnorm(18^2), 18, 18)) + dropped <- suppressWarnings( + net_regression(weight ~ ego(Age), list(good[[1]], bare, good[[2]]), + times = 10, control = list(seed = 1))) + expect_equal(nrow(dropped$pred), 2 * 18 * 17) +}) + +test_that("dropping a network from a list warns", { + good <- qap_net_gaussian(n = 18, seed = 1) + bare <- manynet::as_tidygraph(matrix(stats::rnorm(18^2), 18, 18)) + expect_snet_warning( + net_regression(weight ~ ego(Age), list(good, bare), times = 10, + control = list(seed = 1)), + "Dropping") +}) + +test_that("a missing dyad is dropped from the model", { + g <- qap_net_gaussian(n = 20) + m <- manynet::as_matrix(g) + m[1, 2] <- NA + holed <- manynet::mutate(manynet::as_tidygraph(m), + Age = manynet::node_attribute(g, "Age")) + fit <- net_regression(weight ~ ego(Age), holed, times = 10, + control = list(seed = 1)) + expect_equal(nrow(fit$pred), 20 * 19 - 1) +})