Search code examples
rcurve-fittingnls

NLS Curve Fit Singular matrix error

I need help fitting nls and finding initial estimates that would not lead to singular matrix. I will greatly appreciate any help.

via_data$Concentration <- c(0.197, 0.398, 0.792, 1.575,
                            3.154, 6.270, 12.625, 25.277,
                            25.110, 49.945, 74.680)
via_data$Viability <- c(100, 94.62, 96.21, 87.53, 
                        80,  62.22,  39.11, 
                        30.80,  30, 22, 2.56) 
x <- via_data$Concentration
y <- via_data$Viability
fit <- nls(y ~((1/(1+Epsup/x)^Bup)*(1/(1+Epsdn/x)^Bdn)), start=list(Epsup=0, Bup=1, Epsdn=10, Bdn=-5), trace=T)

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

thanks, Krina

Solution

  • Below st is your starting values except we have used Epsup=1 to avoid degeneracy at 0. fo is the formula. To prevent raising negative numbers to a power we have replaced Epsup with sqrt(Epsup^2) and similarly for Epsdn -- this adds the assumption that Epsup and Espdn cannot be negative. (This would be the same as using abs(Epsup); however, nlxb does not have abs in its derivative table.) Next use nls2 to generate values on a grid between the boundaries st/10 and 10*st. nls2 will generate these and return an "nls" object with the best one found. Now use that as the starting values for nlxb of the nlmrt package. It handles difficult problems bettter than nls. nlxb does not return an "nls" object (although the package does have wrapnls which runs nlxb followed by nls but then we don't get the direct output from nlxb) so feed that through nls2 again to create an "nls" object allowing us to use the fitted method. We plot the resulting fit.

    library(nlmrt)
library(nls2)

st <- c(Epsup=1, Bup=1, Epsdn=10, Bdn=-5)
fo <- y ~ (1/(1+sqrt(Epsup^2)/x)^Bup)*(1/(1+sqrt(Epsdn^2)/x)^Bdn)

fit.nls2 <- nls2(fo, start = data.frame(rbind(st/10, 10*st)), alg = "brute")
fit.nlxb <- nlxb(fo, data = data.frame(x, y), start = coef(fit.nls2))

    giving the following:

    > fit.nlxb
nlmrt class object: x 
residual sumsquares =  171.2  on  11 observations
    after  19    Jacobian and  25 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval   
Epsup            10.7464         10.95     0.9814     0.3591   6.855e-05        1584  
Bup              1.15049        0.5928      1.941    0.09345    0.001839       120.2  
Epsdn            642.754         908.5     0.7075     0.5021  -1.298e-06       1.406  
Bdn             -1.13885        0.6315     -1.804     0.1143    0.004964    0.005443

    and plotting to visually assess the fit:

    fit.nlxb.nls <- nls2(fo, start = coef(fit.nlxb))
plot(y ~ x)
lines(fitted(fit.nlxb.nls) ~ x)

    screenshot

    Note: We used this input:

    via_data <- data.frame(Concentration = c(0.197, 0.398, 0.792, 1.575,
     3.154, 6.270, 12.625, 25.277, 25.110, 49.945, 74.680),
Viability = c(100, 94.62, 96.21, 87.53, 80,  62.22,  39.11, 
                        30.80,  30, 22, 2.56))
x <- via_data$Concentration
y <- via_data$Viability