diff --git a/examples/plot_minimal_pydeseq2_pipeline.py b/examples/plot_minimal_pydeseq2_pipeline.py
index b3a7431f..8141bf08 100644
--- a/examples/plot_minimal_pydeseq2_pipeline.py
+++ b/examples/plot_minimal_pydeseq2_pipeline.py
@@ -381,6 +381,31 @@
 ds_Y_vs_X = DeseqStats(dds, contrast=["group", "Y", "X"], inference=inference)
 ds_Y_vs_X.summary()
 
+# %%
+# .. _lrt_ref:
+#
+# Likelihood ratio test
+# """""""""""""""""""""
+#
+# Instead of the Wald test, we may assess significance with a likelihood ratio test
+# (LRT), as in R DESeq2's ``DESeq(dds, test="LRT", reduced=...)``. The LRT compares the
+# full model to a nested ``reduced`` model that drops the term(s) being tested, and is
+# especially useful to test a factor with several levels jointly, or to compare nested
+# models in general. Here we test the effect of ``condition`` on top of ``group`` by
+# comparing the full ``~group + condition`` model to a reduced ``~group`` model.
+#
+# The p-value reflects the model comparison, while the reported ``log2FoldChange`` still
+# corresponds to the requested ``contrast`` of the full model.
+
+ds_lrt = DeseqStats(
+    dds,
+    contrast=["condition", "B", "A"],
+    test="LRT",
+    reduced="~group",
+    inference=inference,
+)
+ds_lrt.summary()
+
 # %%
 # LFC shrinkage (multifactor)
 # """""""""""""""""""""""""""
diff --git a/pydeseq2/ds.py b/pydeseq2/ds.py
index 6e22fde7..dd9bc5ef 100644
--- a/pydeseq2/ds.py
+++ b/pydeseq2/ds.py
@@ -2,11 +2,14 @@
 import time
 import warnings
 from typing import Literal
+from typing import cast
 
 # import anndata as ad
 import numpy as np
 import pandas as pd
+from formulaic_contrasts import FormulaicContrasts
 from scipy.optimize import root_scalar  # type: ignore
+from scipy.stats import chi2  # type: ignore
 from scipy.stats import false_discovery_control  # type: ignore
 
 from pydeseq2.dds import DeseqDataSet
@@ -14,6 +17,7 @@
 from pydeseq2.inference import Inference
 from pydeseq2.utils import lowess
 from pydeseq2.utils import make_MA_plot
+from pydeseq2.utils import nb_nll
 
 
 class DeseqStats:
@@ -29,6 +33,21 @@ class DeseqStats:
     dds : DeseqDataSet
         DeseqDataSet for which dispersion and LFCs were already estimated.
 
+    test : str
+        The statistical test to use for p-value estimation. One of ``["Wald", "LRT"]``.
+        The Wald test assesses the significance of a single coefficient (or contrast),
+        while the likelihood ratio test (``"LRT"``) compares the full model to a nested
+        ``reduced`` model, as in R DESeq2's ``DESeq(dds, test="LRT", reduced=...)``.
+        (default: ``"Wald"``).
+
+    reduced : str, ndarray, pandas.DataFrame, optional
+        The reduced model to compare against the full model when ``test="LRT"``.
+        Either a formulaic formula string (e.g. ``"~group"``, which must be nested in the
+        design of ``dds``), or an explicit design matrix (as a numpy array or a pandas
+        DataFrame whose columns are a subset of the full design matrix columns).
+        Required when ``test="LRT"``, and must be left to ``None`` for the Wald test.
+        (default: ``None``).
+
     contrast : list or ndarray
         Either a list of three strings or a numpy array.
         If a list of three strings, it must be in the following format:
@@ -102,10 +121,19 @@ class DeseqStats:
         Standard LFC error.
 
     statistics : pandas.Series
-        Wald statistics.
+        Wald statistics (``test="Wald"``) or likelihood ratio statistics
+        (``test="LRT"``).
 
     p_values : pandas.Series
-        P-values estimated from Wald statistics.
+        P-values estimated from the Wald statistics (``test="Wald"``) or from the
+        likelihood ratio statistics using a chi-squared distribution (``test="LRT"``).
+
+    test : str
+        The statistical test used for p-value estimation, one of ``["Wald", "LRT"]``.
+
+    reduced_design_matrix : pandas.DataFrame, optional
+        The reduced-model design matrix used for the likelihood ratio test
+        (``None`` for the Wald test).
 
     padj : pandas.Series
         P-values adjusted for multiple testing.
@@ -132,6 +160,8 @@ def __init__(
         self,
         dds: DeseqDataSet,
         contrast: list[str] | np.ndarray,
+        test: Literal["Wald", "LRT"] = "Wald",
+        reduced: str | np.ndarray | pd.DataFrame | None = None,
         alpha: float = 0.05,
         cooks_filter: bool = True,
         independent_filter: bool = True,
@@ -150,6 +180,10 @@ def __init__(
 
         self.dds = dds
 
+        if test not in ("Wald", "LRT"):
+            raise ValueError(f"test must be one of 'Wald' or 'LRT', got '{test}'.")
+        self.test = test
+
         self.alpha = alpha
         self.cooks_filter = cooks_filter
         self.independent_filter = independent_filter
@@ -166,8 +200,11 @@ def __init__(
 
         # Initialize the design matrix and LFCs. If the chosen reference level are the
         # same as in dds, keep them unchanged. Otherwise, change reference level.
-        self.design_matrix = self.dds.obsm["design_matrix"].copy()
-        self.LFC = self.dds.varm["LFC"].copy()
+        self.design_matrix = cast(pd.DataFrame, self.dds.obsm["design_matrix"].copy())
+        self.LFC = cast(pd.DataFrame, self.dds.varm["LFC"].copy())
+
+        # Build the reduced-model design matrix for the likelihood ratio test.
+        self.reduced_design_matrix = self._build_reduced_design_matrix(reduced)
 
         # Check the validity of the contrast (if provided) or build it.
         self.contrast: list[str] | np.ndarray
@@ -252,16 +289,29 @@ def summary(
                 f"positive lfc_null value (got {lfc_null}).",
             )
 
+        if self.test == "LRT" and (
+            new_lfc_null != "default" or new_alt_hypothesis != "default"
+        ):
+            warnings.warn(
+                "`lfc_null` and `alt_hypothesis` are only supported for the Wald test "
+                "and are ignored when `test='LRT'`.",
+                UserWarning,
+                stacklevel=2,
+            )
+
         if (
             not hasattr(self, "p_values")
             or self.lfc_null != lfc_null
             or self.alt_hypothesis != alt_hypothesis
         ):
-            # Estimate p-values with Wald test
+            # Estimate p-values with the Wald test or the likelihood ratio test
             self.lfc_null = lfc_null
             self.alt_hypothesis = alt_hypothesis
             rerun_summary = True
-            self.run_wald_test()
+            if self.test == "Wald":
+                self.run_wald_test()
+            else:
+                self.run_likelihood_ratio_test()
 
         if self.cooks_filter:
             # Filter p-values based on Cooks outliers
@@ -286,16 +336,20 @@ def summary(
         self.results_df["padj"] = self.padj
 
         if not self.quiet:
+            pval_desc = (
+                "Wald test p-value"
+                if self.test == "Wald"
+                else "likelihood ratio test p-value"
+            )
             if isinstance(self.contrast, np.ndarray):
                 # The contrast vector was directly provided
                 print(
-                    f"Log2 fold change & Wald test p-value, contrast vector: "
-                    f"{self.contrast}"
+                    f"Log2 fold change & {pval_desc}, contrast vector: {self.contrast}"
                 )
             else:
                 # The factor is categorical
                 print(
-                    f"Log2 fold change & Wald test p-value: "
+                    f"Log2 fold change & {pval_desc}: "
                     f"{self.contrast[0]} {self.contrast[1]} vs {self.contrast[2]}"
                 )
             print(self.results_df)
@@ -359,6 +413,123 @@ def run_wald_test(self) -> None:
             self.statistics.loc[self.dds.new_all_zeroes_genes] = 0.0
             self.p_values.loc[self.dds.new_all_zeroes_genes] = 1.0
 
+    def run_likelihood_ratio_test(self) -> None:
+        r"""Perform a likelihood ratio test (LRT).
+
+        Compares the full model (the design of the ``DeseqDataSet``) to a nested
+        ``reduced`` model, using the same gene-wise dispersions for both. This is the
+        Python equivalent of R DESeq2's ``DESeq(dds, test="LRT", reduced=...)``.
+
+        For each gene, the test statistic is
+
+        .. math::
+            \Lambda = 2 \left( \ell_{\text{full}} - \ell_{\text{reduced}} \right),
+
+        where :math:`\ell` is the negative-binomial log-likelihood evaluated at the
+        maximum-likelihood fit of each model. Under the null hypothesis that the extra
+        full-model coefficients are zero, :math:`\Lambda` follows a chi-squared
+        distribution with degrees of freedom equal to the difference in the number of
+        coefficients between the two models. The reported ``log2FoldChange`` and
+        ``lfcSE`` still correspond to the requested ``contrast`` of the full model,
+        exactly as in R DESeq2.
+        """
+        if self.shrunk_LFCs and not self.quiet:
+            print(
+                "Note: running the likelihood ratio test on shrunk LFCs. The LRT is "
+                "usually run on the maximum-likelihood (unshrunk) estimates.",
+                file=sys.stderr,
+            )
+
+        # A reduced design matrix is always set when test="LRT".
+        assert self.reduced_design_matrix is not None
+
+        non_zero_idx = self.dds.non_zero_idx
+        non_zero_genes = self.dds.non_zero_genes
+
+        full_design_matrix = self.design_matrix.values
+        reduced_design_matrix = self.reduced_design_matrix.values
+        df = full_design_matrix.shape[1] - reduced_design_matrix.shape[1]
+
+        size_factors = np.asarray(self.dds.obs["size_factors"], dtype=float)
+        disp = np.asarray(self.dds.var["dispersions"], dtype=float)[non_zero_idx]
+
+        # Counts the full model was fit on: original counts, with Cooks outliers
+        # replaced by imputed values for refitted genes (as in R's
+        # ``counts(dds, replaced=TRUE)``). This keeps the reduced-model fit and the
+        # deviances consistent with the (possibly refitted) full-model LFCs.
+        # ``np.array`` (not ``asarray``) to guarantee a copy, so the in-place
+        # replacement below never mutates ``dds.X``.
+        counts = np.array(self.dds.X, dtype=float)
+        if self.dds.refit_cooks and hasattr(self.dds, "counts_to_refit"):
+            refit_pos = self.dds.var_names.get_indexer(
+                self.dds.counts_to_refit.var_names
+            )
+            counts[:, refit_pos] = np.asarray(self.dds.counts_to_refit.X, dtype=float)
+        counts = counts[:, non_zero_idx]
+
+        # Full-model mean, recomputed from the stored LFCs (natural log scale),
+        # matching the convention used by the Wald test.
+        lfc = self.LFC.values[non_zero_idx]
+        mu_full = np.exp(full_design_matrix @ lfc.T) * size_factors[:, None]
+
+        # Set regularization factors (identical to the Wald test).
+        if self.prior_LFC_var is not None:
+            ridge_factor = np.diag(1 / self.prior_LFC_var**2)
+        else:
+            ridge_factor = np.diag(np.repeat(1e-6, full_design_matrix.shape[1]))
+
+        if not self.quiet:
+            print("Running LRT tests...", file=sys.stderr)
+        start = time.time()
+
+        # Fit the reduced model with the SAME dispersions as the full model.
+        _, mu_reduced, _, _ = self.inference.irls(
+            counts=counts,
+            size_factors=size_factors,
+            design_matrix=reduced_design_matrix,
+            disp=disp,
+            min_mu=self.dds.min_mu,
+            beta_tol=self.dds.beta_tol,
+        )
+
+        # Standard error of the requested contrast, from the full-model fit, so that
+        # the ``lfcSE`` column matches the Wald test output.
+        _, _, se = self.inference.wald_test(
+            design_matrix=full_design_matrix,
+            disp=disp,
+            lfc=lfc,
+            mu=mu_full,
+            ridge_factor=ridge_factor,
+            contrast=self.contrast_vector,
+            lfc_null=np.array(0.0),
+            alt_hypothesis=None,
+        )
+
+        # LRT statistic and chi-squared p-value. ``nb_nll`` is the negative
+        # log-likelihood, so 2 * (ll_full - ll_reduced) = 2 * (nll_reduced - nll_full).
+        stats = 2.0 * (nb_nll(counts, mu_reduced, disp) - nb_nll(counts, mu_full, disp))
+        # Clip tiny negative values that can arise from numerical noise when the two
+        # models fit essentially identically (the true statistic is non-negative).
+        stats = np.maximum(stats, 0.0)
+        pvals = chi2.sf(stats, df)
+
+        end = time.time()
+        if not self.quiet:
+            print(f"... done in {end - start:.2f} seconds.\n", file=sys.stderr)
+
+        self.p_values = pd.Series(np.nan, index=self.dds.var_names)
+        self.statistics = pd.Series(np.nan, index=self.dds.var_names)
+        self.SE = pd.Series(np.nan, index=self.dds.var_names)
+        self.p_values.loc[non_zero_genes] = pvals
+        self.statistics.loc[non_zero_genes] = stats
+        self.SE.loc[non_zero_genes] = se
+
+        # Account for possible all_zeroes due to outlier refitting in DESeqDataSet
+        if self.dds.refit_cooks and self.dds.var["replaced"].sum() > 0:
+            self.SE.loc[self.dds.new_all_zeroes_genes] = 0.0
+            self.statistics.loc[self.dds.new_all_zeroes_genes] = 0.0
+            self.p_values.loc[self.dds.new_all_zeroes_genes] = 1.0
+
     # TODO update this to reflect the new contrast format
     def lfc_shrink(self, coeff: str, adapt: bool = True) -> None:
         """LFC shrinkage with an apeGLM prior :cite:p:`DeseqStats-zhu2019heavy`.
@@ -494,7 +665,7 @@ def _independent_filtering(self) -> None:
         """
         # Check that p-values are available. If not, compute them.
         if not hasattr(self, "p_values"):
-            self.run_wald_test()
+            self._run_test()
 
         lower_quantile = np.mean(self.base_mean == 0)
 
@@ -537,8 +708,8 @@ def _p_value_adjustment(self) -> None:
         This method and the `_independent_filtering` are mutually exclusive.
         """
         if not hasattr(self, "p_values"):
-            # Estimate p-values with Wald test
-            self.run_wald_test()
+            # Estimate p-values with the configured test
+            self._run_test()
 
         self.padj = pd.Series(np.nan, index=self.dds.var_names)
         self.padj.loc[~self.p_values.isna()] = false_discovery_control(
@@ -549,7 +720,7 @@ def _cooks_filtering(self) -> None:
         """Filter p-values based on Cooks outliers."""
         # Check that p-values are available. If not, compute them.
         if not hasattr(self, "p_values"):
-            self.run_wald_test()
+            self._run_test()
 
         self.p_values[self.dds.cooks_outlier()] = np.nan
 
@@ -603,3 +774,92 @@ def _build_contrast_vector(self) -> None:
         self.contrast_vector = self.dds.contrast(
             column=factor, baseline=ref, group_to_compare=alternative
         )
+
+    def _run_test(self) -> None:
+        """Run the configured statistical test (Wald or LRT)."""
+        if self.test == "Wald":
+            self.run_wald_test()
+        else:
+            self.run_likelihood_ratio_test()
+
+    def _build_reduced_design_matrix(
+        self, reduced: str | np.ndarray | pd.DataFrame | None
+    ) -> pd.DataFrame | None:
+        """Validate and build the reduced-model design matrix for the LRT.
+
+        Parameters
+        ----------
+        reduced : str, ndarray, pandas.DataFrame, optional
+            The reduced model, as passed to the constructor.
+
+        Returns
+        -------
+        pandas.DataFrame or None
+            The reduced-model design matrix (``None`` for the Wald test).
+        """
+        if self.test == "Wald":
+            if reduced is not None:
+                raise ValueError(
+                    "`reduced` is only used for the likelihood ratio test "
+                    "(`test='LRT'`); leave it to None for the Wald test."
+                )
+            return None
+
+        # test == "LRT"
+        if reduced is None:
+            raise ValueError("A `reduced` model must be provided when `test='LRT'`.")
+
+        if isinstance(reduced, str):
+            if not isinstance(self.dds.design, str):
+                raise ValueError(
+                    "A formula string can only be used for `reduced` when the "
+                    "DeseqDataSet design is itself a formula. Provide `reduced` as a "
+                    "design matrix (numpy array or pandas DataFrame) instead."
+                )
+            reduced_design_matrix = FormulaicContrasts(
+                self.dds.obs, reduced
+            ).design_matrix
+            # The reduced model must be nested in the full model.
+            full_columns = set(self.design_matrix.columns)
+            if not set(reduced_design_matrix.columns).issubset(full_columns):
+                raise ValueError(
+                    "The reduced model must be nested in the full model: its columns "
+                    f"{list(reduced_design_matrix.columns)} must be a subset of the "
+                    f"full design matrix columns {list(self.design_matrix.columns)}."
+                )
+        elif isinstance(reduced, pd.DataFrame):
+            reduced_design_matrix = reduced.copy()
+        elif isinstance(reduced, np.ndarray):
+            if reduced.ndim != 2 or reduced.shape[0] != self.dds.n_obs:
+                raise ValueError(
+                    "The reduced design matrix must be 2D with one row per sample "
+                    f"({self.dds.n_obs}); got shape {reduced.shape}."
+                )
+            reduced_design_matrix = pd.DataFrame(
+                reduced,
+                index=self.design_matrix.index,
+                columns=[f"reduced_{i}" for i in range(reduced.shape[1])],
+            )
+        else:
+            raise TypeError(
+                "`reduced` must be a formula string, a numpy array, or a pandas "
+                f"DataFrame; got {type(reduced).__name__}."
+            )
+
+        # The reduced model must have strictly fewer coefficients than the full model.
+        n_df = self.design_matrix.shape[1] - reduced_design_matrix.shape[1]
+        if n_df < 1:
+            raise ValueError(
+                "The reduced model must have fewer coefficients than the full model "
+                f"(full: {self.design_matrix.shape[1]}, "
+                f"reduced: {reduced_design_matrix.shape[1]}). "
+                "The likelihood ratio test compares nested models."
+            )
+        if reduced_design_matrix.shape[0] != self.design_matrix.shape[0]:
+            raise ValueError(
+                "The reduced design matrix must have one row per sample "
+                f"({self.design_matrix.shape[0]}); got "
+                f"{reduced_design_matrix.shape[0]}."
+            )
+
+        return reduced_design_matrix
diff --git a/tests/data/README.md b/tests/data/README.md
index b496a61c..1b04902d 100644
--- a/tests/data/README.md
+++ b/tests/data/README.md
@@ -12,3 +12,18 @@ in `/datasets/synthetic/`, respectively using `~condition` and `~condition + gro
 - `r_test_size_factors.csv` contains DESeq2's `estimateSizeFactors` output,
 - `r_vst.csv` contains DESeq2's `varianceStabilizingTransformation` output with `blind=TRUE` and `fitType="parametric"`,
 - `r_vst_with_design.csv` contains DESeq2's `varianceStabilizingTransformation` output with `blind=FALSE` and `fitType="parametric"`.
+
+### Likelihood ratio test (LRT) reference data
+
+The `r_test_res_lrt*.csv` files contain the output of DESeq2's `results()` after running
+`DESeq(dds, test="LRT", reduced=...)`, and are used by the LRT tests. They were generated
+with DESeq2 1.46.0 by `generate_lrt_reference.R` (run `Rscript tests/data/generate_lrt_reference.R .`
+from the repository root to regenerate them). For each case, a `*_no_independent_filtering.csv`
+variant is also stored (i.e. `results(..., independentFiltering=FALSE)`), which isolates the raw
+LRT statistic and p-value from the p-value adjustment step.
+
+- `single_factor/r_test_res_lrt.csv`: full `~condition` vs reduced `~1` (`df=1`),
+- `multi_factor/r_test_res_lrt.csv`: full `~group + condition` vs reduced `~group` (`df=1`),
+- `multi_factor/r_test_res_lrt_reduced_intercept.csv`: full `~group + condition` vs reduced `~1` (`df=2`),
+- `multi_factor/r_test_res_lrt_outliers.csv`: same as above (`~group + condition` vs `~group`) but with
+  injected Cooks outliers that trigger outlier replacement and model refitting.
diff --git a/tests/data/generate_lrt_reference.R b/tests/data/generate_lrt_reference.R
new file mode 100644
index 00000000..86b53171
--- /dev/null
+++ b/tests/data/generate_lrt_reference.R
@@ -0,0 +1,66 @@
+## Generate the R DESeq2 reference outputs used by the PyDESeq2 likelihood ratio
+## test (LRT) suite.
+##
+## These files are the ground truth that `tests/test_pydeseq2.py` compares the
+## PyDESeq2 LRT implementation against, mirroring how the existing Wald-test
+## reference data was produced. They use DESeq2 defaults (parametric dispersion
+## trend, median-of-ratios size factors, Cooks refitting), matching the defaults
+## of `DeseqDataSet`/`DeseqStats`.
+##
+## Usage (from the repository root):
+##     Rscript tests/data/generate_lrt_reference.R .
+##
+## Requires: DESeq2 (reference data was generated with DESeq2 1.46.0).
+
+suppressMessages(library(DESeq2))
+
+repo <- commandArgs(trailingOnly = TRUE)[1]
+if (is.na(repo)) repo <- "."
+
+counts   <- read.csv(file.path(repo, "datasets/synthetic/test_counts.csv"),   row.names = 1)  # genes x samples
+metadata <- read.csv(file.path(repo, "datasets/synthetic/test_metadata.csv"), row.names = 1)  # samples x factors
+
+## Run an LRT and write both the default `results()` output (independent
+## filtering ON) and a version with independent filtering OFF (which isolates the
+## raw LRT statistic and p-value from the p-value adjustment step).
+run_lrt <- function(counts, metadata, full, reduced, name, outfile) {
+  dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = full)
+  dds <- DESeq(dds, test = "LRT", reduced = reduced, quiet = TRUE)
+  write.csv(as.data.frame(results(dds, name = name)), outfile)
+  write.csv(
+    as.data.frame(results(dds, name = name, independentFiltering = FALSE)),
+    sub("\\.csv$", "_no_independent_filtering.csv", outfile)
+  )
+  invisible(dds)
+}
+
+## --- Clean data (no outliers) -------------------------------------------------
+meta <- metadata
+meta$condition <- factor(meta$condition, levels = c("A", "B"))
+meta$group     <- factor(meta$group,     levels = c("X", "Y"))
+
+# Case A: single factor, full = ~condition, reduced = ~1 (df = 1)
+run_lrt(counts, meta, ~condition, ~1, "condition_B_vs_A",
+        file.path(repo, "tests/data/single_factor/r_test_res_lrt.csv"))
+
+# Case B: multi factor, full = ~group + condition, reduced = ~group (df = 1)
+run_lrt(counts, meta, ~group + condition, ~group, "condition_B_vs_A",
+        file.path(repo, "tests/data/multi_factor/r_test_res_lrt.csv"))
+
+# Case C: multi factor, full = ~group + condition, reduced = ~1 (df = 2)
+run_lrt(counts, meta, ~group + condition, ~1, "condition_B_vs_A",
+        file.path(repo, "tests/data/multi_factor/r_test_res_lrt_reduced_intercept.csv"))
+
+## --- Data with Cooks outliers (triggers refitting) ---------------------------
+## Same mild outlier injection as the existing Wald outlier test.
+counts_out <- counts
+counts_out["gene1", "sample1"]  <- 2000L
+counts_out["gene7", "sample11"] <- 1000L
+
+dds_out <- DESeqDataSetFromMatrix(countData = counts_out, colData = meta,
+                                  design = ~group + condition)
+dds_out <- DESeq(dds_out, test = "LRT", reduced = ~group, quiet = TRUE)
+write.csv(as.data.frame(results(dds_out, name = "condition_B_vs_A")),
+          file.path(repo, "tests/data/multi_factor/r_test_res_lrt_outliers.csv"))
+
+cat("Done. n replaced (outlier case):", sum(mcols(dds_out)$replace, na.rm = TRUE), "\n")
diff --git a/tests/data/multi_factor/r_test_res_lrt.csv b/tests/data/multi_factor/r_test_res_lrt.csv
new file mode 100644
index 00000000..f5bafda4
--- /dev/null
+++ b/tests/data/multi_factor/r_test_res_lrt.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,6.31579426514429,0.0119667128198756,0.0299167820496891
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,12.6913815611018,0.000367344758633349,0.00122448252877783
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,5.38895826761268,0.0202645651531322,0.0405291303062643
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,15.8216712349465,6.96007574745871e-05,0.000619867084233537
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,14.7312636318204,0.000123973416846707,0.000619867084233537
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,0.0039664276938538,0.949782763590195,0.949782763590195
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,0.813353148799138,0.36713074461821,0.407923049575788
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,4.25158626412713,0.0392136864907029,0.0653561441511715
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,2.57942335791518,0.108261094491313,0.141449313668176
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,2.50953313686284,0.113159450934541,0.141449313668176
diff --git a/tests/data/multi_factor/r_test_res_lrt_no_independent_filtering.csv b/tests/data/multi_factor/r_test_res_lrt_no_independent_filtering.csv
new file mode 100644
index 00000000..f5bafda4
--- /dev/null
+++ b/tests/data/multi_factor/r_test_res_lrt_no_independent_filtering.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,6.31579426514429,0.0119667128198756,0.0299167820496891
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,12.6913815611018,0.000367344758633349,0.00122448252877783
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,5.38895826761268,0.0202645651531322,0.0405291303062643
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,15.8216712349465,6.96007574745871e-05,0.000619867084233537
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,14.7312636318204,0.000123973416846707,0.000619867084233537
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,0.0039664276938538,0.949782763590195,0.949782763590195
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,0.813353148799138,0.36713074461821,0.407923049575788
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,4.25158626412713,0.0392136864907029,0.0653561441511715
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,2.57942335791518,0.108261094491313,0.141449313668176
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,2.50953313686284,0.113159450934541,0.141449313668176
diff --git a/tests/data/multi_factor/r_test_res_lrt_outliers.csv b/tests/data/multi_factor/r_test_res_lrt_outliers.csv
new file mode 100644
index 00000000..c9237c93
--- /dev/null
+++ b/tests/data/multi_factor/r_test_res_lrt_outliers.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.55034480088269,0.760136902014561,0.285643269329525,6.8291022584383,0.00896843310671871,0.0224210827667968
+"gene2",21.3782230626323,0.53636303894464,0.152196255941874,12.3502395356451,0.000440929894170484,0.00146976631390161
+"gene3",5.02959691052961,-0.674061326257139,0.287293062831826,5.39886648506354,0.02014983406417,0.04029966812834
+"gene4",101.066398705215,-0.424994832954691,0.110699439617865,14.6717658138381,0.000127948581466862,0.000911728459446636
+"gene5",27.3066986731155,0.588386623184279,0.156762787257669,14.0047881763476,0.000182345691889327,0.000911728459446636
+"gene6",5.446988786503,-0.019685001214121,0.309645507173678,0.00399799057129258,0.949583622558906,0.949583622558906
+"gene7",28.4098070748334,0.13780930645167,0.152926774840609,0.809310627127047,0.36832414249836,0.409249047220401
+"gene8",40.5694778066866,-0.270864009796395,0.135240979433579,3.9986201979707,0.0455375283850031,0.0758958806416718
+"gene9",37.3644308033897,-0.212456964616852,0.136834847771096,2.40207465321987,0.121174453507509,0.151468066884386
+"gene10",11.6436906897954,0.386870211256403,0.245746587087819,2.45919232559606,0.116838846893463,0.151468066884386
diff --git a/tests/data/multi_factor/r_test_res_lrt_reduced_intercept.csv b/tests/data/multi_factor/r_test_res_lrt_reduced_intercept.csv
new file mode 100644
index 00000000..a469f519
--- /dev/null
+++ b/tests/data/multi_factor/r_test_res_lrt_reduced_intercept.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,11.650996821073,0.00295133281036829,0.00491888801728048
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,13.4478350080053,0.00120182084187058,0.00295217659418973
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,13.8294233735875,0.000993067731985363,0.00295217659418973
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,42.2876068754728,6.5669368585673e-10,6.5669368585673e-09
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,17.8105083122929,0.000135674197059866,0.000678370985299332
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,2.6977762507496,0.259528663617795,0.28836518179755
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,1.86785410620894,0.393007315373388,0.393007315373388
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,13.0367194654262,0.00147608829709487,0.00295217659418973
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,3.20467970692152,0.201424661980554,0.251780827475693
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,5.63505781320873,0.0597534166962626,0.0853620238518037
diff --git a/tests/data/multi_factor/r_test_res_lrt_reduced_intercept_no_independent_filtering.csv b/tests/data/multi_factor/r_test_res_lrt_reduced_intercept_no_independent_filtering.csv
new file mode 100644
index 00000000..a469f519
--- /dev/null
+++ b/tests/data/multi_factor/r_test_res_lrt_reduced_intercept_no_independent_filtering.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,11.650996821073,0.00295133281036829,0.00491888801728048
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,13.4478350080053,0.00120182084187058,0.00295217659418973
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,13.8294233735875,0.000993067731985363,0.00295217659418973
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,42.2876068754728,6.5669368585673e-10,6.5669368585673e-09
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,17.8105083122929,0.000135674197059866,0.000678370985299332
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,2.6977762507496,0.259528663617795,0.28836518179755
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,1.86785410620894,0.393007315373388,0.393007315373388
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,13.0367194654262,0.00147608829709487,0.00295217659418973
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,3.20467970692152,0.201424661980554,0.251780827475693
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,5.63505781320873,0.0597534166962626,0.0853620238518037
diff --git a/tests/data/single_factor/r_test_res_lrt.csv b/tests/data/single_factor/r_test_res_lrt.csv
new file mode 100644
index 00000000..a72b91a6
--- /dev/null
+++ b/tests/data/single_factor/r_test_res_lrt.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.632812315254828,0.289103638926867,4.73285162924537,0.0295917841865267,0.0663276489248662
+"gene2",21.2812387436366,0.538551973774898,0.149962947101771,12.8392275422336,0.000339427708341099,0.0016971385417055
+"gene3",5.01012348853472,-0.632831858568865,0.295221453480919,4.53730195229139,0.0331638244624331,0.0663276489248662
+"gene4",100.51796142035,-0.412101731455792,0.118627877254952,12.0212585890596,0.000525971731194201,0.00175323910398067
+"gene5",27.1424502740787,0.582066403700141,0.15473043602656,14.0806522268775,0.000175136067078282,0.0016971385417055
+"gene6",5.4130427476525,0.00145730049165907,0.310309910140032,0.00819390042352097,0.927873869713322,0.927873869713322
+"gene7",28.2940230404605,0.134335684024563,0.149986853732723,0.805396982860088,0.369484643110571,0.410538492345078
+"gene8",40.3583444203556,-0.270655780074565,0.136401659050961,3.92714778918253,0.0475124627494821,0.0791874379158035
+"gene9",37.1661826339853,-0.212714657425418,0.133243837041424,2.54476500804708,0.110660374187379,0.144411939198395
+"gene10",11.5893249023836,0.386011089345352,0.244586104807745,2.47690282625604,0.115529551358716,0.144411939198395
diff --git a/tests/data/single_factor/r_test_res_lrt_no_independent_filtering.csv b/tests/data/single_factor/r_test_res_lrt_no_independent_filtering.csv
new file mode 100644
index 00000000..a72b91a6
--- /dev/null
+++ b/tests/data/single_factor/r_test_res_lrt_no_independent_filtering.csv
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.632812315254828,0.289103638926867,4.73285162924537,0.0295917841865267,0.0663276489248662
+"gene2",21.2812387436366,0.538551973774898,0.149962947101771,12.8392275422336,0.000339427708341099,0.0016971385417055
+"gene3",5.01012348853472,-0.632831858568865,0.295221453480919,4.53730195229139,0.0331638244624331,0.0663276489248662
+"gene4",100.51796142035,-0.412101731455792,0.118627877254952,12.0212585890596,0.000525971731194201,0.00175323910398067
+"gene5",27.1424502740787,0.582066403700141,0.15473043602656,14.0806522268775,0.000175136067078282,0.0016971385417055
+"gene6",5.4130427476525,0.00145730049165907,0.310309910140032,0.00819390042352097,0.927873869713322,0.927873869713322
+"gene7",28.2940230404605,0.134335684024563,0.149986853732723,0.805396982860088,0.369484643110571,0.410538492345078
+"gene8",40.3583444203556,-0.270655780074565,0.136401659050961,3.92714778918253,0.0475124627494821,0.0791874379158035
+"gene9",37.1661826339853,-0.212714657425418,0.133243837041424,2.54476500804708,0.110660374187379,0.144411939198395
+"gene10",11.5893249023836,0.386011089345352,0.244586104807745,2.47690282625604,0.115529551358716,0.144411939198395
diff --git a/tests/test_pydeseq2.py b/tests/test_pydeseq2.py
index 1a59afbe..881426d3 100644
--- a/tests/test_pydeseq2.py
+++ b/tests/test_pydeseq2.py
@@ -929,6 +929,247 @@ def test_vst_transform_no_fit(train_counts, train_metadata, test_counts):
         train_dds.vst_transform(test_counts.to_numpy())
 
 
+# Likelihood ratio test (LRT)
+
+
+@pytest.mark.parametrize(
+    "design, reduced, ref_file",
+    [
+        ("~condition", "~1", "data/single_factor/r_test_res_lrt.csv"),
+        ("~group + condition", "~group", "data/multi_factor/r_test_res_lrt.csv"),
+        (
+            "~group + condition",
+            "~1",
+            "data/multi_factor/r_test_res_lrt_reduced_intercept.csv",
+        ),
+    ],
+)
+def test_lrt_matches_r(counts_df, metadata, design, reduced, ref_file):
+    """The likelihood ratio test matches R DESeq2's ``DESeq(dds, test="LRT")``.
+
+    Covers a single-factor model (``df=1``), a multi-factor model testing one
+    factor adjusting for another (``df=1``), and a joint test of both factors
+    (``df=2``).
+    """
+    test_path = str(Path(os.path.realpath(tests.__file__)).parent.resolve())
+    r_res = pd.read_csv(os.path.join(test_path, ref_file), index_col=0)
+
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design=design)
+    dds.deseq2()
+
+    ds = DeseqStats(dds, contrast=["condition", "B", "A"], test="LRT", reduced=reduced)
+    ds.summary()
+
+    assert_lrt_res_almost_equal(ds.results_df, r_res)
+
+
+def test_lrt_no_independent_filtering(counts_df, metadata):
+    """The raw LRT statistic and p-value match R with independent filtering off."""
+    test_path = str(Path(os.path.realpath(tests.__file__)).parent.resolve())
+    r_res = pd.read_csv(
+        os.path.join(
+            test_path, "data/multi_factor/r_test_res_lrt_no_independent_filtering.csv"
+        ),
+        index_col=0,
+    )
+
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~group + condition")
+    dds.deseq2()
+
+    ds = DeseqStats(
+        dds,
+        contrast=["condition", "B", "A"],
+        test="LRT",
+        reduced="~group",
+        independent_filter=False,
+    )
+    ds.summary()
+
+    assert_lrt_res_almost_equal(ds.results_df, r_res)
+
+
+def test_lrt_with_outliers(counts_df, metadata, tol=0.04):
+    """The LRT matches R when Cooks outliers are replaced and models refitted.
+
+    Uses the same mild outlier injection as the Wald outlier test. The outlier
+    counts are replaced and both the full and reduced models are refitted on the
+    imputed counts, matching R's ``counts(dds, replaced=TRUE)`` behaviour.
+    """
+    test_path = str(Path(os.path.realpath(tests.__file__)).parent.resolve())
+    r_res = pd.read_csv(
+        os.path.join(test_path, "data/multi_factor/r_test_res_lrt_outliers.csv"),
+        index_col=0,
+    )
+
+    counts_df.loc["sample1", "gene1"] = 2000
+    counts_df.loc["sample11", "gene7"] = 1000
+
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~group + condition")
+    dds.deseq2()
+
+    # At least one gene should have been flagged as an outlier and refitted.
+    assert dds.var["replaced"].sum() > 0
+
+    counts_before = dds.X.copy()
+
+    ds = DeseqStats(dds, contrast=["condition", "B", "A"], test="LRT", reduced="~group")
+    ds.summary()
+
+    # The LRT must not mutate the raw counts when reconstructing replaced counts.
+    np.testing.assert_array_equal(dds.X, counts_before)
+
+    assert_lrt_res_almost_equal(ds.results_df, r_res, tol=tol)
+
+
+def test_lrt_reduced_as_design_matrix(counts_df, metadata):
+    """A reduced model given as a design matrix matches the equivalent formula."""
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~group + condition")
+    dds.deseq2()
+
+    ds_formula = DeseqStats(
+        dds, contrast=["condition", "B", "A"], test="LRT", reduced="~group"
+    )
+    ds_formula.summary()
+
+    # Build the same reduced design as an explicit matrix.
+    reduced_matrix = model_matrix("~group", metadata).to_numpy()
+    ds_matrix = DeseqStats(
+        dds, contrast=["condition", "B", "A"], test="LRT", reduced=reduced_matrix
+    )
+    ds_matrix.summary()
+
+    pd.testing.assert_series_equal(
+        ds_formula.statistics, ds_matrix.statistics, check_names=False
+    )
+    pd.testing.assert_series_equal(
+        ds_formula.p_values, ds_matrix.p_values, check_names=False
+    )
+
+
+def test_lrt_statistic_definition(counts_df, metadata):
+    """The LRT is internally consistent, independent of the R reference.
+
+    Checks that the statistic is non-negative, that the p-value is exactly the
+    chi-squared survival function of the statistic, and that for a single removed
+    coefficient (``df=1``) the LRT statistic is close to the squared Wald
+    statistic (their asymptotic relationship) for genes with signal.
+    """
+    from scipy.stats import chi2
+
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
+    dds.deseq2()
+
+    ds_lrt = DeseqStats(dds, contrast=["condition", "B", "A"], test="LRT", reduced="~1")
+    ds_lrt.summary()
+
+    stat = ds_lrt.statistics.dropna()
+    # Statistic is non-negative.
+    assert (stat >= 0).all()
+    # p-value is the chi2 survival function of the statistic (df = 1 here).
+    df = ds_lrt.design_matrix.shape[1] - ds_lrt.reduced_design_matrix.shape[1]
+    assert df == 1
+    np.testing.assert_allclose(
+        ds_lrt.p_values.dropna().values, chi2.sf(stat.values, df), rtol=1e-10
+    )
+
+    # LRT statistic ~ (Wald statistic)^2 for genes with signal (df = 1).
+    ds_wald = DeseqStats(dds, contrast=["condition", "B", "A"], test="Wald")
+    ds_wald.summary()
+    signal = ds_wald.p_values < 0.1
+    ratio = (stat[signal] / ds_wald.statistics[signal] ** 2).dropna()
+    assert ((ratio > 0.9) & (ratio < 1.1)).all()
+
+
+def test_lrt_lfc_matches_wald(counts_df, metadata):
+    """The LRT reports the same LFC as the Wald test (both use the full model)."""
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~group + condition")
+    dds.deseq2()
+
+    ds_wald = DeseqStats(dds, contrast=["condition", "B", "A"], test="Wald")
+    ds_wald.summary()
+
+    ds_lrt = DeseqStats(
+        dds, contrast=["condition", "B", "A"], test="LRT", reduced="~group"
+    )
+    ds_lrt.summary()
+
+    pd.testing.assert_series_equal(
+        ds_wald.results_df["log2FoldChange"],
+        ds_lrt.results_df["log2FoldChange"],
+        check_names=False,
+    )
+
+
+def test_lrt_input_validation(counts_df, metadata):
+    """Invalid ``test``/``reduced`` combinations raise informative errors."""
+    dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~group + condition")
+    dds.deseq2()
+
+    # LRT requires a reduced model.
+    with pytest.raises(ValueError, match="reduced"):
+        DeseqStats(dds, contrast=["condition", "B", "A"], test="LRT")
+
+    # reduced is not allowed for the Wald test.
+    with pytest.raises(ValueError, match="reduced"):
+        DeseqStats(dds, contrast=["condition", "B", "A"], test="Wald", reduced="~group")
+
+    # Unknown test.
+    with pytest.raises(ValueError, match="test"):
+        DeseqStats(dds, contrast=["condition", "B", "A"], test="chi2")
+
+    # Reduced model not smaller than the full model.
+    with pytest.raises(ValueError, match="fewer coefficients"):
+        DeseqStats(
+            dds,
+            contrast=["condition", "B", "A"],
+            test="LRT",
+            reduced="~group + condition",
+        )
+
+    # Reduced model not nested in the full model (a factor absent from the design).
+    dds_single = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
+    dds_single.deseq2()
+    with pytest.raises(ValueError, match="nested"):
+        DeseqStats(
+            dds_single, contrast=["condition", "B", "A"], test="LRT", reduced="~group"
+        )
+
+
+def assert_lrt_res_almost_equal(py_res, r_res, tol=0.02, stat_atol=0.1, pval_tol=0.03):
+    """Compare PyDESeq2 LRT results to an R DESeq2 reference.
+
+    The full-model log-fold changes match R tightly. The LRT statistic matches
+    with a relative tolerance for genes with signal and an absolute tolerance for
+    near-null genes, whose statistic is ~0 and therefore dominated by IRLS
+    convergence noise (in R as well as here). For genes R deems differentially
+    expressed, p-values and adjusted p-values match R; very small p-values may
+    differ by a few percent because the chi-squared tail amplifies tiny
+    differences in the statistic.
+    """
+    # Same genes flagged NaN (Cooks / all-zero filtering).
+    assert (py_res.pvalue.isna() == r_res.pvalue.isna()).all()
+    assert (py_res.padj.isna() == r_res.padj.isna()).all()
+
+    # Full-model LFCs.
+    assert (
+        abs(r_res.log2FoldChange - py_res.log2FoldChange) / abs(r_res.log2FoldChange)
+    ).max() < tol
+
+    # LRT statistics.
+    common = py_res.stat.dropna().index
+    np.testing.assert_allclose(
+        py_res.stat.loc[common].values,
+        r_res.stat.loc[common].values,
+        rtol=tol,
+        atol=stat_atol,
+    )
+
+    # p-values / adjusted p-values for genes R deems significant.
+    signal = r_res.pvalue < 0.1
+    assert (abs(r_res.pvalue - py_res.pvalue) / r_res.pvalue)[signal].max() < pval_tol
+    assert (abs(r_res.padj - py_res.padj) / r_res.padj)[signal].max() < pval_tol
+
+
 def assert_res_almost_equal(py_res, r_res, tol=0.02):
     # check that the same p-values are NaN
     assert (py_res.pvalue.isna() == r_res.pvalue.isna()).all()