From 69f208b5685c332ffca80ecbc7122e3a1ec182aa Mon Sep 17 00:00:00 2001
From: Nicolas Jouvin <nicolas.jouvin@inrae.fr>
Date: Wed, 29 Oct 2025 11:53:09 +0100
Subject: [PATCH] Clearer error when using singular models

In this commit
 1. Use singular.ok=FALSE in lm.fit(log(Y) ~ X) which causes a traceable R error
 2. Update documentation.
 3. Write tests.
---
 R/utils.R                         |  4 ++--
 man/compute_PLN_starting_point.Rd |  2 +-
 tests/testthat/test-pln.R         | 18 ++++++++++++++++++
 3 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/R/utils.R b/R/utils.R
index a0ae3bc7..c36a44d0 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -269,7 +269,7 @@ create_parameters <- function(
 #' Barebone function to compute starting points for B, M and S when fitting a PLN. Mostly intended for internal use.
 #'
 #' @param Y Response count matrix
-#' @param X Covariate matrix
+#' @param X Covariate matrix. Note that initialization will fail if the model matrix is singular.
 #' @param O Offset matrix (in log-scale)
 #' @param w Weight vector (defaults to 1)
 #' @param s Scale parameter for S (defaults to 0.1)
@@ -293,7 +293,7 @@ create_parameters <- function(
 compute_PLN_starting_point <- function(Y, X, O, w, s = 0.1) {
   # Y = responses, X = covariates, O = offsets (in log scale), w = weights
   n <- nrow(Y); p <- ncol(Y); d <- ncol(X)
-  fits <- lm.fit(w * X, w * log((1 + Y)/exp(O)))
+  fits <- lm.fit(w * X, w * log((1 + Y)/exp(O)), singular.ok = FALSE)
   list(B = matrix(fits$coefficients, d, p),
        M = matrix(fits$residuals, n, p),
        S = matrix(s, n, p))
diff --git a/man/compute_PLN_starting_point.Rd b/man/compute_PLN_starting_point.Rd
index 03fb3ad7..0390059b 100644
--- a/man/compute_PLN_starting_point.Rd
+++ b/man/compute_PLN_starting_point.Rd
@@ -9,7 +9,7 @@ compute_PLN_starting_point(Y, X, O, w, s = 0.1)
 \arguments{
 \item{Y}{Response count matrix}
 
-\item{X}{Covariate matrix}
+\item{X}{Covariate matrix. Note that initialization will fail if the model matrix is singular.}
 
 \item{O}{Offset matrix (in log-scale)}
 
diff --git a/tests/testthat/test-pln.R b/tests/testthat/test-pln.R
index 54f38e5a..a4afbe2a 100644
--- a/tests/testthat/test-pln.R
+++ b/tests/testthat/test-pln.R
@@ -203,3 +203,21 @@ test_that("PLN: Check that all univariate PLN models are equivalent with the mul
   ))
 
 })
+
+test_that("PLN: check initialization fails when the covariate model matrix is singular",  {
+  n = 10; d = 1; p = 10
+  Y <- matrix(rpois(n*p, 1), n, p)
+
+  # Artifically build singular model matrix with continuous covariates
+  x <- matrix(rnorm(n*d), n, d)
+  X_singular <- cbind(x, x)
+
+  # Artificially build correlated factors (discrete covariates)
+  f1 <- gl(2, n/2, labels = c("1.1", "1.2"))
+  f2 <- gl(2, n/2, labels = c("2.1", "2.2"))
+
+  # In both cases, model.matrix(formula) is singular
+  expect_error(PLN(Y ~ X_singular))
+  expect_error(PLN(Y ~ f1 + f2))
+
+})