Clearer error when using singular models

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))

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))
+
+})