1

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

Krina M
  • 135
  • 2
  • 13
  • Your formula is identical to `y ~ 1/(1+Epsup/x)^(Bup+Bdn)` Because of that you can only estimate `Bup+Bdn`. Also identical: `y ~ (1+Epsup/x)^(-Bup-Bdn)` – jogo Feb 28 '17 at 14:38
  • Ths is not the same problem as http://stackoverflow.com/questions/34201377/r-and-nls-singular-gradient-matrix-at-initial-parameter. In that question the problem was overparameterization. This problem is not, in general, overparameterized but only if `Epsup = Epsdn` or if any of the parameters are close to zero and as well there is the problem that we must be sure that `1+Epsup/x` and `1+Epsdn/x` remain positive. Thus if we can find starting values sufficiently good to stay away from such bad areas it can be solved.. – G. Grothendieck Mar 02 '17 at 02:01

1 Answers1

2

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
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341