Metadata-Version: 2.4
Name: biplotpy
Version: 0.4.0
Summary: Biplots for Python, ported from R and validated against it: classical and HJ, canonical/CVA and HE plots, logistic (incl. constrained), STATIS, CA/MCA/CCA/RDA, PCoA/MDS with bootstrap, FA, PLS (incl. binary), GGE/AMMI and more, with their graphics
Author-email: Carlos Freire-Cadme <cfreirec@unemi.edu.ec>
License: GPL-3.0-or-later
Project-URL: Homepage, https://github.com/cfreirec/biplotpy
Project-URL: Source, https://github.com/cfreirec/biplotpy
Project-URL: Issues, https://github.com/cfreirec/biplotpy/issues
Classifier: Programming Language :: Python :: 3
Classifier: License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.8
Requires-Dist: pandas>=1.4
Requires-Dist: matplotlib>=3.5
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# biplotpy

Biplots for Python, with their graphics, ported from R and validated
numerically against the R packages they come from: **MultBiplotR**
(Vicente-Villardón, Vicente-González & Frutos-Bernal), GGEBiplots,
agricolae, smacof, vegan, ade4, FactoMineR, ca, pls, psych, GPArotation,
DistatisR, biplotbootGUI, MANOVA.RM and heplots.

| Method | Function |
|---|---|
| Classical biplots: GH / CMP (α=0), JK / RMP (α=1), SQRT (α=0.5) | `gh_biplot`, `jk_biplot`, `sqrt_biplot`, `pca_biplot(alpha=)` |
| HJ-Biplot (Galindo, 1986) | `hj_biplot` |
| Quality of representation, contributions (relative and to the factor) | `.summary()`, `.qualities()` |
| Canonical / MANOVA biplot with confidence circles (t, Bonferroni, Hotelling, χ²), MANOVA (Wilks, Pillai, Hotelling-Lawley, Roy), ANOVAs, LDA classification | `canonical_biplot`, `canonical_variate_analysis` |
| Bootstrap of the PCA biplot | `pca_bootstrap` |
| Supplementary individuals and continuous variables | `project_rows`, `add_continuous_variables` |
| Calibrated (predictive) axes | `biplot_scales`, `calibrated_axis`, `plot(calibrated=True)` |
| 12 initial transformations | `initial_transform` |
| Logistic biplots: binary, ordinal (cumulative logits), nominal | `binary_logistic_biplot`, `ordinal_logistic_biplot`, `nominal_logistic_biplot` |
| STATIS-ACT, Dual STATIS, Canonical STATIS-ACT | `statis`, `dual_statis`, `canonical_statis` |
| Dynamic biplot (Egido & Galindo) | `dynamic_biplot` |
| Metric (ratio) unfolding by SMACOF | `unfolding` |
| GGE biplot and AMMI | `gge_model`, `ammi` |
| DISTATIS (DistatisR) | `distatis` |
| Clustering and disjoint PCA (CDPCA, biplotbootGUI) | `cdpca` |
| Correspondence analysis (MultBiplotR CA, ca, FactoMineR) and multiple CA (FactoMineR MCA) | `ca`, `mca` |
| Canonical correspondence analysis (MultBiplotR CCA, vegan cca) and redundancy analysis (vegan rda, scalings 1-3) | `cca`, `rda` |
| Principal coordinates (MultBiplotR, ade4 dudi.pco), 17 binary coefficients, continuous and Gower distances | `pcoa`, `binary_distances`, `continuous_distances`, `gower_distances` |
| Metric and non-metric MDS (smacof::mds: ratio, interval, ordinal with primary/secondary/tertiary ties) | `mds` |
| External binary logistic biplot and ridge logistic regression (MultBiplotR) | `external_logistic_biplot`, `ridge_logistic` |
| Canonical analysis of distances (MultBiplotR) | `canonical_distance_analysis` |
| Factor-analysis biplot with GPA rotations (varimax, quartimax, equamax, parsimax, entropy, quartimin, oblimin) | `fa_biplot`, `rotate` |
| NIPALS biplot with missing values (ade4 nipals), sparse and truncated NIPALS | `nipals`, `nipals_biplot` |
| PLS regression biplot (pls::plsr) | `pls_regression`, `pls_biplot` |
| Co-inertia analysis (ade4) | `coinertia` |
| Polychoric correlations (psych) and polychoric ordinal biplot (MultBiplotR PolyBiplot) | `polychoric`, `polychoric_biplot` |
| Mardia's test; heteroscedastic MANOVA with WTS/MATS and parametric or wild bootstrap (MANOVA.RM); HE plots (heplots) | `mardia`, `manova_wts`, `he_matrices`, `he_ellipses`, `he_plot` |
| Multi-factor MANOVA (type II and III SSP, car::Anova) and HE plots with one ellipse per term | `manova_terms`, `he_plot_terms` |
| Constrained (restricted) logistic biplots, binary and ordinal (MultBiplotR) | `constrained_logistic_biplot`, `constrained_ordinal_logistic_biplot` |
| External ordinal logistic biplot and ridge cumulative-logit regression (MultBiplotR AddOrdVars2Biplot, OrdinalLogisticFit) | `external_ordinal_logistic_biplot`, `ridge_ordinal_logistic` |
| PLS regression with a binary response (PLSR1Bin) and PLS for binary data (BinaryPLSR) | `plsr1_bin`, `binary_plsr`, `plot_plsbin` |
| Log-frequency biplot and weighted ALS biplot (CrissCross) | `log_frequency_biplot`, `criss_cross` |
| Bootstrap of principal coordinates (BootstrapDistance) and of SMACOF MDS | `bootstrap_pcoa`, `bootstrap_mds`, `plot_bootstrap_ordination` |
| R-compatible random numbers and R's `optim(method = "CG")` | `RRNG`, `optim_cg` |

```python
import biplotpy as bp
iris = bp.load_iris()
hj = bp.hj_biplot(iris.iloc[:, :4], dim=3)
hj.summary()["col_qualities"]
hj.plot(groups=iris["Species"])
cb = bp.canonical_biplot(iris.iloc[:, :4], iris["Species"])
cb.manova; cb.plot(radius="bonferroni")
```

## Validation

`tests/test_vs_multbiplotr.py` (fixtures generated by `tests/fixtures/make_fixtures.R`)
checks, up to the sign of each axis: coordinates, eigenvalues, inertia, scale
factor, structure correlations and contributions of 12 classical biplots
(α = 0, 0.5, 1, HJ; five transformations); 11 initial transformations; the
canonical biplot (group, variable, individual markers, weights, eigenvalues,
Wilks F, ANOVAs, four radii, contributions, qualities); supplementary
variables; calibrated scales (port of R's `pretty()`); MANOVA tables against
R's `summary.manova`; LDA accuracy against `MASS::lda`.

Phase 2 (`tests/test_phase2.py`, fixtures from `fixtures/make_phase2.R`):

| Method | Reference | Max. abs. difference |
|---|---|---|
| `statis` | MultBiplotR::StatisBiplot; RV also vs ade4::statis | < 1e-13 |
| `canonical_statis` | MultBiplotR::CanonicalStatisBiplot | < 1e-14 |
| `dual_statis` | core lines of MultBiplotR::DualStatisBiplot (the function itself is unfinished) | < 1e-13 |
| `unfolding` (ratio) | smacof::unfolding 2.1-7 (same iterations, stress-1, configurations) | < 1e-13 |
| `gge_model` | GGEBiplots::GGEModel, 4 SVP x 3 centerings | < 1e-13 |
| `ammi` | agricolae::AMMI (ANOVA, Gollob table, scores) | < 1e-10 |
| `ordinal_logistic_fit` | MASS::polr | < 3e-7 |
| `nominal_logistic_fit` | nnet::multinom | < 2e-6 |
| `binary_logistic_biplot` | R optim (BFGS) on the same penalized objective | same optimum (eta < 1e-3) |
| `dynamic_biplot` | algebraic identities (reference occasion reproduces the biplot; mean of trajectories = reference) | < 1e-10 |

Phase 3 (`tests/test_phase3.py`, fixtures from `fixtures/make_phase3.R`):

| Method | Reference | Max. abs. difference |
|---|---|---|
| `mardia` | psych::mardia | < 1e-12 |
| `manova_wts` | MANOVA.RM 0.6.0 (statistics sourced from the package) | WTS, MATS, chi-square p < 1e-12; bootstrap p within Monte Carlo error |
| `he_matrices`, `roy_crit` | car::Anova (SSP), heplots Roy.crit | < 1e-12 |
| `distatis` | DistatisR 1.1.2 (MFA, SUMPCA, NUCLEAR norms) | < 1e-14 (NUCLEAR ~1e-8, see note) |
| `cdpca` | biplotbootGUI::CDpca from the same start | identical partitions, F < 1e-12 |

Note on DISTATIS with the nuclear norm: DistatisR sums the square roots of the
eigenvalues larger than 2 eps; a numerically zero eigenvalue of a
double-centred matrix can pass that threshold in one library and not in the
other, which moves the results in the 8th decimal.

Simulated recovery: STATIS recovers a common configuration and down-weights a
pure-noise table; unfolding recovers an exact configuration (stress < 1e-3);
logistic biplots recover the true probabilities, with correlation increasing
with n and p (0.81 at 100 x 10, 0.97 at 400 x 80 for the binary model).

Caveat: the deviance "tests" of the logistic biplots use df = dim as in
MultBiplotR, but the row scores are estimated from the same data, so the
p-values are optimistic; use them as descriptive fit indices. Joint ML with a
very small penalization overfits (near separation); use penalization >= 0.5
for small tables.

Phase 4 (`tests/test_phase4.py`, fixtures from `fixtures/make_phase4.R`):

| Method | Reference | Max. abs. difference |
|---|---|---|
| `ca` | MultBiplotR::CA (alpha 0, 0.5, 1), ca::ca | < 1e-14 |
| `mca` | FactoMineR::MCA (eigenvalues, individuals, categories, eta2) | < 1e-13 |
| `cca` | MultBiplotR::CCA (all scores), vegan::cca (eigenvalues, total inertia) | < 1e-13 |
| `rda` | vegan::rda + scores (scalings 1, 2, 3; scaled and unscaled) | < 1e-12 |
| `pcoa`, `binary_distances` | MultBiplotR::PrincipalCoordinates, ade4::dudi.pco, BinaryDistances (8 coefficients) | < 1e-14 |
| `mds` | smacof::mds (ratio, interval, ordinal primary and secondary): same iterations | < 1e-10 |
| `external_logistic_biplot` | MultBiplotR::ExternalBinaryLogisticBiplot | < 1e-10 |
| `canonical_distance_analysis` | MultBiplotR::CanonicalDistanceAnalysis | < 1e-12 |
| `rotate` | GPArotation (GPForth, GPFoblq) | < 1e-8 |
| `fa_biplot` | MultBiplotR::FA.Biplot loadings; psych::factor.scores (Bartlett, Thurstone) | < 1e-10 |
| `nipals` | ade4::nipals with missing values | < 1e-10 |
| `pls_regression` | pls::plsr (oscorespls) | < 1e-6 (iterative) |
| `coinertia` | ade4::coinertia | < 1e-10 |
| `polychoric` | psych::polychoric | < 1.2e-4 (psych's optimizer tolerance) |

Phase 5, version 0.4.0 (`tests/test_phase5.py`, fixtures from `fixtures/make_phase5.R`):

| Method | Reference | Max. abs. difference |
|---|---|---|
| `constrained_logistic_biplot` | MultBiplotR::ConstrainedLogisticBiplot (coefficients, coordinates, deviances, pseudo R2, loadings, sensitivity, specificity, supplementary variables) | < 1e-13 |
| `constrained_ordinal_logistic_biplot` | MultBiplotR::ConstrainedOrdinalLogisticBiplot | < 1e-12 |
| `external_ordinal_logistic_biplot` | MultBiplotR::AddOrdVars2Biplot (see note 1) | < 1e-12 |
| `plsr1_bin` | MultBiplotR::PLSR1Bin, algorithm 1 (see note 2) | < 1e-14 |
| `binary_plsr` | MultBiplotR::BinaryPLSR, same seed: same random starts and same CG optimizer | intercepts < 1e-15; scores and loadings < 1e-4 (note 3) |
| `optim_cg` | R's optim(method = "CG") | same path (Rosenbrock: same counts, parameters < 1e-12) |
| `log_frequency_biplot` | MultBiplotR::LogFrequencyBiplot (CrissCross) with the same random start | < 1e-13 |
| `bootstrap_pcoa` | MultBiplotR::BootstrapDistance, same seed, sampling and permutation | < 1e-12 |
| `bootstrap_mds` | the same residual bootstrap in R on smacof::mds (note 4) | < 1e-13 |
| `manova_terms` | car::Anova on a two-factor unbalanced MANOVA, types II and III | < 1e-12 |

Notes. (1) `AddOrdVars2Biplot` converts the factors to integers and
`RidgeOrdinalLogistic` then stops ("The dependent variable must be
ordinal"); the reference is computed with that conversion disabled. (2) In
R >= 4.2 `PLSR1Bin` stops inside `RidgeBinaryLogistic` (`class(y) ==
"formula"` on a matrix); the reference uses `class(y)[1]`. Its algorithm 2
refers to undefined objects and is not ported. (3) BinaryPLSR alternates
conjugate-gradient fits from random starts; biplotpy draws the same starts
(R's generator) and runs R's CG algorithm, but the path amplifies
differences in the last digits of the arithmetic (BLAS summation order), so
results agree to about 1e-4. (4) `BootstrapSmacof` stops (it reads `Sol$Dhat`,
which `SMACOF` does not return); `bootstrap_mds` resamples the residuals
between disparities and configuration distances, truncates negative
dissimilarities at 0, refits from the original configuration and aligns by
orthogonal Procrustes; the R reference runs the same scheme on smacof::mds.
`MultBiplotR::SimpleProcrustes`, used by `BootstrapDistance`, applies the
transpose of the optimal rotation (it fixes reflections but turns small
rotations the wrong way); `bootstrap_pcoa` uses the optimal rotation by
default and `procrustes="MultBiplotR"` reproduces R. In the constrained
binary logistic biplot the deviance tests use `dim` degrees of freedom (R
uses 1) and the percentages of correct classification threshold the
probability at 0.5 (R thresholds the linear predictor).

## Graphics

Every method has a plot (`examples/gallery.py` renders all 49 figures into
`examples/gallery/`): classical biplots with arrows or calibrated axes,
group concentration ellipses, convex hulls, alpha-bags and densities
(ggbiplot, factoextra, biplotEZ), quality colouring, correlation circle,
scree and contribution plots, 3-D biplot (bpca), bootstrap regions;
canonical/CVA biplot with confidence circles, classification regions and
calibrated axes (biplotEZ); HE plots (heplots); binary/ordinal logistic
biplots and nominal prediction regions; STATIS interstructure, compromise
and trajectories; dynamic biplot; unfolding; DISTATIS; CDPCA; the ten GGE
views of GGEBiplots (which-won-where, mean vs stability,
discriminativeness vs representativeness, rankings...); AMMI1/AMMI2;
CA/MCA maps, CCA and RDA triplots; PCoA/MDS with supplementary variables;
external logistic biplot; canonical analysis of distances; co-inertia;
FA, NIPALS and PLS biplots; polychoric biplot; constrained and external
ordinal logistic biplots; PLSR1Bin and BinaryPLSR biplots; log-frequency
biplot; bootstrap clouds of PCoA and MDS; HE plots of several factors. Their geometry is checked by
tests (calibrated axes read the CVA group means exactly; ellipses and
alpha-bags have their nominal coverage) and every figure is rendered in the
test suite.

## Differences from MultBiplotR (corrections)

See `biplotpy.FIXES`: within-groups standardization divides by the standard
deviation (MultBiplotR divides by the variance); the initial logit is applied
(MultBiplotR leaves the data untransformed); supplementary rows are projected
after transforming them (MultBiplotR projects them untransformed); STATIS
honours the requested transformation; Dual STATIS is complete; the binary
logistic objective and gradient are consistent; NIPALS, Coinertia, PLSR,
FA scores, PolyBiplot, BinaryDistances and weighted PCoA are corrected (16
entries in total).

## Not ported

Already available in Python, so not duplicated: CA, MCA, MFA, FAMD
(`prince`); CCA, RDA, PCoA, PERMDISP/betadisper (`scikit-bio`); PPML with
high-dimensional fixed effects (`pyfixest`); Ledoit-Wolf shrinkage and LDA
(`scikit-learn`); nearest PSD matrix (`statsmodels.stats.correlation_tools.cov_nearest`);
orthogonal Procrustes and generalized eigenproblems (`scipy`).

Not included, and why:

| Method (R) | Reason |
|---|---|
| Graphical interfaces (BiplotGUI, dynBiplotGUI, multibiplotGUI, biplotbootGUI's GUI) | Interactive Tcl/Tk applications, not statistical methods; every method behind them is in biplotpy with a matplotlib plot. |
| Functional PCA biplot (`BiplotFPCA`) | Needs a functional-data layer (basis expansions and smoothing of curves) that biplotpy does not have; it would be a separate library. |
| TetraDualStatis (Dual STATIS on tetrachoric correlations of binary tables) | Its help page has no example ("Not yet") and it depends on psych's tetrachoric estimates; there is no documented use to validate against. |
| Tucker and Parafac (ThreeWay) | Three-way decompositions with their own estimation and identification problems; available in Python in `tensorly`, so not duplicated. |
| STATICO | Depends on ade4's `ktab`/`statico` machinery (partial triadic analysis of co-inertia tables); porting ade4's k-table classes is outside the scope of a biplot library. |
| Estimation with mirt (`BinaryLogBiplotMirt`, ordinal factor analysis) | Item-response estimation by marginal ML with numerical integration (mirt); Python has no validated equivalent. biplotpy fits logistic biplots by penalized joint ML, which is MultBiplotR's default. |
| EM estimation of logistic biplots (`BinaryLogBiplotEM`, `OrdLogBipEM`) | Marginal ML with Gauss-Hermite quadrature, the same model as mirt; not ported for the same reason. |
| `PLSR1Bin` algorithm 2 | It refers to objects that are never defined (`c0`, `c`) and cannot run in R. |

A tutorial in Spanish with an example of each family of methods is in
`examples/tutorial_es.ipynb`.
