pySpATS

Spatial correction of field-trial phenotypes in Python. Fit genetic effects, experimental factors and a two-dimensional spline surface jointly, then export plot diagnostics and genotype effects for downstream analysis.

Version 0.3 replaces the 2025 numerical implementation. It follows the supplied R SpATS 1.0-20 mathematics and is tested against frozen outputs generated by that package. See the methodology review, validation report and migration notes.

Install

Python 3.10 or newer:

pip install -e .
# Optional plots and development tests:
pip install -e '.[plot,test]'

NumPy, SciPy and pandas are the runtime dependencies. R, scikit-learn and CHOLMOD are not required. Matplotlib is imported only when plotting.

Fit one field

from pyspats import fit_trial, PSANOVA
from pyspats.datasets import load_wheatdata

plots = load_wheatdata()  # Actual 330-plot wheat trial supplied with R SpATS
result = fit_trial(
    data=plots,
    response="yield",
    genotype="geno",
    genotype_as_random=True,
    spatial=PSANOVA("col", "row", nseg=(10, 10)),
    random=["R", "C"],  # Categorical row/column effects in addition to the smooth
)

print(result)                         # Includes convergence status
print(result.summary())              # Variances and exact effective dimensions
print(result.heritability)            # Available for random genotypes only
plot_results = result.to_frame()      # Original order/index, plus diagnostics
varieties = result.genotype_predictions()  # Effects, conditional SEs, BLUE/BLUP label
plot_results.to_csv("wheat_plots.csv", index=False)
varieties.to_csv("wheat_genotypes.csv", index=False)

SpATS(...) remains an equivalent constructor that fits immediately. The keyword-only fit_trial(...) entry point makes scripts easier to read.

Choose what the model means

Input	Meaning
genotype_as_random=False	Estimate the tested genotype effects without shrinkage (BLUEs). No heritability calculation.
genotype_as_random=True	Estimate genotype variance and shrink noisy genotype effects toward zero (BLUPs).
fixed=["treatment", "flowering_date"]	Additive covariates. Numeric columns are continuous; categorical/string columns are factors.
random=["block"]	Independent categorical random effects. Numeric IDs are still treated as categories.
weights="precision"	Gaussian observation variance is dispersion/precision. Zero excludes a plot from estimation.
offset="known_effect"	Known contribution on the link scale, not a covariate whose slope should be estimated.
geno_decomp="population"	Separate random genotype variances by population; each genotype belongs to one population.

Cast numeric treatment codes to category if they represent categories. Coordinates must be numeric with variation along both axes. Do not combine independent fields by inventing adjacent coordinates. Use physical distances if plot spacing differs and distance-based diagnostics need physical meaning. The spline treats each axis separately; its basis uses the supplied coordinate range and the specified segment count.

Use fixed=["treatment"] for treatment effects; a covariate perfectly confounded with genotype or the spatial polynomial is not estimable. The fitter raises an error rather than silently dropping terms or adding a ridge.

Configure the surface

from pyspats import PSANOVA, SAP, SpATSControl

surface = PSANOVA("col", "row", nseg=(20, 30), nest_div=(2, 3), center=True)
control = SpATSControl(tolerance=1e-7, max_iter=500)
# Supply spatial=surface and control=control to fit_trial.

nseg is the number of equal-width segments, not the number of coefficients. nest_div reduces the smooth-by-smooth interaction basis only and must divide nseg. This is R's explicit speed/complexity tradeoff. It changes the model; choose it deliberately. PSANOVA has five variance components and second-order penalties. SAP(...) selects the R two-penalty alternative; SAP(..., ANOVA=True) selects the four-penalty alternative. A coordinate tuple ("col", "row") selects default PSANOVA. No setting is silently ignored.

Understand the output

• fitted_values: expected phenotype including genotype, spatial, fixed and random contributions (and offset).
• residuals: observed minus fitted, on the response scale.
• spatial_trend: fitted spatial polynomial and smooth contributions, centered over fitted plots using precision weights; on the link scale for GLMMs.
• adjusted_values: Gaussian observed phenotype minus that centered spatial trend. It retains genotype, treatments, blocks and noise. It is not a genotype BLUE/BLUP or an independent error-free input to a second analysis.
• genotype_predictions(): coefficients labeled BLUE or BLUP, with conditional SEs. Random effects exclude the intercept; fixed effects are evaluated at zero for other model columns. These are not marginal treatment-adjusted means.
• effective_dim: actual effective dimensions, not numbers of parameters.
• psi: residual dispersion. var_comp: random-component variances, not smoothing penalties. residual_df: n minus fixed rank minus random ED.
• deviance: R-compatible restricted working objective, omitting a constant; not Gaussian residual sum of squares or a general-purpose AIC statistic.
• converged, history, outer_history: fit status and numerical history.

The generalized heritability denominator accounts for confounding with all fixed effects. It is often the number of genotypes minus one, but need not be. It is a measure for this trial and fitted model, not a universal trait constant.

Missing responses and zero-weight plots do not enter estimation; they receive predictions if predictors and fitted factor levels are available. Missing predictors yield NaN. No input DataFrame or global random seed is modified. The used_for_fit output column makes exclusions explicit.

Three standard tables

plot_values = result.adjusted_plots()
genotype_values = result.genotype_means()
variance = result.variance_partition()

# Use identical overrides when comparing plot and genotype outputs:
reference = {"treatment": "irrigated", "flowering_time": 65}
# These must be fitted fixed covariates with valid levels/values.
# plot_values = result.adjusted_plots(reference=reference)
# genotype_values = result.genotype_means(reference=reference)

The shared default reference averages fixed categorical levels equally, uses numeric covariate means over observed positive-weight plots, averages the spatial surface equally over those plots, and sets other random effects to zero. Known offsets use their mean; reference_offset=0 (or another value) overrides this. Every genotype uses the same reference, regardless of its occupied locations. These outputs currently support Gaussian models.

genotype_means() returns environment, trait, genotype, predicted, se, n_obs, and type. Predictions include the reference phenotype level; random genotypes retain shrinkage. SEs include covariance with the reference, conditional on estimated variance parameters. They describe mean estimation/prediction error, not the spread of future plot observations.

adjusted_plots() returns original-order rows and index, with environment, trait, plot_id (zero-based input position), genotype, observed, fitted, residual, predicted_at_reference, adjusted, and used_for_fit. Adjusted = predicted_at_reference + observed - fitted. Thus all nuisance effects are standardized while genotype differences and plot residuals remain. Missing observations remain missing in adjusted; predictions are separate. Original input columns remain available in result.data and can be joined by row position even when the original index is duplicated.

Both tables store the resolved reference and convergence status in .attrs. CSV does not preserve attrs; save that metadata alongside exported tables. Environment and trait columns allow concatenation across independently fitted fields without implying that they were fitted jointly. Existing to_frame() continues to provide its spatial-only adjustment.

Variance partitioning

result = fit_trial(
    data=plots, response="yield", genotype="geno",
    spatial=PSANOVA("col", "row"), genotype_as_random=True,
    environment_id="wheat_trial",
)
partition = result.variance_partition()
spatial_details = result.variance_partition(spatial="components")
partition.to_csv("variance_partition.csv", index=False)

Both tables contain environment, trait, component, role, variance, and percent. Random-effect and residual percentages sum to 100. Ordinary random intercept factors contribute their estimated variance; splines contribute average plot-level variance, trace(Z G Z') / n, using the fitted prior coefficient covariance. This puts them on the same squared-phenotype scale. The calculation uses bounded chunks and never constructs a plot-by-plot covariance matrix. Detailed spatial rows replace the combined spatial row.

Fixed effects are listed with missing variance and percentage: like fixed terms in an lme4 model, they have no estimated variance component. This includes fixed genotypes and the unpenalized spatial polynomial. Thus the denominator is random effects plus residual, not the sample variance of the observations. The report currently supports Gaussian models. Each observed positive-weight plot counts equally; weighted residual variance is the average of psi / weight. Spatial contributions use the fitted basis and its centering setting. partition.attrs records total variance, plot count and convergence status.

environment_id labels a single fitted field; it does not introduce an environment effect. fit_trials supplies this label from its grouping key (a tuple when grouping by multiple columns), so reports can be concatenated across fields and traits. Joint multi-environment estimation is not yet supported.

Prediction and plotting

new_plots = plots.iloc[:10].copy()
predictions = result.predict(new_plots, return_se=True)
# Reuse known offsets explicitly: result.predict(new_plots, offset="known_effect")
fig = result.plot(show=False)
fig.savefig("wheat_diagnostics.png", dpi=150)

Training knots, centering and categorical encodings are reused. Prediction rejects unseen factor levels and coordinates outside the fitted field. se_link measures uncertainty in the latent predictor, conditional on fitted variance parameters; it is not a prediction interval for a future observation. Objects can be pickled for reuse within a compatible Python/package environment.

pyspats.variogram(result) computes an exact empirical variogram of response residuals, with bounded temporary memory. Pair enumeration still takes O(n²) time; use this diagnostic deliberately on large fields.

Multiple fields and traits

from pyspats import fit_trials, PSANOVA

for job in fit_trials(
    data=all_plots, by=["site", "year"], responses=["yield", "height"],
    genotype="genotype", genotype_as_random=True,
    spatial=PSANOVA("col", "row", nseg=10, nest_div=2),
    workers=1, on_error="record",
):
    if job.ok:
        print(job.key, job.response, job.model.summary())
    else:
        print(job.key, job.response, job.error or "Did not converge")

Each job is an independent single-field fit. The iterator preserves job identity and reports errors. workers>1 uses bounded process submission; place script calls inside if __name__ == "__main__":. Set BLAS thread counts before starting Python to avoid multiplying process and BLAS parallelism. Use column names for per-plot weights and offsets so grouping aligns them. See the batch CLI example for CSV exports and a machine-readable success/failure manifest.

Counts and binary traits

Pass family=poisson() or family=binomial() from pyspats.families for counts or binary/proportion responses. Binomial weights represent trial counts; Poisson inputs must be nonnegative integers. These follow R's approximate working mixed-model method. Matching R does not establish accurate uncertainty for rare events or separation. No Gaussian-style adjusted phenotype is reported for these families. Monitor both inner and outer convergence.

Development and provenance

python -m pytest -q
python examples/wheat.py
OPENBLAS_NUM_THREADS=1 PYTHONPATH=. python benchmarks/benchmark.py
# Optional regeneration; requires the supplied R SpATS 1.0-20 installed:
Rscript tests/reference/generate.R

GPL-2.0-or-later; translated from R SpATS by María Xosé Rodríguez-Álvarez, Martin Boer, Paul Eilers and Fred van Eeuwijk. Cite Rodríguez-Álvarez et al. (2018) for the spatial methodology. The wheat data are attributed in the data provenance notes.