# SciPy vs. statsmodels for Bioinformatics Statistics SciPy gives you fast, self-contained statistical functions that operate directly on arrays, while statsmodels fits full statistical models with covariates, interactions, and detailed diagnostic output. The rule of thumb for choosing between them in bioinformatics is simple. If you need a quick two-group test, a rank test across thousands of features, or a distance calculation, reach for SciPy. If you need to model expression as a function of genotype, batch, age, and treatment together, or you need [confidence intervals](/blog/data-analysis/confidence-interval-formula-examples), post-hoc comparisons, and proper multiple-testing correction, reach for statsmodels. Most real Python for bioinformatics projects end up using both, often in the same notebook. This guide explains what each library does on its own terms, compares them on the criteria that actually decide the choice, works through concrete bioinformatics scenarios, and lists the mistakes that show up most often in real analyses. ## SciPy vs. statsmodels for Bioinformatics Statistics at a Glance ![Decision path asking whether covariates or multiple testing are needed, leading to SciPy or statsmodels.](/article-images/illustrations/scipy-vs-statsmodels-bioinformatics-library-choice.webp) The difference between SciPy and statsmodels for bioinformatics comes down to scope. SciPy answers narrow questions about arrays. statsmodels answers broad questions about data-generating processes. The table below is the core of the comparison, and the rest of the article expands each row. | Criterion | SciPy | statsmodels | |--|--|--| | Core role | Scientific computing toolkit with a large statistics module | Econometrics-grade [statistical modeling](/blog/guides/understanding-model-assumptions-in-statistical-analysis) library | | Typical unit of work | A single function call on arrays | A fitted model object with a full summary | | Input format | NumPy arrays, lists, pandas columns | DataFrames plus a model formula or design matrices | | Regression with covariates | Not supported in `scipy.stats` | Central feature (OLS, GLM, mixed models) | | Categorical predictors and interactions | Not supported | First-class via the formula API | | Simple hypothesis tests | Extensive (`ttest_ind`, `mannwhitneyu`, `fisher_exact`, `chi2_contingency`) | Available but usually more than you need | | Post-hoc comparisons | Not provided | Tukey HSD and others in `statsmodels.stats.multicomp` | | Multiple-testing correction | Basic (false discovery control in recent versions) | Full toolkit (`multipletests` with many methods) | | Probability distributions | Very large catalog (200-plus distributions) | Distribution-related helpers tied to models | | Speed on simple tests | Very fast, easy to vectorize | Slower per call, designed for model fitting | | Statistical output | Test statistic and p-value | Coefficients, standard errors, p-values, confidence intervals, fit diagnostics | | Best fit in bioinformatics | Per-gene screens, enrichment tests, distances, clustering inputs | Covariate-adjusted models, ANOVA designs, mixed models, power analysis | ## What SciPy Does SciPy is a general scientific computing library built on NumPy, the array library that underlies nearly all numerical work in Python. Its statistics module, `scipy.stats`, contains a large collection of probability distributions and hypothesis tests designed to run on plain arrays. That design choice explains almost everything about how SciPy behaves in a bioinformatics pipeline. ### The shape of a SciPy analysis A SciPy [statistical analysis](/blog/guides/statistical-analysis) is a function call. You hand it data, it hands back a statistic and a p-value. A p-value is the probability of observing a result at least as extreme as yours if the null hypothesis, the assumption of no effect, is true. There is no model object, no stored design, and no built-in summary beyond the core numbers. ```python from scipy import stats ko = [7.1, 7.4, 6.9, 7.3, 7.0] # log2 expression, knockout wt = [6.2, 6.5, 6.1, 6.4, 6.3] # log2 expression, wild type result = stats.ttest_ind(ko, wt, equal_var=False) print(result.statistic, result.pvalue) ``` This is a Welch two-sample t-test, which compares the means of two independent groups without assuming equal variances. Note the explicit `equal_var=False`. The default in `ttest_ind` is the Student version, which assumes equal variances, so you must opt in to Welch when variances might differ. This default is one of the most consequential details in `scipy.stats` and one of the most commonly missed. ### The breadth of scipy.stats The module covers far more than t-tests: - **Distribution tests and transformations**: `shapiro` for normality, `levene` for variance homogeneity, `boxcox` for variance-stabilizing transforms. - **Nonparametric tests**: `mannwhitneyu` (Wilcoxon rank-sum), `wilcoxon` (signed-rank for paired data), `kruskal` (nonparametric one-way ANOVA). - **Count and categorical methods**: `chisquare`, `chi2_contingency`, `fisher_exact`, `binomtest`. - **Correlation**: `pearsonr`, `spearmanr`, `kendalltau`. - **Pairwise distances and linkage**: `pdist`, `squareform`, and `scipy.cluster.hierarchy.linkage`, which feed directly into dendrograms for sample clustering. - **Sparse matrices**: `scipy.sparse` handles cell-by-gene matrices in single-cell work, where most entries are zero. - **Multivariate helpers**: `scipy.spatial.distance` and singular value decomposition via `scipy.linalg` and `sklearn`-adjacent workflows for dimensionality reduction. In single-cell RNA-seq, the scanpy ecosystem leans on SciPy heavily. Sparse expression matrices, k-nearest-neighbor graphs, and the per-cluster Wilcoxon rank-sum tests used to rank marker genes all run through SciPy primitives. If you load a [10x Genomics](/blog/guides/10x-genomics) matrix with tens of millions of nonzero counts, it is almost certainly stored as a `scipy.sparse` matrix under the hood. ### When to use scipy Use SciPy when the statistical question is about the relationship between two or three arrays and nothing else. Concretely: - Testing whether a single gene differs between two conditions, after log transformation. - Running a rank-sum test per gene across thousands of genes, vectorized over a NumPy array. - Computing a Fisher exact test for a 2x2 enrichment table. - Calculating pairwise Euclidean or correlation distances between samples for a heatmap. - Fitting a distribution to a peptide intensity or fragment-length distribution and pulling quantiles from it. Because each call is cheap and stateless, SciPy scales beautifully when you repeat the same simple test many times. A loop or an `apply` over 20,000 genes calling `mannwhitneyu` is a normal pattern. The same loop calling a statsmodels regression per gene is possible but noticeably heavier. One more reason to know SciPy deeply: many p-values elsewhere in Python are computed by SciPy. When a downstream tool reports a p-value, understanding `scipy.stats` distributions and their survival functions (the probability of exceeding a value, available as `.sf` on any distribution object) helps you verify what the number means. ## What statsmodels Does statsmodels is a statistical modeling library in the tradition of R's model formula interface. Where SciPy asks "are these two arrays different," statsmodels asks "how does an outcome depend on several variables at once, and how certain are we about each contribution". ### The formula API and design matrices The heart of statsmodels for bioinformatics work is the formula interface, powered by the patsy formula syntax. You write a model the way you would in R: ```python import statsmodels.formula.api as smf model = smf.ols("log2_expr ~ genotype + batch + age", data=df).fit() print(model.summary()) ``` Here `log2_expr` is the response variable, and `genotype`, `batch`, and `age` are covariates, the explanatory variables you want to adjust for. A batch effect is a systematic technical difference between processing groups, such as sequencing lane or library prep date, and adjusting for it is one of the most common reasons bioinformaticians move from SciPy to statsmodels. The formula API automatically expands categorical variables into indicator terms and can express interactions, for example `genotype * treatment`, which tests whether the treatment effect differs by genotype. The fitted object returned by `.fit()` is where statsmodels earns its keep. `model.summary()` prints coefficients, standard errors, t-values, p-values, 95 percent confidence intervals, R-squared, and diagnostic statistics in one table. With SciPy you would have to compute most of those quantities yourself. ### Regression families that matter for biological data Biological data rarely meets the assumptions of ordinary least squares, the linear model that assumes continuous, roughly normally distributed residuals. statsmodels provides [generalized linear models](/blog/data-analysis/generalized-linear-models-explained) (GLMs), which extend regression to response variables with other distributions: - **Gaussian**: standard linear regression, appropriate for log-transformed intensities and continuous phenotypes. - **Binomial with a logit link**: [logistic regression](/blog/data-analysis/logistic-regression-definition-formula) for binary outcomes, such as responder versus nonresponder, or variant present versus absent. - **Poisson and negative binomial**: for count data. Sequencing read counts per gene are counts, and the negative binomial family accommodates the overdispersion (more variance than a Poisson model predicts) that is characteristic of RNA-seq data. ```python import statsmodels.api as sm import statsmodels.formula.api as smf nb = smf.glm( "counts ~ genotype + batch", data=df, family=sm.families.NegativeBinomial() ).fit() ``` For a single gene or a handful of genes, this is a legitimate way to fit a count model with covariates. For genome-wide differential expression, dedicated tools handle the hard parts better, a point covered in the limitations section below. statsmodels also carries tools that have no SciPy equivalent: - **ANOVA with multiple factors**: `anova_lm` decomposes sums of squares for multi-factor designs with interactions, while SciPy's `f_oneway` handles only a one-way design. - **Post-hoc comparisons**: `pairwise_tukeyhsd` runs Tukey's honestly significant difference test across all pairs of group means, controlling the family-wise error rate, the chance of at least one false positive among the comparisons. - **Mixed models**: `mixedlm` fits linear mixed-effects models with random effects, which is how you analyze repeated measurements on the same animal, plate, or donor without pretending the observations are independent. - **Power analysis**: `statsmodels.stats.power` computes the sample size needed to detect an effect. For example, `TTestIndPower().solve_power(effect_size=0.8, alpha=0.05, power=0.8)` returns the sample size per group needed to detect a large standardized effect (Cohen's d of 0.8) with 80 percent power at a significance level of 0.05. - **Robust and quantile regression**: alternatives to least squares when residuals are heavy-tailed, which happens with proteomics and metabolomics data. ### The multiple-testing toolkit The single most cited statsmodels feature in genomics is `multipletests` in `statsmodels.stats.multitest`. When you test thousands of genes at a conventional alpha of 0.05, you expect thousands of false positives by chance alone. The Benjamini-Hochberg procedure controls the false discovery rate (FDR), the expected proportion of false positives among your called significant results, and `multipletests` implements it plus Bonferroni, Holm, Benjamini-Yekutieli, and others in one call: ```python from statsmodels.stats.multitest import multipletests reject, q_values, _, _ = multipletests(raw_pvalues, alpha=0.05, method="fdr_bh") ``` The returned `q_values` are adjusted p-values, often called q-values, which you can threshold directly. Recent versions of SciPy added `scipy.stats.false_discovery_control` for Benjamini-Hochberg and related adjustments, but statsmodels offers the wider menu of methods and is the more standard reference point for this step. ## Head-to-Head on the Criteria That Decide Real Choices ### Hypothesis testing depth Both libraries perform standard tests, but they are not interchangeable. `scipy.stats` gives you a broader catalog of individual tests with minimal ceremony. statsmodels gives you deeper machinery around the tests it supports, including Wald tests, likelihood ratio tests, and contrast testing on fitted coefficients. If your question maps onto a named test, SciPy is usually the shorter path. If your question requires defining a linear contrast across coefficients, for example testing whether the knockout effect equals the rescue effect, statsmodels is the tool that supports it. ### Models, covariates, and design This is the sharpest dividing line. SciPy's hypothesis tests have no concept of a covariate. If your two groups differ in sex composition, [sequencing depth](/knowledge/diagnostics/molecular/how-to-calculate-sequencing-depth-and-coverage-for-your-ngs-run), or age, a plain SciPy test cannot absorb those differences, and confounding, where a variable influences both the grouping and the outcome, will bias the result. statsmodels puts covariates into the model explicitly. In practice, most animal studies and human cohort analyses have at least one variable that must be adjusted for, which is why regression-based analysis is the default in mature genomics workflows. ### Data structure and ergonomics SciPy functions take arrays. statsmodels models take DataFrames plus formulas. If your data lives in a pandas DataFrame with mixed types and categorical columns (both libraries integrate cleanly with pandas), statsmodels formulas eliminate a class of ordering and encoding bugs, because column names travel with the data into the summary output. If your data lives in a NumPy matrix and the operation is elementwise or pairwise, SciPy's interface is faster to write and faster to run. ### Speed and scale For simple tests, SciPy is markedly faster and easier to vectorize. Calling `stats.mannwhitneyu` across 20,000 genes or `stats.pearsonr` across thousands of pairs is routine. statsmodels model fitting has overhead per call: formula parsing, design matrix construction, and iterative fitting for GLMs. Fitting 20,000 negative binomial GLMs in a naive loop can take minutes to hours depending on the data, and the common workaround is to run such jobs in parallel or to move the problem to a dedicated method-of-moments or shrunken-dispersion tool. Choose per-feature statsmodels fits only when the number of features is small or you genuinely need covariate adjustment for each one. ### Output, diagnostics, and reproducibility A SciPy call returns a statistic and a p-value. A statsmodels fit returns a result object with coefficient estimates, standard errors, confidence intervals, and fit diagnostics, plus methods for predictions, residuals, and contrasts. When a reviewer asks for an effect size with a confidence interval, statsmodels hands it to you. When a collaborator asks why a sample is an outlier, statsmodels residual methods help. SciPy will not object to any of this, but you will be doing more arithmetic yourself. A second, quieter benefit of statsmodels is reproducibility of intent. A formula string like `outcome ~ genotype + batch + (1 | litter)` documents the analysis design in the code itself. A stack of SciPy calls requires separate documentation of what was adjusted for and how. ## Worked Scenarios These four scenarios cover the decisions that come up most often. All code is illustrative and uses small placeholder-style values in structure, with realistic magnitudes for the field. ### Scenario 1: One gene, two groups, log-scale expression You measured log2-transformed expression of a candidate gene in five knockout and five wild-type samples and want a straightforward comparison. This is the textbook case for when to use scipy. ```python from scipy import stats ko = [7.1, 7.4, 6.9, 7.3, 7.0] wt = [6.2, 6.5, 6.1, 6.4, 6.3] res = stats.ttest_ind(ko, wt, equal_var=False) # Welch's t-test print(res.statistic, res.pvalue) ``` Three details matter. First, work on the log scale, because expression intensities are right-skewed and log transformation makes t-test assumptions far more plausible. Second, use Welch by passing `equal_var=False` unless you have a reason to assume equal variances. Third, a SciPy t-test gives you no effect size, so report the mean difference on the log2 scale yourself, which corresponds to a fold change after back-transformation. If the same question must be answered per gene for the whole transcriptome, this exact call is what you vectorize, followed by FDR correction from statsmodels. ### Scenario 2: Pathway enrichment from a 2x2 table Your differentially expressed gene list contains 12 genes from a pathway of interest, and 88 of your 4,000 detected genes in that pathway were tested. Wait, let me set the table cleanly. Suppose the background is 20,000 detected genes, of which 200 belong to the pathway. Your significant list contains 100 genes, of which 12 are pathway members. The 2x2 [contingency table](/blog/data-analysis/what-is-a-contingency-table) is pathway membership crossed with significance. ```python from scipy import stats # rows: in pathway, not in pathway # columns: significant, not significant table = [[12, 188], [88, 19712]] odds_ratio, p = stats.fisher_exact(table, alternative="greater") print(odds_ratio, p) ``` `fisher_exact` handles 2x2 tables exactly, meaning it computes probabilities from the hypergeometric distribution rather than relying on large-sample approximations. Use it when any expected cell count is small, which is the norm in enrichment testing when a pathway is rare in the background. For larger tables, `scipy.stats.chi2_contingency` applies the chi-square test of independence, and `scipy.stats.hypergeom` gives you direct access to the underlying distribution for gene-set arithmetic beyond 2x2. ### Scenario 3: Expression with batch effects and covariates You suspect a gene differs by genotype, but samples were run in two sequencing batches and donors range widely in age. A plain two-group test would conflate genotype with batch. This is the canonical case for statsmodels. ```python import statsmodels.formula.api as smf model = smf.ols("log2_expr ~ genotype + batch + age + sex", data=df).fit() print(model.summary()) # The genotype coefficient now estimates the genotype effect # holding batch, age, and sex constant, with a 95 percent CI. ci = model.conf_int().loc["genotype[T.KO]"] ``` The formula syntax expands categorical variables automatically, and the summary reports the genotype coefficient with its standard error and confidence interval. If instead of log expression your response is raw read counts, switch to a negative binomial GLM as shown earlier. If samples were measured repeatedly over time, for example weekly weights on the same animals, use a mixed model so that within-animal correlation is modeled rather than ignored: ```python import statsmodels.formula.api as smf mixed = smf.mixedlm("weight ~ time + treatment", df, groups=df["animal_id"]).fit() ``` Mixed models, also called random-effects models, include a per-cluster (here, per-animal) random intercept, which prevents the inflated significance that comes from treating repeated measures as independent samples. SciPy cannot do this at all. ### Scenario 4: A genome-wide screen that combines both libraries ![Flowchart from expression matrix through vectorized SciPy t-tests, FDR correction, and a sorted results table.](/article-images/illustrations/scipy-vs-statsmodels-bioinformatics-screen-workflow.webp) The most common real pattern uses both libraries in sequence. Vectorized SciPy tests generate raw p-values for every feature, then statsmodels handles the correction and the modeling around the screen. ```python import numpy as np import pandas as pd from scipy import stats from statsmodels.stats.multitest import multipletests # gene x sample matrix of log2 expression, plus a group label vector t_stats, p_vals = stats.ttest_ind(expr_ko, expr_wt, axis=1, equal_var=False) reject, qvals, _, _ = multipletests(p_vals, alpha=0.05, method="fdr_bh") results = pd.DataFrame({ "gene": gene_names, "log2FC": expr_ko.mean(axis=1) - expr_wt.mean(axis=1), "p": p_vals, "q": qvals, }).sort_values("q") ``` Here SciPy does the per-gene testing at speed, statsmodels does the FDR correction, and pandas assembles the result table. This combined pattern, simple test plus multipletests, is a workhorse in labs that analyze omics data in Python, and it also mirrors how rank-based marker screens in single-cell pipelines are organized. ## Where These Libraries Sit in a Python for Bioinformatics Stack Neither library operates alone. A typical Python for bioinformatics environment stacks them like this: 1. **NumPy** provides the array foundation. 2. **pandas** provides labeled tables for sample sheets, feature annotation, and results. 3. **Biopython** handles the biological file formats and online resources. A standard Biopython tutorial covers parsing FASTA and GenBank records, running and parsing BLAST output, working with alignments, and retrieving records from NCBI Entrez. Biopython itself is not a statistics library, but everything it parses eventually becomes the arrays and DataFrames that SciPy and statsmodels consume. 4. **SciPy** supplies tests, distributions, distances, sparse matrices, and clustering primitives. 5. **statsmodels** supplies covariate-adjusted modeling, post-hoc tests, mixed models, and multiple-testing correction. 6. **Domain frameworks** such as scanpy for single-cell and PyDESeq2 for [RNA-seq differential expression](/knowledge/bioinformatics/rna-seq-differential-expression-deseq2-edger) build on these primitives and add field-specific methods like dispersion shrinkage. The practical workflow is often: parse sequences or import a count matrix with Biopython or a domain tool, assemble a sample-sheet DataFrame with pandas, screen features with vectorized SciPy tests, correct and model with statsmodels, and report from a pandas DataFrame. Learning order matters less than knowing which layer solves which problem. If you are starting from scratch, learn the SciPy tests first because they get you analyzing data within a day, then add statsmodels as soon as your questions involve covariates or multiple testing. ## Common Mistakes and Limitations These are the errors that appear repeatedly in real projects, with the fix for each. **Using the default t-test when variances differ.** `scipy.stats.ttest_ind` defaults to Student's t-test, which pools variance across groups. If group variances differ, which is common when one condition has a subset of genes or cells strongly induced, Welch's test is the safer choice. Pass `equal_var=False` explicitly. **Treating paired data as independent.** Measurements taken on the same animal before and after treatment, or the same cell line across passages, are paired. The correct tools are `scipy.stats.ttest_rel` or `scipy.stats.wilcoxon` for paired designs, or a mixed model with the subject as a random effect when there are more than two time points. Ignoring pairing inflates significance because it treats correlated observations as independent evidence. **Applying the chi-square test to sparse contingency tables.** The chi-square approximation requires adequate expected counts in each cell. When cells are small, which is typical in pathway and variant enrichment problems, use `fisher_exact` (2x2) or exact and Monte Carlo-based methods instead. A chi-square p-value from a table with expected counts under 5 is not trustworthy. **Skipping multiple-testing correction in screens.** Running 20,000 uncorrected tests at alpha 0.05 yields roughly 1,000 nominal hits under the null by chance. Every genome-wide screen needs FDR control, most commonly Benjamini-Hochberg via `multipletests`. Report q-values, not raw p-values, in any screen output. **Interpreting a p-value as an effect size.** A p-value measures evidence against the null hypothesis, not the magnitude of the effect. With large sample sizes, trivial differences become significant. statsmodels summaries report coefficients and confidence intervals, and you should always report the effect (a fold change, an odds ratio, a difference in means) alongside the p-value. **Fitting raw counts with a linear model.** Ordinary least squares assumes continuous, homoscedastic errors. Raw sequencing counts are counts, with variance tied to the mean and often overdispersed. Log-transform plus a linear model is an acceptable approximation for many intensity-style measurements, but for counts a negative binomial GLM in statsmodels, or a dedicated count-based differential expression tool, is the more defensible route. **Expecting statsmodels to be a full RNA-seq differential expression engine.** statsmodels can fit a negative binomial model per gene, but dedicated tools estimate gene-wise dispersion with shrinkage across genes, handle zero inflation and size-factor normalization, and apply empirical Bayes moderation. A per-gene statsmodels GLM without those refinements will generally be less well calibrated for genome-wide inference. Use statsmodels for covariate modeling and multiple testing within a pipeline, and use a dedicated differential expression method for the genome-wide count modeling step. **Ignoring diagnostics.** A significant p-value from a model that violates its assumptions is noise with a p-value attached. Plot residuals against fitted values, check for influential points, and verify that batch and group are not perfectly confounded. No p-value can rescue a design where every knockout was sequenced in one lane and every control in another, because genotype and batch are then algebraically inseparable. statsmodels' residual and influence methods make these checks straightforward, and SciPy's normality and variance tests help with the simpler checks. **Vectorizing incorrectly.** A Python loop over 20,000 genes calling a test one at a time can be dramatically slower than the array-oriented call with `axis=1` or an equivalent vectorized formulation. When a loop is unavoidable, at minimum keep it in NumPy arrays rather than pandas `apply` on object columns, and consider multiprocessing across features. **Reinventing known results without checking the platform.** NCBI maintains extensive documentation and tooling for sequence data standards and retrieval, and established statistical practice for omics data is codified in standard textbooks of [biostatistics](/blog/data-analysis/what-is-biostatistics) and genomics. Before writing a novel statistical procedure, confirm that a maintained method already exists, because most common omics questions have a well-validated answer with published calibration. A final limitation that applies to this whole article: statistical libraries produce numbers, and the interpretation of those numbers for a specific experiment, animal study, or clinical dataset requires domain judgment about the design and the biology, ideally from a collaborator who has seen the raw data. ## Frequently Asked Questions ### Do I need statsmodels if I already use SciPy? Yes, as soon as your analysis involves covariates, interactions, post-hoc comparisons, mixed models, or multiple-testing correction beyond basic FDR. SciPy remains the faster tool for simple array-level tests, and the two libraries coexist in virtually every serious bioinformatics environment. ### Which library should I learn first for Python for bioinformatics? Learn SciPy first. Its tests run directly on arrays, require no formula syntax, and cover the everyday tasks of comparing groups, computing correlations, and testing enrichment. Add statsmodels when you begin adjusting for batch effects or fitting regression models. ### Can statsmodels replace DESeq2 for RNA-seq differential expression? No. statsmodels can fit negative binomial GLMs per gene, but dedicated tools add dispersion shrinkage, size-factor normalization, and empirical Bayes moderation that are important for well-calibrated genome-wide inference. Use a dedicated differential expression tool for the genome-wide step and statsmodels for covariate-heavy or custom modeling questions. ### How do I correct p-values when testing thousands of genes? Use `statsmodels.stats.multitest.multipletests` with `method="fdr_bh"` for Benjamini-Hochberg false discovery rate control, which is the standard choice for genome-wide screens. Bonferroni is more conservative and appropriate when any false positive is unacceptable. ### Does Biopython do statistics? Biopython is focused on sequence parsing, alignment interfaces, and NCBI data retrieval, not statistical modeling. Statistical work in a Biopython-driven project happens in SciPy and statsmodels after Biopython hands over parsed data. ### Which library is faster for testing thousands of genes? SciPy. Its test functions are lightweight array operations that vectorize cleanly with `axis` arguments, while statsmodels fits a full model object per call. For genome-wide screens, run vectorized SciPy tests first, then apply statsmodels correction and modeling to the results.