Error in NLME function: Singularity in backsolve at level 0, block 1

Question

I'm trying to fit ELISA plate data to a non-linear mixed effect model, using the nlme function. Previously, I asked about how to make the syntax work for nested groups, though responses made it seem like that was a lost cause. However, I'm still having difficulty producing a model object when just applying a single group.

Sample data can be found here: https://pastebin.com/UzVa04TY

This is my original envisioning of running the NLME function, using separate variables for fixed and random effects (upper- and lower-case letters, respectively). x is the sample dilution factor, while y is ELISA signal density.

# Define the 3-parameter logistic function
logistic3PL <- function(x, A, B, C, a, b, c) {
  (A + a) / (1 + (x / (C + c))^ (B + b))
}

# Define the nonlinear mixed effects model
model <- nlme(y ~ logistic3PL(x, A, B, C, a, b, c),
              data = dataSet,
              fixed = A + B + C ~ 1,
              random = a + b + c ~ 1,
              groups =  ~ Plate,
              start = c(A = 1, B = 1, C = 1))

This produced the following error:

Error in nlme.formula(y ~ logistic3PL(x, A, B, C, a, b, c), data = dataSet,  : 
  Singularity in backsolve at level 0, block 1

A commenter from my prior question suggested that he would actually phrase the 3-parameter function to have both the fixed and random effects denoted by the same set of letters.

# Define the 3-parameter logistic function
logistic3PL <- function(x, A, B, C) {
  A / (1 + (x / C) ^B)
}

# Define the nonlinear mixed effects model
model <- nlme(y ~ logistic3PL(x, A, B, C),
              data = dataSet,
              fixed = A + B + C ~ 1,
              random = A + B + C ~ 1,
              groups =  ~ Plate,
              start = c(A = 1, B = 1, C = 1))

This produced an identical error.

I should note that I also attempted replacing random = A + B + C ~ 1, with ... ~ Plate, in both methods. Each time, the program run would think for about 30 seconds before giving the exact same error, as opposed to giving the error immediately when the line ended with ~ 1.

Happy to answer any questions that may help solve the issue. Thank you in advance for your help!

Answer 1

I suggest using the parameterization given as Equation 1 in the reference [1].

#plot the data
library(ggplot2)
p <- ggplot(dataSet, aes(x, y)) +
  geom_point(aes(color = factor(Plate))) +
  scale_x_continuous(trans = "log2")
print(p)

#we use the other parametrization from the reference
foo <- function(x, A, B, C) A / (1 + exp((C - log2(x))/B))

#or with gradient attribute as suggested by Ben Bolker in the comments
foo <- deriv(~ A / (1 + exp((C - log2(x))/B)),
         namevec = c("A", "B", "C"),
         function.arg = c("x", "A", "B", "C"))

#plot the function to see if the starting values are sensible 
p + geom_function(fun = foo, args = list(A = 1, B = 1, C = 6))

library(nlme)
#first use nlsList
fit1 <- nlsList(y ~ foo(x, A, B, C) | Plate, data = dataSet, 
                start = list(A = 1, B = 1, C = 6))

#then use that fit as input for nlme
#specify no correlation between random effects with pdDiag
fit2 <- nlme(fit1, random = pdDiag(list(A ~ 1, B ~ 1, C ~ 1)), 
             groups = ~ Plate)
summary(fit2)

#plot predictions
p + lapply(1:3, 
           \(i) geom_function(fun = \(x) predict(fit2, newdata = data.frame(x = x, Plate = i)),
                              aes(color = factor(i))))

Most people would agree that three subjects are not enough to estimate a random effect. You should dramatically increased the number of plates. Then maybe you could even estimate random effects with an unstructured covariance matrix (like you attempted in your question), i.e., including correlations.

[1]: https://www.aphis.usda.gov/animal_health/vet_biologics/publications/STATWI0006.pdf

---

Data:

dataSet <- structure(list(Plate = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3
), x = c(5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 5, 
         10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120), y = c(0.0128, 0.0269, 0.0611, 0.1783, 0.4333, 0.6237, 0.846, 0.8952, 0.9145, 
        0.9443, 0.9831, 0.0332, 0.042, 0.0755, 0.172, 0.4082, 0.6905, 0.8783, 0.9687, 1.0235, 1.0798, 1.0744, 0.0139, 0.0307, 0.0592, 
        0.1592, 0.4117, 0.641, 0.8177, 0.8915, 0.9636, 1.0324, 0.9909)), row.names = c(NA, 33L), class = "data.frame")