# ANOVA in Python for Comparing Multiple Groups: A Practical Guide By the end of this guide you will be able to load a multi-group dataset, check the assumptions that ANOVA depends on, run a one-way ANOVA with `scipy` and `statsmodels`, quantify the effect size, apply a post-hoc test to find which groups differ, and correct p values across thousands of features in a typical omics table. Everything runs in a plain Python environment with a handful of open-source libraries. You will need Python 3.9 or newer, `numpy`, `pandas`, `scipy`, `statsmodels`, and optionally `pingouin` and `matplotlib`. No prior statistics training beyond a basic grasp of p values is assumed. ## Why ANOVA Shows Up Everywhere in Bioinformatics Analysis of variance (ANOVA) answers one question: do the means of two or more groups differ more than you would expect from within-group variability alone? When a study compares [a control group](/blog/guides/understanding-control-groups-in-clinical-trials-types-and-selection) against several treatment groups, ANOVA is usually the first inferential test applied to continuous measurements such as protein abundances, metabolite concentrations, methylation percentages, or qPCR expression values. The test appears constantly in published bioinformatics workflows. A proteomics package built for mass spectrometry data lists ANOVA among its core statistical algorithms alongside t tests, [principal component analysis](/blog/data-analysis/principal-component-analysis-explained), and hierarchical clustering [1]. A [DNA methylation](/blog/guides/dna-methylation) toolkit includes ANOVA among the [statistical methods](/blog/guides/statistical-methods) it offers for identifying differentially methylated CpG sites between groups [2]. A metabolomics suite built on Python and Galaxy exposes ANOVA and the Kruskal-Wallis test as its basic group-comparison methods [3]. In applied studies, researchers routinely use one-way ANOVA followed by Tukey's HSD test and Cohen's d effect sizes to compare treatment groups, running the analysis and graphics in Python [4]. A comparative metabolomics study of infertility phenotypes used one-way ANOVA in R and Python with false discovery rate correction to find altered metabolites across three patient groups [5]. The common pattern is clear. When you have one continuous measurement and three or more groups, ANOVA is the standard omnibus test, and Python has mature tools to run it. ## What You Need Installed Install the core stack with pip: ```bash python -m venv anova-env source anova-env/bin/activate pip install numpy pandas scipy statsmodels pingouin matplotlib ``` On Windows, replace the `source` line with `anova-env\Scripts\activate`. If you prefer conda: ```bash conda create -n anova-env python=3.11 numpy pandas scipy statsmodels pingouin matplotlib conda activate anova-env ``` Verify the installation: ```python import numpy, pandas, scipy, statsmodels, pingouin print(numpy.__version__, pandas.__version__, scipy.__version__, statsmodels.__version__) ``` Any recent stable versions work. If you later want to read sequence files or fetch records from NCBI as part of a larger workflow, the Biopython package (`pip install biopython`) handles sequence input and retrieval, while the statistics stack covered here handles the analysis. A full Biopython tutorial is outside the scope of this article, but the two stacks complement each other in a typical [Python for bioinformatics setup](/knowledge/bioinformatics/python-for-bioinformatics-setup). ## Step-by-Step: One-Way ANOVA in Python This walkthrough uses a simulated proteomics-style example: one protein measured across three treatment groups (control, low dose, high dose) with six biological replicates per group. The same code works unchanged for [gene expression](/blog/guides/gene-expression), metabolite concentrations, or any continuous measurement. ### Step 1: Build the dataset in long format ANOVA tools in Python expect long-format data, meaning one row per observation with a group label column and a value column. ```python import numpy as np import pandas as pd rng = np.random.default_rng(42) groups = ["control", "low_dose", "high_dose"] params = {"control": (10.0, 0.8), "low_dose": (11.2, 0.9), "high_dose": (13.5, 1.1)} rows = [] for g in groups: mean, sd = params[g] for value in rng.normal(mean, sd, size=6): rows.append({"group": g, "abundance": value}) df = pd.DataFrame(rows) print(df.head()) ``` Expected output: ``` group abundance 0 control 9.87 1 control 10.42 2 control 9.95 3 control 11.21 4 control 10.09 ``` Exact numbers depend on your NumPy version and the random draw, but the structure should match. In real work you would replace this block with `pd.read_csv("results.tsv", sep="\t")` or an import from a proteomics or methylation pipeline. ### Step 2: Look at the data before testing anything Never run ANOVA on a table you have not inspected. Compute group summaries first: ```python summary = df.groupby("group")["abundance"].agg(["count", "mean", "std"]) print(summary.round(3)) ``` Expected output (values will differ slightly): ``` count mean std group control 6.0 10.112 0.742 high_dose 6.0 13.480 1.023 low_dose 6.0 11.345 0.877 ``` Check three things. Every group should have a reasonable number of observations, the standard deviations should be broadly similar across groups, and no value should be an obvious data-entry error. A quick boxplot catches most problems: ```python import matplotlib.pyplot as plt df.boxplot(column="abundance", by="group") plt.title("Protein abundance by group") plt.suptitle("") plt.ylabel("Log2 abundance") plt.savefig("boxplot.png", dpi=150) ``` ### Step 3: Test the assumptions Classical one-way ANOVA assumes that observations are independent, that residuals are approximately normally distributed within each group, and that the groups have similar variances. Two standard checks cover the distributional assumptions. GeneQuantify, a [qPCR analysis](/blog/guides/qpcr-analysis-from-plate-layout-to-a-defensible-result) platform implemented in Python, applies exactly this pair of checks (Shapiro-Wilk for normality and Levene for variance homogeneity) to ΔCt values before further testing [6], and the practice generalizes to any continuous endpoint. Normality per group with Shapiro-Wilk: ```python from scipy import stats for g in groups: stat, p = stats.shapiro(df.loc[df["group"] == g, "abundance"]) print(f"{g}: W = {stat:.3f}, p = {p:.3f}") ``` Expected output looks like this: ``` control: W = 0.912, p = 0.461 high_dose: W = 0.958, p = 0.795 low_dose: W = 0.934, p = 0.617 ``` A p value above 0.05 means you cannot reject normality, which is what you want here. With only six observations per group the Shapiro-Wilk test has low power, so treat it as a screen rather than a guarantee. Homogeneity of variance with Levene's test, using the median-centered version, which is more robust to non-normality: ```python samples = [df.loc[df["group"] == g, "abundance"] for g in groups] stat, p = stats.levene(*samples, center="median") print(f"Levene: W = {stat:.3f}, p = {p:.3f}") ``` A p value above 0.05 supports the equal-variance assumption. If this check fails, skip ahead to Step 8 for the Welch and Kruskal-Wallis alternatives. ### Step 4: Run the ANOVA two ways The fastest route is a single function call in SciPy: ```python f_stat, p_value = stats.f_oneway(*samples) print(f"F = {f_stat:.3f}, p = {p_value:.3g}") ``` Expected output: ``` F = 26.4, p = 8.3e-06 ``` Because the simulated group means are far apart relative to within-group noise, the F statistic is large and the p value is far below 0.05. The statsmodels route gives you a full ANOVA table and is the better choice when you need sums of squares for effect sizes or want to extend to multi-factor designs: ```python import statsmodels.api as sm from statsmodels.formula.api import ols model = ols("abundance ~ C(group)", data=df).fit() table = sm.stats.anova_lm(model, typ=2) print(table) ``` Expected output: ``` sum_sq df F PR(>F) C(group) 36.871 2.0 26.410009 0.000008 Residual 10.472 15.0 NaN NaN ``` The parameters that matter here: - `abundance ~ C(group)`: the model formula. `C(group)` tells statsmodels to treat the group column as categorical. If your column name contains spaces or starts with a number, rename it first, because the formula parser will choke on it. - `typ=2`: the type of sums of squares. Type 2 is the default recommendation for models without interactions. It matters mainly in multi-factor designs. - `anova_lm(model, typ=2)`: produces the classic ANOVA table from a fitted linear model. The `PR(>F)` column is the p value. The `F` column is the ratio of between-group to within-group variance. Both routes should give the same F and p, which is a useful cross-check covered below. ### Step 5: Quantify the effect size ![Eta squared computed as 36.871 over the total sum of squares, giving 0.779, a large effect by Cohen's benchmarks.](/article-images/illustrations/python-anova-bioinformatics-eta-squared.webp) A significant p value tells you the groups differ. It does not tell you how much. Eta squared (η²) is the fraction of total variance explained by group membership, and you can compute it directly from the ANOVA table: ```python ss_group = table["sum_sq"] ss_total = table["sum_sq"].sum() eta_sq = ss_group / ss_total print(f"eta squared = {eta_sq:.3f}") ``` Expected output: ``` eta squared = 0.779 ``` Under Cohen's conventional benchmarks, η² near 0.01 is a small effect, near 0.06 a medium effect, and near 0.14 or larger a large effect. The simulated example produces a large effect because the group means were set far apart on purpose. In applied work, researchers often report Cohen's d for pairwise comparisons alongside ANOVA, as done in a recent treatment study analyzed with Python [4]. ### Step 6: Find which groups differ with Tukey's HSD A significant omnibus F test does not tell you which pairs of groups differ. Tukey's honestly significant difference (HSD) test compares every pair while controlling the family-wise error rate: ```python from statsmodels.stats.multicomp import pairwise_tukeyhsd tukey = pairwise_tukeyhsd(df["abundance"], df["group"], alpha=0.05) print(tukey) ``` Expected output: ``` Multiple Comparison of Means - Tukey HSD, FWER=0.05 ========================================================= group1 group2 meandiff p-adj lower upper reject --------------------------------------------------------- high_dose control -3.3678 0.0001 -4.4244 -2.3111 True high_dose low_dose -2.1352 0.0007 -3.1918 -1.0785 True control low_dose 1.2326 0.0238 0.1759 2.2893 True --------------------------------------------------------- ``` Read the `reject` column: `True` means that pair differs significantly at the 0.05 family-wise level. The `meandiff` column gives the difference in group means, which is the number your biology colleagues will actually care about. A common error here is passing columns in the wrong order. `pairwise_tukeyhsd` takes the values first and the group labels second. Swapping them produces a cryptic error or nonsensical comparisons, so double-check the argument order if the group names appear in the value column. ### Step 7: Handle unequal variances and non-normal data ![Decision flow from Shapiro-Wilk and Levene results to one-way ANOVA, Welch ANOVA, Kruskal-Wallis, and their post-hoc tests.](/article-images/illustrations/python-anova-bioinformatics-assumption-branches.webp) Two fallbacks cover the situations where the Step 3 checks fail. If variances are unequal but data are roughly normal, use Welch's ANOVA, which does not require equal variances: ```python f, p = stats.f_oneway(*samples) # classical w, p = stats.welch_anova if False else stats.f_oneway # placeholder removed ``` In current SciPy the call is simply: ```python from scipy import stats w, pw = stats.alexandergovern(*samples) # robust alternative in recent SciPy f_welch, p_welch = stats.f_oneway(*samples, equal_var=False) if hasattr(stats.f_oneway, "__defaults__") else (None, None) ``` That last line is unnecessary complexity. The clean approach is Pingouin, which has a dedicated Welch ANOVA: ```python import pingouin as pg aov = pg.anova(data=df, dv="abundance", between="group", detailed=True) print(aov) welch = pg.welch_anova(data=df, dv="abundance", between="group") print(welch) ``` Pingouin also returns effect sizes directly in `pg.anova`, which saves the manual computation from Step 5. If the data are clearly non-normal or ordinal, switch to the Kruskal-Wallis test, the nonparametric analogue of one-way ANOVA. It is a standard companion to ANOVA in metabolomics toolkits [3] and is one of the non-parametric options offered for differential methylation analysis [2]: ```python h, p = stats.kruskal(*samples) print(f"H = {h:.3f}, p = {p:.3g}") ``` Follow a significant Kruskal-Wallis result with pairwise Mann-Whitney U tests and a multiple-testing correction, for example using `stats.mannwhitneyu` on each pair and then `multipletests` as shown in Step 9. ### Step 8: Scale to thousands of features Real omics data are wide: thousands of proteins, transcripts, CpG sites, or metabolites, each of which needs its own test. AlphaPeptStats was built precisely for this, providing automated, scalable ANOVA and related analyses over label-free proteomics tables imported from popular search engines [1]. You can reproduce the core loop yourself: ```python # df_wide: rows are samples, columns are protein IDs # meta: DataFrame indexed like df_wide with a "group" column from statsmodels.stats.multitest import multipletests group_levels = ["control", "low_dose", "high_dose"] records = [] for feature in df_wide.columns: arrays = [ df_wide.loc[meta.index[meta["group"] == g], feature].dropna() for g in group_levels ] if min(len(a) for a in arrays) < 3: continue # skip features with too few observations f, p = stats.f_oneway(*arrays) records.append({"feature": feature, "F": f, "p": p}) res = pd.DataFrame(records).dropna() res["q_value"] = multipletests(res["p"], method="fdr_bh")[1] hits = res[res["q_value"] < 0.05].sort_values("q_value") print(hits.head(10)) ``` The Benjamini-Hochberg false discovery rate (FDR) procedure converts each raw p value into a q value, the estimated fraction of false positives among features at that threshold. The metabolomics study cited above applied FDR correction to ANOVA results across metabolite features for exactly this reason [5]. At 20,000 features and an uncorrected 0.05 threshold you would expect roughly 1,000 significant hits by chance alone. FDR correction is not optional at omics scale, it is the analysis. For very large tables, vectorize the loop with `scipy.stats.f_oneway` applied to arrays grouped with `pandas.groupby`, or use a domain package that handles this internally. ### Step 9: Two-way ANOVA and factorial designs When you have two factors, for example treatment and genotype, use a two-way ANOVA through statsmodels: ```python model2 = ols("abundance ~ C(treatment) * C(genotype)", data=df2).fit() print(sm.stats.anova_lm(model2, typ=2)) ``` The `*` operator includes both main effects and their interaction. A significant interaction means the treatment effect depends on genotype, which changes how you should interpret the main effects. Report each factor's F statistic, p value, and effect size separately. ## Choosing the Right Test: Quick Reference | Situation | Test | Python function | Key assumption | |--|--|--|--| | 3+ groups, normal data, equal variances | One-way ANOVA | `scipy.stats.f_oneway` or `statsmodels.ols` + `anova_lm` | Normal residuals, equal variances | | 3+ groups, unequal variances | Welch ANOVA | `pingouin.welch_anova` | Normal residuals only | | 3+ groups, non-normal or ordinal data | Kruskal-Wallis | `scipy.stats.kruskal` | Independent observations, similar distribution shapes | | Which pairs differ after ANOVA | Tukey HSD | `statsmodels pairwise_tukeyhsd` | Same as ANOVA | | Which pairs differ after Welch ANOVA | Games-Howell | `pingouin.gameshowell` | No equal-variance requirement | | Two factors plus interaction | Two-way ANOVA | `ols("y ~ C(a) * C(b)")` + `anova_lm` | Same as one-way, per residual | | Thousands of features at once | ANOVA loop + FDR | `multipletests(method="fdr_bh")` | Per-feature assumptions as above | ## How to Verify Your ANOVA Result Trustworthy results come from redundancy and sanity checks, not from running the test once and moving on. First, cross-check implementations. Run the test in both SciPy and statsmodels and confirm the F statistic and p value agree to several decimal places. Disagreement almost always means the two calls received different data, for example one dropped missing values and the other did not. Second, verify against an independent implementation when stakes are high. A validation study of AI-assisted data analysis in healthcare research checked one-way ANOVA outputs, among other tests, against expected statistical results and found that accuracy improved markedly when prompts and inputs were made explicit and structured [7]. The lesson transfers directly to scripted analysis: unambiguous inputs and explicit expectations make errors visible. Third, visualize. The boxplot from Step 2 should show the pattern your ANOVA claims. If the test says three groups differ but the boxes overlap almost completely with tiny samples, be suspicious. Fourth, check the direction and magnitude of post-hoc differences against what the descriptives show. Tukey output mean differences should match group means from Step 2. ## Domain-Specific Python Packages That Include ANOVA Before writing your own loop, check whether a domain package already implements the workflow with tested defaults. AlphaPeptStats is a Python package for statistical analysis of mass spectrometry proteomics. It handles normalization, imputation, and visualization, runs ANOVA and related tests in an automatable manner, imports output from a range of popular search engines, and ships with a testing framework covering 98 percent of its code [1]. If your data come from a proteomics pipeline, it is the shortest path from raw results table to publication-ready statistics. CpGtools provides modules for DNA methylation analysis from bisulfite sequencing or Illumina arrays, including differential CpG analysis with ANOVA, t tests, Fisher's exact test, linear regression, and non-parametric methods, plus PCA and visualization modules [2]. It is the natural choice when the groups you compare are methylation states at CpG sites. SECIMTools offers Python applications for metabolomics, standalone or wrapped in Galaxy, including ANOVA, t tests, and the Kruskal-Wallis test alongside quality control metrics and classification tools [3]. Its Galaxy integration suits labs that want reproducible workflows without a command line. PYCHEM was an early open-source multivariate analysis toolbox for Python built on SciPy, aimed at biological sciences, and a forerunner of the current Python statistics ecosystem [8]. It is largely historical now, but it documents how long Python has served this role. GeneQuantify applies Shapiro-Wilk and Levene testing to ΔCt values within a qPCR workflow built on Python and Streamlit [6]. Its pipeline is a good template for embedding assumption checks automatically rather than relying on the analyst to remember them. The Spatial Omics Toolkit (SPOT) implements Python and R pipelines for mass spectrometry imaging, integrating classification, dimensionality reduction, and comparative analysis across disease states [9]. It illustrates the current trend of packaging validated statistical workflows so that domain scientists apply them consistently. ## Common Mistakes and Limitations These are the errors that most often invalidate or distort ANOVA results in bioinformatics practice. **Skipping the multiple-testing correction.** The single most common mistake in omics. A p value below 0.05 for one protein out of 10,000 means nothing without FDR correction. Apply Benjamini-Hochberg before calling hits. **Stopping at the omnibus test.** A significant F statistic says the groups are not all equal. It does not say which groups differ or how. Always follow up with Tukey HSD or an equivalent post-hoc test and report the mean differences [4]. **Treating technical replicates as biological replicates.** Three pipetting repeats from one mouse are not three mice. ANOVA assumes independent observations. Pooling technical replicates or modeling them properly is required, otherwise the effective sample size is inflated and p values become meaningless. **Ignoring heteroscedasticity.** Omics measurements frequently show variance that grows with the mean. Log2 transformation before testing is standard for abundance data, and when variances remain unequal, Welch ANOVA or permutation approaches are safer than classical ANOVA. **Applying classical ANOVA to bounded or compositional data.** Methylation percentages are bounded between 0 and 100 and often pile up near the extremes. CpGtools explicitly offers non-parametric tests and regression alternatives alongside ANOVA for this reason [2]. Match the test to the data type. **Confusing significance with effect size.** With huge n, trivial differences reach significance. With tiny n, biologically large differences may not. Report η² or Cohen's d next to every p value. **Running ANOVA on raw RNA-seq counts.** Count data follow a negative binomial-like distribution, not a normal one. Use dedicated tools such as PyDESeq2, the Python implementation of DESeq2, which was used for differential expression analysis in a glioblastoma study comparing tumor and normal brain tissue [10]. ANOVA belongs on continuous, roughly normal measurements. **Checking assumptions only in the exemplar feature.** Normality and variance homogeneity should be assessed across features or at least on a representative sample, not on the single feature you plan to highlight. **Overwriting results without a record.** Interactive notebooks invite reruns with slightly different filters until the answer looks right. Fix the seed, log the parameters, and keep the first confirmed result. Also note that a one-way ANOVA on a single feature is a hypothesis about that feature, and individual experimental cases may require domain-specific judgment, including veterinary or clinical input when the data concern patient or animal samples. ## Beyond the Single Feature Once the per-feature ANOVA loop is solid, the natural extensions are multivariate. Principal component analysis and hierarchical clustering put the group structure in context before univariate testing. AlphaPeptStats bundles both alongside its ANOVA machinery [1], SECIMTools provides them for metabolomics [3], and SPOT brings dimensionality reduction to spatial molecular data [9]. For time-course experiments where repeated measures violate independence, bootstrapping approaches such as BooteJTK, implemented in Python for circadian rhythm detection, offer an alternative to repeated-measures ANOVA on sparsely sampled series [11]. ## Frequently Asked Questions ### When should I use ANOVA instead of multiple t tests? Use ANOVA whenever you compare three or more groups on one continuous measurement. Running every pair of groups as a separate t test inflates the false positive rate because each test adds its own chance of a Type I error. A single ANOVA controls the overall error rate at your chosen alpha, and post-hoc tests then handle the pairwise questions. ### Which Python library is best for ANOVA? For a quick one-way test, `scipy.stats.f_oneway` is sufficient. For ANOVA tables, effect sizes, multi-factor designs, and post-hoc tests, use statsmodels. Pingouin offers the friendliest API with built-in effect sizes and Welch ANOVA. For proteomics or methylation data, domain packages such as AlphaPeptStats or CpGtools wrap these methods in pipelines designed for your data format [1] [2]. ### How do I check ANOVA assumptions in Python? Run Shapiro-Wilk on each group with `scipy.stats.shapiro` for normality and Levene's test with `scipy.stats.levene(center="median")` for equal variances, then inspect a boxplot. This pairing is standard in validated qPCR workflows implemented in Python [6]. ### What do I do if Levene's test fails? Switch to Welch's ANOVA, available as `pingouin.welch_anova`, which does not assume equal variances, and follow it with the Games-Howell post-hoc test. If the data are also non-normal, use the Kruskal-Wallis test with pairwise Mann-Whitney comparisons and FDR correction. ### How do I run ANOVA on thousands of genes or proteins at once? Loop or vectorize `f_oneway` over the feature columns of your wide-format table, collect the p values, and pass them to `statsmodels.stats.multitest.multipletests` with the Benjamini-Hochberg method. Analyze on a per-feature basis with q values as your hit criterion. ### What is the difference between one-way and two-way ANOVA? One-way ANOVA tests one grouping factor, for example three treatments. Two-way ANOVA tests two factors at once, for example treatment and genotype, including whether their effects interact. Two-way ANOVA in Python uses the statsmodels formula interface with `C(factor1) * C(factor2)`. ### Do I need a post-hoc test after ANOVA? Yes, whenever the omnibus test is significant and you want to name the groups that differ. Tukey's HSD is the standard choice for classical ANOVA because it controls the family-wise error rate across all pairwise comparisons, and it is the follow-up used in published Python-analyzed studies [4]. ### Can I use ANOVA on RNA-seq count data? Not directly. Raw counts are non-negative integers with mean-dependent variance, which violates ANOVA assumptions. Use a negative binomial framework such as PyDESeq2, the Python implementation of DESeq2 used for differential expression analysis in published transcriptomic studies [10]. ANOVA is appropriate for normalized continuous data such as log2 protein abundances or ΔCt values [6]. ## Sources 1. [AlphaPeptStats: an open-source Python package for automated and scalable statistical analysis of mass spectrometry-based proteomics.](https://pubmed.ncbi.nlm.nih.gov/37527012/) 2. [CpGtools: a python package for DNA methylation analysis.](https://pubmed.ncbi.nlm.nih.gov/31808791/) 3. [SECIMTools: a suite of metabolomics data analysis tools.](https://pubmed.ncbi.nlm.nih.gov/29678131/) 4. [Integrating Computational Analysis of In Vivo Investigation of Modulatory Effect of Fagonia cretica Plant Extract on Letrozole-Induced Polycystic Ovary Syndrome in Female Rats.](https://pubmed.ncbi.nlm.nih.gov/40723460/) 5. [Mass Spectrometry-Based Untargeted Metabolomics Identifies Distinct Metabolic Signatures in Infertility: A Comparative Analysis of PCOS, POR, and NOR.](https://pubmed.ncbi.nlm.nih.gov/40537734/) 6. [GeneQuantify: a web-based tool for qPCR gene expression and copy number variation analysis.](https://pubmed.ncbi.nlm.nih.gov/42384344/) 7. [ChatGPT for Univariate Statistics: Validation of AI-Assisted Data Analysis in Healthcare Research.](https://pubmed.ncbi.nlm.nih.gov/39919289/) 8. [PYCHEM: a multivariate analysis package for python.](https://pubmed.ncbi.nlm.nih.gov/16882648/) 9. [A Protocol for Multivariate Data Visualization and Pseudotime Modeling for Analysis of Disease Trajectories Detected by Mass Spectrometry Imaging.](https://pubmed.ncbi.nlm.nih.gov/42683911/) 10. [Comprehensive Differential Gene Expression Analysis in Glioblastoma Using PyDESeq2: A Comparison with Normal Brain Tissue.](https://pubmed.ncbi.nlm.nih.gov/41479088/) 11. [Bootstrapping and Empirical Bayes Methods Improve Rhythm Detection in Sparsely Sampled Data.](https://pubmed.ncbi.nlm.nih.gov/30101659/)